![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: DOUBLE PRECISION FUNCTION DLA_GERCOND ( TRANS, N, A, LDA, AF, 2: $ LDAF, IPIV, CMODE, C, 3: $ INFO, WORK, IWORK ) 4: * 5: * -- LAPACK routine (version 3.2.1) -- 6: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- 7: * -- Jason Riedy of Univ. of California Berkeley. -- 8: * -- April 2009 -- 9: * 10: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 11: * -- Univ. of California Berkeley and NAG Ltd. -- 12: * 13: IMPLICIT NONE 14: * .. 15: * .. Scalar Arguments .. 16: CHARACTER TRANS 17: INTEGER N, LDA, LDAF, INFO, CMODE 18: * .. 19: * .. Array Arguments .. 20: INTEGER IPIV( * ), IWORK( * ) 21: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ), 22: $ C( * ) 23: * .. 24: * 25: * Purpose 26: * ======= 27: * 28: * DLA_GERCOND estimates the Skeel condition number of op(A) * op2(C) 29: * where op2 is determined by CMODE as follows 30: * CMODE = 1 op2(C) = C 31: * CMODE = 0 op2(C) = I 32: * CMODE = -1 op2(C) = inv(C) 33: * The Skeel condition number cond(A) = norminf( |inv(A)||A| ) 34: * is computed by computing scaling factors R such that 35: * diag(R)*A*op2(C) is row equilibrated and computing the standard 36: * infinity-norm condition number. 37: * 38: * Arguments 39: * ========== 40: * 41: * TRANS (input) CHARACTER*1 42: * Specifies the form of the system of equations: 43: * = 'N': A * X = B (No transpose) 44: * = 'T': A**T * X = B (Transpose) 45: * = 'C': A**H * X = B (Conjugate Transpose = Transpose) 46: * 47: * N (input) INTEGER 48: * The number of linear equations, i.e., the order of the 49: * matrix A. N >= 0. 50: * 51: * A (input) DOUBLE PRECISION array, dimension (LDA,N) 52: * On entry, the N-by-N matrix A. 53: * 54: * LDA (input) INTEGER 55: * The leading dimension of the array A. LDA >= max(1,N). 56: * 57: * AF (input) DOUBLE PRECISION array, dimension (LDAF,N) 58: * The factors L and U from the factorization 59: * A = P*L*U as computed by DGETRF. 60: * 61: * LDAF (input) INTEGER 62: * The leading dimension of the array AF. LDAF >= max(1,N). 63: * 64: * IPIV (input) INTEGER array, dimension (N) 65: * The pivot indices from the factorization A = P*L*U 66: * as computed by DGETRF; row i of the matrix was interchanged 67: * with row IPIV(i). 68: * 69: * CMODE (input) INTEGER 70: * Determines op2(C) in the formula op(A) * op2(C) as follows: 71: * CMODE = 1 op2(C) = C 72: * CMODE = 0 op2(C) = I 73: * CMODE = -1 op2(C) = inv(C) 74: * 75: * C (input) DOUBLE PRECISION array, dimension (N) 76: * The vector C in the formula op(A) * op2(C). 77: * 78: * INFO (output) INTEGER 79: * = 0: Successful exit. 80: * i > 0: The ith argument is invalid. 81: * 82: * WORK (input) DOUBLE PRECISION array, dimension (3*N). 83: * Workspace. 84: * 85: * IWORK (input) INTEGER array, dimension (N). 86: * Workspace. 87: * 88: * ===================================================================== 89: * 90: * .. Local Scalars .. 91: LOGICAL NOTRANS 92: INTEGER KASE, I, J 93: DOUBLE PRECISION AINVNM, TMP 94: * .. 95: * .. Local Arrays .. 96: INTEGER ISAVE( 3 ) 97: * .. 98: * .. External Functions .. 99: LOGICAL LSAME 100: EXTERNAL LSAME 101: * .. 102: * .. External Subroutines .. 103: EXTERNAL DLACN2, DGETRS, XERBLA 104: * .. 105: * .. Intrinsic Functions .. 106: INTRINSIC ABS, MAX 107: * .. 108: * .. Executable Statements .. 109: * 110: DLA_GERCOND = 0.0D+0 111: * 112: INFO = 0 113: NOTRANS = LSAME( TRANS, 'N' ) 114: IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T') 115: $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN 116: INFO = -1 117: ELSE IF( N.LT.0 ) THEN 118: INFO = -2 119: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 120: INFO = -4 121: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 122: INFO = -6 123: END IF 124: IF( INFO.NE.0 ) THEN 125: CALL XERBLA( 'DLA_GERCOND', -INFO ) 126: RETURN 127: END IF 128: IF( N.EQ.0 ) THEN 129: DLA_GERCOND = 1.0D+0 130: RETURN 131: END IF 132: * 133: * Compute the equilibration matrix R such that 134: * inv(R)*A*C has unit 1-norm. 135: * 136: IF (NOTRANS) THEN 137: DO I = 1, N 138: TMP = 0.0D+0 139: IF ( CMODE .EQ. 1 ) THEN 140: DO J = 1, N 141: TMP = TMP + ABS( A( I, J ) * C( J ) ) 142: END DO 143: ELSE IF ( CMODE .EQ. 0 ) THEN 144: DO J = 1, N 145: TMP = TMP + ABS( A( I, J ) ) 146: END DO 147: ELSE 148: DO J = 1, N 149: TMP = TMP + ABS( A( I, J ) / C( J ) ) 150: END DO 151: END IF 152: WORK( 2*N+I ) = TMP 153: END DO 154: ELSE 155: DO I = 1, N 156: TMP = 0.0D+0 157: IF ( CMODE .EQ. 1 ) THEN 158: DO J = 1, N 159: TMP = TMP + ABS( A( J, I ) * C( J ) ) 160: END DO 161: ELSE IF ( CMODE .EQ. 0 ) THEN 162: DO J = 1, N 163: TMP = TMP + ABS( A( J, I ) ) 164: END DO 165: ELSE 166: DO J = 1, N 167: TMP = TMP + ABS( A( J, I ) / C( J ) ) 168: END DO 169: END IF 170: WORK( 2*N+I ) = TMP 171: END DO 172: END IF 173: * 174: * Estimate the norm of inv(op(A)). 175: * 176: AINVNM = 0.0D+0 177: 178: KASE = 0 179: 10 CONTINUE 180: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 181: IF( KASE.NE.0 ) THEN 182: IF( KASE.EQ.2 ) THEN 183: * 184: * Multiply by R. 185: * 186: DO I = 1, N 187: WORK(I) = WORK(I) * WORK(2*N+I) 188: END DO 189: 190: IF (NOTRANS) THEN 191: CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, 192: $ WORK, N, INFO ) 193: ELSE 194: CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV, 195: $ WORK, N, INFO ) 196: END IF 197: * 198: * Multiply by inv(C). 199: * 200: IF ( CMODE .EQ. 1 ) THEN 201: DO I = 1, N 202: WORK( I ) = WORK( I ) / C( I ) 203: END DO 204: ELSE IF ( CMODE .EQ. -1 ) THEN 205: DO I = 1, N 206: WORK( I ) = WORK( I ) * C( I ) 207: END DO 208: END IF 209: ELSE 210: * 211: * Multiply by inv(C'). 212: * 213: IF ( CMODE .EQ. 1 ) THEN 214: DO I = 1, N 215: WORK( I ) = WORK( I ) / C( I ) 216: END DO 217: ELSE IF ( CMODE .EQ. -1 ) THEN 218: DO I = 1, N 219: WORK( I ) = WORK( I ) * C( I ) 220: END DO 221: END IF 222: 223: IF (NOTRANS) THEN 224: CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV, 225: $ WORK, N, INFO ) 226: ELSE 227: CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, 228: $ WORK, N, INFO ) 229: END IF 230: * 231: * Multiply by R. 232: * 233: DO I = 1, N 234: WORK( I ) = WORK( I ) * WORK( 2*N+I ) 235: END DO 236: END IF 237: GO TO 10 238: END IF 239: * 240: * Compute the estimate of the reciprocal condition number. 241: * 242: IF( AINVNM .NE. 0.0D+0 ) 243: $ DLA_GERCOND = ( 1.0D+0 / AINVNM ) 244: * 245: RETURN 246: * 247: END