![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, 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 DLACN2 in place of DLACON, 5 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: INTEGER IWORK( * ) 18: DOUBLE PRECISION A( LDA, * ), WORK( * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * DGECON estimates the reciprocal of the condition number of a general 25: * real matrix A, in either the 1-norm or the infinity-norm, using 26: * the LU factorization computed by DGETRF. 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) DOUBLE PRECISION array, dimension (LDA,N) 45: * The factors L and U from the factorization A = P*L*U 46: * as computed by DGETRF. 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) DOUBLE PRECISION array, dimension (4*N) 60: * 61: * IWORK (workspace) INTEGER array, dimension (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: * .. 79: * .. Local Arrays .. 80: INTEGER ISAVE( 3 ) 81: * .. 82: * .. External Functions .. 83: LOGICAL LSAME 84: INTEGER IDAMAX 85: DOUBLE PRECISION DLAMCH 86: EXTERNAL LSAME, IDAMAX, DLAMCH 87: * .. 88: * .. External Subroutines .. 89: EXTERNAL DLACN2, DLATRS, DRSCL, XERBLA 90: * .. 91: * .. Intrinsic Functions .. 92: INTRINSIC ABS, MAX 93: * .. 94: * .. Executable Statements .. 95: * 96: * Test the input parameters. 97: * 98: INFO = 0 99: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' ) 100: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN 101: INFO = -1 102: ELSE IF( N.LT.0 ) THEN 103: INFO = -2 104: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 105: INFO = -4 106: ELSE IF( ANORM.LT.ZERO ) THEN 107: INFO = -5 108: END IF 109: IF( INFO.NE.0 ) THEN 110: CALL XERBLA( 'DGECON', -INFO ) 111: RETURN 112: END IF 113: * 114: * Quick return if possible 115: * 116: RCOND = ZERO 117: IF( N.EQ.0 ) THEN 118: RCOND = ONE 119: RETURN 120: ELSE IF( ANORM.EQ.ZERO ) THEN 121: RETURN 122: END IF 123: * 124: SMLNUM = DLAMCH( 'Safe minimum' ) 125: * 126: * Estimate the norm of inv(A). 127: * 128: AINVNM = ZERO 129: NORMIN = 'N' 130: IF( ONENRM ) THEN 131: KASE1 = 1 132: ELSE 133: KASE1 = 2 134: END IF 135: KASE = 0 136: 10 CONTINUE 137: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 138: IF( KASE.NE.0 ) THEN 139: IF( KASE.EQ.KASE1 ) THEN 140: * 141: * Multiply by inv(L). 142: * 143: CALL DLATRS( 'Lower', 'No transpose', 'Unit', NORMIN, N, A, 144: $ LDA, WORK, SL, WORK( 2*N+1 ), INFO ) 145: * 146: * Multiply by inv(U). 147: * 148: CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 149: $ A, LDA, WORK, SU, WORK( 3*N+1 ), INFO ) 150: ELSE 151: * 152: * Multiply by inv(U'). 153: * 154: CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A, 155: $ LDA, WORK, SU, WORK( 3*N+1 ), INFO ) 156: * 157: * Multiply by inv(L'). 158: * 159: CALL DLATRS( 'Lower', 'Transpose', 'Unit', NORMIN, N, A, 160: $ LDA, WORK, SL, WORK( 2*N+1 ), INFO ) 161: END IF 162: * 163: * Divide X by 1/(SL*SU) if doing so will not cause overflow. 164: * 165: SCALE = SL*SU 166: NORMIN = 'Y' 167: IF( SCALE.NE.ONE ) THEN 168: IX = IDAMAX( N, WORK, 1 ) 169: IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 170: $ GO TO 20 171: CALL DRSCL( N, SCALE, WORK, 1 ) 172: END IF 173: GO TO 10 174: END IF 175: * 176: * Compute the estimate of the reciprocal condition number. 177: * 178: IF( AINVNM.NE.ZERO ) 179: $ RCOND = ( ONE / AINVNM ) / ANORM 180: * 181: 20 CONTINUE 182: RETURN 183: * 184: * End of DGECON 185: * 186: END