Previous Next Contents Generated Index Doc Set Home



Solution to a Linear System in a UDU- or LDL-Factored Hermitian Matrix in Packed Storage

The subroutines described in this section solve the linear system AX = B for a Hermitian matrix A in packed storage, which has been UDU-factored or LDL-factored by xHPTRF, and general matrices B and X.

Calling Sequence

CALL ZHPTRS 
(UPLO, N, NRHS, ZA, IPIVOT, ZB, LDB, INFO)
CALL CHPTRS 
(UPLO, N, NRHS, CA, IPIVOT, CB, LDB, INFO)






void zhptrs 
(char uplo, int n, int nrhs, doublecomplex *za, int 
*ipivot, doublecomplex *zb, int ldb, int *info)
void chptrs 
(char uplo, int n, int nrhs, complex *ca, int *ipivot, 
complex *cb, int ldb, int *info)

Arguments

UPLO

Indicates whether xA contains the upper or lower triangle of the matrix. The legal values for UPLO are listed below. Any values not listed below are illegal.

'U' or 'u'

xA contains the upper triangle.

'L' or 'l'

xA contains the lower triangle.

N

Order of the matrix A. N 0.

NRHS

Number of right-hand sides, equal to the number of columns of the matrix B. NRHS 0.

xA

UDU or LDL factorization of the matrix A, as computed by xHPTRF. The dimension of xA is (N × N + N) / 2.

IPIVOT

Pivot indices as computed by xHPTRF.

xB

On entry, the N×NRHS right-hand side matrix B.
On exit, the N×NRHS solution matrix X.

LDB

Leading dimension of the array B as specified in a dimension or type statement. LDB max(1, N).

INFO

On exit:

INFO = 0

Subroutine completed normally.

INFO < 0

The ith argument, where i = |INFO|, had an illegal value.

Sample Program




      PROGRAM TEST
      IMPLICIT NONE
C
      DOUBLE PRECISION  ZERO
      INTEGER           LENGTA, LDB, LDWORK, N, NRHS
      PARAMETER        (N = 3)
      PARAMETER        (LENGTA = (N * N + N) / 2)
      PARAMETER        (LDB = N)
      PARAMETER        (LDWORK = 2 * N)
      PARAMETER        (NRHS = 1)
      PARAMETER        (ZERO = 0.0D0)
C
      DOUBLE PRECISION  ANORM, RCOND
      COMPLEX*16        A(LENGTA), B(LDB,NRHS), WORK(LDWORK)
      INTEGER           ICOL, INFO, IPIVOT(N), IROW
C
      EXTERNAL          ZHPCON, ZHPTRF, ZHPTRS
      INTRINSIC         ABS, CMPLX, CONJG, DBLE
C
C     Initialize the array A to store the coefficient array A
C     shown below.  Initialize the array B to store the right
C     hand side vector b shown below.
C
C         1    1-i  1-i        4+i
C     A = 1+i   3   3-i    b = 4+6i
C         1+i  3+i   5         4+7i
C
      DATA A / (1.0D0,8D8), (1.0D0,-1.0D0), (3.0D0,8D8),
     $         (1.0D0,-1.0D0), (3.0D0,-1.0D0), (5.0D0,8D8) /
      DATA B / (4.0D0,1.0D0), (4.0D0,6.0D0), (4.0D0,7.0D0) /
C
C     Print the initial value of the arrays.
C
      PRINT 1000
      PRINT 1010, CMPLX(DBLE(A(1)),ZERO), A(2), A(4)
      PRINT 1010, CONJG(A(2)), CMPLX(DBLE(A(3)),ZERO), A(5)
      PRINT 1010, CONJG(A(4)), CONJG(A(5)), CMPLX(DBLE(A(6)),
     $   ZERO)
      PRINT 1020
      DO 130, ICOL = 1, NRHS
        DO 120, IROW = 1, N
          PRINT 1010, B(IROW,ICOL)
  120   CONTINUE
  130 CONTINUE
C
C     Factor the matrix.
C
      CALL ZHPTRF ('UPPER TRIANGLE OF A STORED', N, A, IPIVOT,
     $             INFO)
      IF (INFO .LT. 0) THEN
        PRINT 1030, ABS(INFO)
        STOP 1
      ELSE IF (INFO .GT. 0) THEN
        PRINT 1040
        STOP 2
      END IF
C
C     Estimate the condition number.  Print a warning if it is
C     unreasonably high.
C
      ANORM = 1.0D1
      CALL ZHPCON ('UPPER TRIANGULAR FACTOR OF A', N, A, IPIVOT,
     $             ANORM, RCOND, WORK, INFO)
      IF (INFO .LT. 0) THEN
        PRINT 1050, ABS(INFO)
        STOP 3
      END IF
      PRINT 1060, RCOND
      IF ((1.0D0 + RCOND) .EQ. 1.0D0) THEN
        PRINT 1070
      END IF
C
C     Compute and print the solution.
C
      CALL ZHPTRS ('UPPER TRIANGULAR FACTOR OF A', N, NRHS, A,
     $             IPIVOT, B, LDB, INFO)
      IF (INFO .LT. 0) THEN
        PRINT 1080, ABS(INFO)
        STOP 4
      END IF
      PRINT 1090
      DO 210, ICOL = 1, NRHS
        DO 200, IROW = 1, N
          PRINT 1010, B(IROW,ICOL)
  200   CONTINUE
  210 CONTINUE
C
 1000 FORMAT (1X, 'A:')
 1010 FORMAT (1X, 10(: 2X, '(', F4.1, ',', F4.1, ')'))
 1020 FORMAT (/1X, 'b:')
 1030 FORMAT (1X, 'Illegal argument to ZHPTRF, argument #', I2)
 1040 FORMAT (1X, 'A is singular to working precision.')
 1050 FORMAT (1X, 'Illegal argument to ZHPCON, argument #', I2)
 1060 FORMAT (/1X,
     $        'Estimated reciprocal condition number of A: ',
     $        F5.3)
 1070 FORMAT (1X, 'A may be singular to working precision.')
 1080 FORMAT (1X, 'Illegal argument to ZHPTRS, argument #', I2)
 1090 FORMAT (/1X, 'x:')
C
      END
 

Sample Output

 
 A:
   ( 1.0, 0.0)  ( 1.0,-1.0)  ( 1.0,-1.0)
   ( 1.0, 1.0)  ( 3.0, 0.0)  ( 3.0,-1.0)
   ( 1.0, 1.0)  ( 3.0, 1.0)  ( 5.0, 0.0)



 b:
   ( 4.0, 1.0)
   ( 4.0, 6.0)
   ( 4.0, 7.0)



 Estimated reciprocal condition number of A: 0.010



 x:
   ( 1.0, 0.0)
   ( 0.0, 2.0)
   ( 1.0, 0.0)






Previous Next Contents Generated Index Doc Set Home