001:       SUBROUTINE ZHPCON( UPLO, N, AP, IPIV, ANORM, RCOND, WORK, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          UPLO
011:       INTEGER            INFO, N
012:       DOUBLE PRECISION   ANORM, RCOND
013: *     ..
014: *     .. Array Arguments ..
015:       INTEGER            IPIV( * )
016:       COMPLEX*16         AP( * ), WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  ZHPCON estimates the reciprocal of the condition number of a complex
023: *  Hermitian packed matrix A using the factorization A = U*D*U**H or
024: *  A = L*D*L**H computed by ZHPTRF.
025: *
026: *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
027: *  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
028: *
029: *  Arguments
030: *  =========
031: *
032: *  UPLO    (input) CHARACTER*1
033: *          Specifies whether the details of the factorization are stored
034: *          as an upper or lower triangular matrix.
035: *          = 'U':  Upper triangular, form is A = U*D*U**H;
036: *          = 'L':  Lower triangular, form is A = L*D*L**H.
037: *
038: *  N       (input) INTEGER
039: *          The order of the matrix A.  N >= 0.
040: *
041: *  AP      (input) COMPLEX*16 array, dimension (N*(N+1)/2)
042: *          The block diagonal matrix D and the multipliers used to
043: *          obtain the factor U or L as computed by ZHPTRF, stored as a
044: *          packed triangular matrix.
045: *
046: *  IPIV    (input) INTEGER array, dimension (N)
047: *          Details of the interchanges and the block structure of D
048: *          as determined by ZHPTRF.
049: *
050: *  ANORM   (input) DOUBLE PRECISION
051: *          The 1-norm of the original matrix A.
052: *
053: *  RCOND   (output) DOUBLE PRECISION
054: *          The reciprocal of the condition number of the matrix A,
055: *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
056: *          estimate of the 1-norm of inv(A) computed in this routine.
057: *
058: *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
059: *
060: *  INFO    (output) INTEGER
061: *          = 0:  successful exit
062: *          < 0:  if INFO = -i, the i-th argument had an illegal value
063: *
064: *  =====================================================================
065: *
066: *     .. Parameters ..
067:       DOUBLE PRECISION   ONE, ZERO
068:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
069: *     ..
070: *     .. Local Scalars ..
071:       LOGICAL            UPPER
072:       INTEGER            I, IP, KASE
073:       DOUBLE PRECISION   AINVNM
074: *     ..
075: *     .. Local Arrays ..
076:       INTEGER            ISAVE( 3 )
077: *     ..
078: *     .. External Functions ..
079:       LOGICAL            LSAME
080:       EXTERNAL           LSAME
081: *     ..
082: *     .. External Subroutines ..
083:       EXTERNAL           XERBLA, ZHPTRS, ZLACN2
084: *     ..
085: *     .. Executable Statements ..
086: *
087: *     Test the input parameters.
088: *
089:       INFO = 0
090:       UPPER = LSAME( UPLO, 'U' )
091:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
092:          INFO = -1
093:       ELSE IF( N.LT.0 ) THEN
094:          INFO = -2
095:       ELSE IF( ANORM.LT.ZERO ) THEN
096:          INFO = -5
097:       END IF
098:       IF( INFO.NE.0 ) THEN
099:          CALL XERBLA( 'ZHPCON', -INFO )
100:          RETURN
101:       END IF
102: *
103: *     Quick return if possible
104: *
105:       RCOND = ZERO
106:       IF( N.EQ.0 ) THEN
107:          RCOND = ONE
108:          RETURN
109:       ELSE IF( ANORM.LE.ZERO ) THEN
110:          RETURN
111:       END IF
112: *
113: *     Check that the diagonal matrix D is nonsingular.
114: *
115:       IF( UPPER ) THEN
116: *
117: *        Upper triangular storage: examine D from bottom to top
118: *
119:          IP = N*( N+1 ) / 2
120:          DO 10 I = N, 1, -1
121:             IF( IPIV( I ).GT.0 .AND. AP( IP ).EQ.ZERO )
122:      $         RETURN
123:             IP = IP - I
124:    10    CONTINUE
125:       ELSE
126: *
127: *        Lower triangular storage: examine D from top to bottom.
128: *
129:          IP = 1
130:          DO 20 I = 1, N
131:             IF( IPIV( I ).GT.0 .AND. AP( IP ).EQ.ZERO )
132:      $         RETURN
133:             IP = IP + N - I + 1
134:    20    CONTINUE
135:       END IF
136: *
137: *     Estimate the 1-norm of the inverse.
138: *
139:       KASE = 0
140:    30 CONTINUE
141:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
142:       IF( KASE.NE.0 ) THEN
143: *
144: *        Multiply by inv(L*D*L') or inv(U*D*U').
145: *
146:          CALL ZHPTRS( UPLO, N, 1, AP, IPIV, WORK, N, INFO )
147:          GO TO 30
148:       END IF
149: *
150: *     Compute the estimate of the reciprocal condition number.
151: *
152:       IF( AINVNM.NE.ZERO )
153:      $   RCOND = ( ONE / AINVNM ) / ANORM
154: *
155:       RETURN
156: *
157: *     End of ZHPCON
158: *
159:       END
160: