Previous Next Contents Generated Index Doc Set Home


UDU Factorization and Condition Number of a Hermitian Matrix in Packed Storage

The subroutines described in this section compute the UDU factorization and condition number of a Hermitian matrix A in packed storage. If the condition number is not needed then xHPFA is slightly faster. It is typical to follow a call to xHPCO with a call to xHPSL to solve Ax = b or to xHPDI to compute the determinant, inverse, and inertia of A.

Calling Sequence

CALL ZHPCO

(ZA, N, IPIVOT, DRCOND, ZWORK)

CALL CHPCO

(CA, N, IPIVOT, SRCOND, CWORK)

void zhpco

(doublecomplex *za, int n, int *ipivot, double *drcond)

void chpco

(complex *ca, int n, int *ipivot, float *srcond)

Arguments

xA

On entry, the upper triangle of the matrix A.
On exit, a UDU factorization of the matrix A.

N

Order of the matrix A. N 0.

IPIVOT

On exit, a vector of pivot indices.

xRCOND

On exit, an estimate of the reciprocal condition number of A.
0.0 RCOND 1.0. As the value of RCOND gets smaller, operations with A such as solving Ax = b may become less stable. If RCOND satisfies RCOND + 1.0 = 1.0 then A may be singular to working precision.

xWORK

Scratch array with a dimension of N.

Sample Program

 

      PROGRAM TEST

      IMPLICIT NONE

C

      INTEGER    LENGTA, N

      PARAMETER (N = 3)

      PARAMETER (LENGTA = (N * N + N) / 2)

C

      REAL       RCOND

      COMPLEX    A(LENGTA), B(N), WORK(N)

      INTEGER    IPIVOT(N)

C

      EXTERNAL   CHPCO, CHPSL

      INTRINSIC  CONJG

C

C     Initialize the array A to store the matrix A shown below.

C     Initialize the array B to store the vector b shown below.

C

C          1    1+2i  1+2i         95-180i

C     A = 1+2i   6   -2+6i    b = 545-118i

C         1+2i -2+6i   11         865+ 62i

C

      DATA A / (1.0,0.0), (1.0,-2.0), (6.0,0.0),

     $         (1.0,-2.0), (6.0,-2.0), (11.0,0.0) /

      DATA B / (95.0,-180.0), (545.0,-118.0), (865.0,62.0) /

C

      PRINT 1000

      PRINT 1010, A(1), A(2), A(4)

      PRINT 1010, CONJG(A(2)), A(3), A(5)

      PRINT 1010, CONJG(A(4)), CONJG(A(5)), A(6)

      PRINT 1020

      PRINT 1030, B

      CALL CHPCO (A, N, IPIVOT, RCOND, WORK)

      PRINT 1040, RCOND

      IF ((RCOND + 1.0) .EQ. 1.0) THEN

        PRINT 1050

      END IF

      CALL CHPSL (A, N, IPIVOT, B)

      PRINT 1060

      PRINT 1030, B

C

 1000 FORMAT (1X, 'A in full form:')

 1010 FORMAT (3(3X, '(', F4.1, ',', F4.1, ')'))

 1020 FORMAT (/1X, 'b:')

 1030 FORMAT (3X, '(', F6.1, ',', F6.1, ')')

 1040 FORMAT (/1X, 'Reciprocal condition number of A:', F6.3)

 1050 FORMAT (1X, 'A may be singular to working precision.')

 1060 FORMAT (/1X, 'A**(-1) * b:')

C

      END

Sample Output

 

 A in full form:

   ( 1.0, 0.0)   ( 1.0,-2.0)   ( 1.0,-2.0)

   ( 1.0, 2.0)   ( 6.0, 0.0)   ( 6.0,-2.0)

   ( 1.0, 2.0)   ( 6.0, 2.0)   (11.0, 0.0)










 b:










   (  95.0,-180.0)










   ( 545.0,-118.0)










   ( 865.0,  62.0)



















 Reciprocal condition number of A: 0.001



















 A**(-1) * b:










   (   5.0,   0.0)










   (  26.0,   0.0)










   (  64.0,   0.0)



















 























Previous


Next

Contents

Generated Index

Doc Set

Home