1 SUBROUTINE ZPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, 2 $ INFO ) 3* 4* -- LAPACK routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* March 31, 1993 8* 9* .. Scalar Arguments .. 10 CHARACTER UPLO 11 INTEGER INFO, LDA, N 12 DOUBLE PRECISION ANORM, RCOND 13* .. 14* .. Array Arguments .. 15 DOUBLE PRECISION RWORK( * ) 16 COMPLEX*16 A( LDA, * ), WORK( * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* ZPOCON estimates the reciprocal of the condition number (in the 23* 1-norm) of a complex Hermitian positive definite matrix using the 24* Cholesky factorization A = U**H*U or A = L*L**H computed by ZPOTRF. 25* 26* An estimate is obtained for norm(inv(A)), and the reciprocal of the 27* condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). 28* 29* Arguments 30* ========= 31* 32* UPLO (input) CHARACTER*1 33* = 'U': Upper triangle of A is stored; 34* = 'L': Lower triangle of A is stored. 35* 36* N (input) INTEGER 37* The order of the matrix A. N >= 0. 38* 39* A (input) COMPLEX*16 array, dimension (LDA,N) 40* The triangular factor U or L from the Cholesky factorization 41* A = U**H*U or A = L*L**H, as computed by ZPOTRF. 42* 43* LDA (input) INTEGER 44* The leading dimension of the array A. LDA >= max(1,N). 45* 46* ANORM (input) DOUBLE PRECISION 47* The 1-norm (or infinity-norm) of the Hermitian matrix A. 48* 49* RCOND (output) DOUBLE PRECISION 50* The reciprocal of the condition number of the matrix A, 51* computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an 52* estimate of the 1-norm of inv(A) computed in this routine. 53* 54* WORK (workspace) COMPLEX*16 array, dimension (2*N) 55* 56* RWORK (workspace) DOUBLE PRECISION array, dimension (N) 57* 58* INFO (output) INTEGER 59* = 0: successful exit 60* < 0: if INFO = -i, the i-th argument had an illegal value 61* 62* ===================================================================== 63* 64* .. Parameters .. 65 DOUBLE PRECISION ONE, ZERO 66 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 67* .. 68* .. Local Scalars .. 69 LOGICAL UPPER 70 CHARACTER NORMIN 71 INTEGER IX, KASE 72 DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM 73 COMPLEX*16 ZDUM 74* .. 75* .. External Functions .. 76 LOGICAL LSAME 77 INTEGER IZAMAX 78 DOUBLE PRECISION DLAMCH 79 EXTERNAL LSAME, IZAMAX, DLAMCH 80* .. 81* .. External Subroutines .. 82 EXTERNAL XERBLA, ZDRSCL, ZLACON, ZLATRS 83* .. 84* .. Intrinsic Functions .. 85 INTRINSIC ABS, DBLE, DIMAG, MAX 86* .. 87* .. Statement Functions .. 88 DOUBLE PRECISION CABS1 89* .. 90* .. Statement Function definitions .. 91 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 92* .. 93* .. Executable Statements .. 94* 95* Test the input parameters. 96* 97 INFO = 0 98 UPPER = LSAME( UPLO, 'U' ) 99 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 100 INFO = -1 101 ELSE IF( N.LT.0 ) THEN 102 INFO = -2 103 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 104 INFO = -4 105 ELSE IF( ANORM.LT.ZERO ) THEN 106 INFO = -5 107 END IF 108 IF( INFO.NE.0 ) THEN 109 CALL XERBLA( 'ZPOCON', -INFO ) 110 RETURN 111 END IF 112* 113* Quick return if possible 114* 115 RCOND = ZERO 116 IF( N.EQ.0 ) THEN 117 RCOND = ONE 118 RETURN 119 ELSE IF( ANORM.EQ.ZERO ) THEN 120 RETURN 121 END IF 122* 123 SMLNUM = DLAMCH( 'Safe minimum' ) 124* 125* Estimate the 1-norm of inv(A). 126* 127 KASE = 0 128 NORMIN = 'N' 129 10 CONTINUE 130 CALL ZLACON( N, WORK( N+1 ), WORK, AINVNM, KASE ) 131 IF( KASE.NE.0 ) THEN 132 IF( UPPER ) THEN 133* 134* Multiply by inv(U'). 135* 136 CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', 137 $ NORMIN, N, A, LDA, WORK, SCALEL, RWORK, INFO ) 138 NORMIN = 'Y' 139* 140* Multiply by inv(U). 141* 142 CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 143 $ A, LDA, WORK, SCALEU, RWORK, INFO ) 144 ELSE 145* 146* Multiply by inv(L). 147* 148 CALL ZLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, 149 $ A, LDA, WORK, SCALEL, RWORK, INFO ) 150 NORMIN = 'Y' 151* 152* Multiply by inv(L'). 153* 154 CALL ZLATRS( 'Lower', 'Conjugate transpose', 'Non-unit', 155 $ NORMIN, N, A, LDA, WORK, SCALEU, RWORK, INFO ) 156 END IF 157* 158* Multiply by 1/SCALE if doing so will not cause overflow. 159* 160 SCALE = SCALEL*SCALEU 161 IF( SCALE.NE.ONE ) THEN 162 IX = IZAMAX( N, WORK, 1 ) 163 IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 164 $ GO TO 20 165 CALL ZDRSCL( N, SCALE, WORK, 1 ) 166 END IF 167 GO TO 10 168 END IF 169* 170* Compute the estimate of the reciprocal condition number. 171* 172 IF( AINVNM.NE.ZERO ) 173 $ RCOND = ( ONE / AINVNM ) / ANORM 174* 175 20 CONTINUE 176 RETURN 177* 178* End of ZPOCON 179* 180 END 181