![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZGECON( NORM, 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 NORM 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: * ZGECON estimates the reciprocal of the condition number of a general 25: * complex matrix A, in either the 1-norm or the infinity-norm, using 26: * the LU factorization computed by ZGETRF. 27: * 28: * An estimate is obtained for norm(inv(A)), and the reciprocal of the 29: * condition number is computed as 30: * RCOND = 1 / ( norm(A) * norm(inv(A)) ). 31: * 32: * Arguments 33: * ========= 34: * 35: * NORM (input) CHARACTER*1 36: * Specifies whether the 1-norm condition number or the 37: * infinity-norm condition number is required: 38: * = '1' or 'O': 1-norm; 39: * = 'I': Infinity-norm. 40: * 41: * N (input) INTEGER 42: * The order of the matrix A. N >= 0. 43: * 44: * A (input) COMPLEX*16 array, dimension (LDA,N) 45: * The factors L and U from the factorization A = P*L*U 46: * as computed by ZGETRF. 47: * 48: * LDA (input) INTEGER 49: * The leading dimension of the array A. LDA >= max(1,N). 50: * 51: * ANORM (input) DOUBLE PRECISION 52: * If NORM = '1' or 'O', the 1-norm of the original matrix A. 53: * If NORM = 'I', the infinity-norm of the original matrix A. 54: * 55: * RCOND (output) DOUBLE PRECISION 56: * The reciprocal of the condition number of the matrix A, 57: * computed as RCOND = 1/(norm(A) * norm(inv(A))). 58: * 59: * WORK (workspace) COMPLEX*16 array, dimension (2*N) 60: * 61: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) 62: * 63: * INFO (output) INTEGER 64: * = 0: successful exit 65: * < 0: if INFO = -i, the i-th argument had an illegal value 66: * 67: * ===================================================================== 68: * 69: * .. Parameters .. 70: DOUBLE PRECISION ONE, ZERO 71: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 72: * .. 73: * .. Local Scalars .. 74: LOGICAL ONENRM 75: CHARACTER NORMIN 76: INTEGER IX, KASE, KASE1 77: DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU 78: COMPLEX*16 ZDUM 79: * .. 80: * .. Local Arrays .. 81: INTEGER ISAVE( 3 ) 82: * .. 83: * .. External Functions .. 84: LOGICAL LSAME 85: INTEGER IZAMAX 86: DOUBLE PRECISION DLAMCH 87: EXTERNAL LSAME, IZAMAX, DLAMCH 88: * .. 89: * .. External Subroutines .. 90: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLATRS 91: * .. 92: * .. Intrinsic Functions .. 93: INTRINSIC ABS, DBLE, DIMAG, MAX 94: * .. 95: * .. Statement Functions .. 96: DOUBLE PRECISION CABS1 97: * .. 98: * .. Statement Function definitions .. 99: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 100: * .. 101: * .. Executable Statements .. 102: * 103: * Test the input parameters. 104: * 105: INFO = 0 106: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' ) 107: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN 108: INFO = -1 109: ELSE IF( N.LT.0 ) THEN 110: INFO = -2 111: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 112: INFO = -4 113: ELSE IF( ANORM.LT.ZERO ) THEN 114: INFO = -5 115: END IF 116: IF( INFO.NE.0 ) THEN 117: CALL XERBLA( 'ZGECON', -INFO ) 118: RETURN 119: END IF 120: * 121: * Quick return if possible 122: * 123: RCOND = ZERO 124: IF( N.EQ.0 ) THEN 125: RCOND = ONE 126: RETURN 127: ELSE IF( ANORM.EQ.ZERO ) THEN 128: RETURN 129: END IF 130: * 131: SMLNUM = DLAMCH( 'Safe minimum' ) 132: * 133: * Estimate the norm of inv(A). 134: * 135: AINVNM = ZERO 136: NORMIN = 'N' 137: IF( ONENRM ) THEN 138: KASE1 = 1 139: ELSE 140: KASE1 = 2 141: END IF 142: KASE = 0 143: 10 CONTINUE 144: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 145: IF( KASE.NE.0 ) THEN 146: IF( KASE.EQ.KASE1 ) THEN 147: * 148: * Multiply by inv(L). 149: * 150: CALL ZLATRS( 'Lower', 'No transpose', 'Unit', NORMIN, N, A, 151: $ LDA, WORK, SL, RWORK, INFO ) 152: * 153: * Multiply by inv(U). 154: * 155: CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 156: $ A, LDA, WORK, SU, RWORK( N+1 ), INFO ) 157: ELSE 158: * 159: * Multiply by inv(U'). 160: * 161: CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', 162: $ NORMIN, N, A, LDA, WORK, SU, RWORK( N+1 ), 163: $ INFO ) 164: * 165: * Multiply by inv(L'). 166: * 167: CALL ZLATRS( 'Lower', 'Conjugate transpose', 'Unit', NORMIN, 168: $ N, A, LDA, WORK, SL, RWORK, INFO ) 169: END IF 170: * 171: * Divide X by 1/(SL*SU) if doing so will not cause overflow. 172: * 173: SCALE = SL*SU 174: NORMIN = 'Y' 175: IF( SCALE.NE.ONE ) THEN 176: IX = IZAMAX( N, WORK, 1 ) 177: IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 178: $ GO TO 20 179: CALL ZDRSCL( N, SCALE, WORK, 1 ) 180: END IF 181: GO TO 10 182: END IF 183: * 184: * Compute the estimate of the reciprocal condition number. 185: * 186: IF( AINVNM.NE.ZERO ) 187: $ RCOND = ( ONE / AINVNM ) / ANORM 188: * 189: 20 CONTINUE 190: RETURN 191: * 192: * End of ZGECON 193: * 194: END