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