![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: DOUBLE PRECISION FUNCTION DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB, 2: $ AFB, LDAFB, IPIV, CMODE, C, 3: $ INFO, WORK, IWORK ) 4: * 5: * -- LAPACK routine (version 3.2.2) -- 6: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- 7: * -- Jason Riedy of Univ. of California Berkeley. -- 8: * -- June 2010 -- 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, LDAB, LDAFB, INFO, KL, KU, CMODE 18: * .. 19: * .. Array Arguments .. 20: INTEGER IWORK( * ), IPIV( * ) 21: DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ), 22: $ C( * ) 23: * .. 24: * 25: * Purpose 26: * ======= 27: * 28: * DLA_GBRCOND 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: * KL (input) INTEGER 52: * The number of subdiagonals within the band of A. KL >= 0. 53: * 54: * KU (input) INTEGER 55: * The number of superdiagonals within the band of A. KU >= 0. 56: * 57: * AB (input) DOUBLE PRECISION array, dimension (LDAB,N) 58: * On entry, the matrix A in band storage, in rows 1 to KL+KU+1. 59: * The j-th column of A is stored in the j-th column of the 60: * array AB as follows: 61: * AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) 62: * 63: * LDAB (input) INTEGER 64: * The leading dimension of the array AB. LDAB >= KL+KU+1. 65: * 66: * AFB (input) DOUBLE PRECISION array, dimension (LDAFB,N) 67: * Details of the LU factorization of the band matrix A, as 68: * computed by DGBTRF. U is stored as an upper triangular 69: * band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, 70: * and the multipliers used during the factorization are stored 71: * in rows KL+KU+2 to 2*KL+KU+1. 72: * 73: * LDAFB (input) INTEGER 74: * The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1. 75: * 76: * IPIV (input) INTEGER array, dimension (N) 77: * The pivot indices from the factorization A = P*L*U 78: * as computed by DGBTRF; row i of the matrix was interchanged 79: * with row IPIV(i). 80: * 81: * CMODE (input) INTEGER 82: * Determines op2(C) in the formula op(A) * op2(C) as follows: 83: * CMODE = 1 op2(C) = C 84: * CMODE = 0 op2(C) = I 85: * CMODE = -1 op2(C) = inv(C) 86: * 87: * C (input) DOUBLE PRECISION array, dimension (N) 88: * The vector C in the formula op(A) * op2(C). 89: * 90: * INFO (output) INTEGER 91: * = 0: Successful exit. 92: * i > 0: The ith argument is invalid. 93: * 94: * WORK (input) DOUBLE PRECISION array, dimension (5*N). 95: * Workspace. 96: * 97: * IWORK (input) INTEGER array, dimension (N). 98: * Workspace. 99: * 100: * ===================================================================== 101: * 102: * .. Local Scalars .. 103: LOGICAL NOTRANS 104: INTEGER KASE, I, J, KD, KE 105: DOUBLE PRECISION AINVNM, TMP 106: * .. 107: * .. Local Arrays .. 108: INTEGER ISAVE( 3 ) 109: * .. 110: * .. External Functions .. 111: LOGICAL LSAME 112: EXTERNAL LSAME 113: * .. 114: * .. External Subroutines .. 115: EXTERNAL DLACN2, DGBTRS, XERBLA 116: * .. 117: * .. Intrinsic Functions .. 118: INTRINSIC ABS, MAX 119: * .. 120: * .. Executable Statements .. 121: * 122: DLA_GBRCOND = 0.0D+0 123: * 124: INFO = 0 125: NOTRANS = LSAME( TRANS, 'N' ) 126: IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T') 127: $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN 128: INFO = -1 129: ELSE IF( N.LT.0 ) THEN 130: INFO = -2 131: ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN 132: INFO = -3 133: ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 134: INFO = -4 135: ELSE IF( LDAB.LT.KL+KU+1 ) THEN 136: INFO = -6 137: ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 138: INFO = -8 139: END IF 140: IF( INFO.NE.0 ) THEN 141: CALL XERBLA( 'DLA_GBRCOND', -INFO ) 142: RETURN 143: END IF 144: IF( N.EQ.0 ) THEN 145: DLA_GBRCOND = 1.0D+0 146: RETURN 147: END IF 148: * 149: * Compute the equilibration matrix R such that 150: * inv(R)*A*C has unit 1-norm. 151: * 152: KD = KU + 1 153: KE = KL + 1 154: IF ( NOTRANS ) THEN 155: DO I = 1, N 156: TMP = 0.0D+0 157: IF ( CMODE .EQ. 1 ) THEN 158: DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 159: TMP = TMP + ABS( AB( KD+I-J, J ) * C( J ) ) 160: END DO 161: ELSE IF ( CMODE .EQ. 0 ) THEN 162: DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 163: TMP = TMP + ABS( AB( KD+I-J, J ) ) 164: END DO 165: ELSE 166: DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 167: TMP = TMP + ABS( AB( KD+I-J, J ) / C( J ) ) 168: END DO 169: END IF 170: WORK( 2*N+I ) = TMP 171: END DO 172: ELSE 173: DO I = 1, N 174: TMP = 0.0D+0 175: IF ( CMODE .EQ. 1 ) THEN 176: DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 177: TMP = TMP + ABS( AB( KE-I+J, I ) * C( J ) ) 178: END DO 179: ELSE IF ( CMODE .EQ. 0 ) THEN 180: DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 181: TMP = TMP + ABS( AB( KE-I+J, I ) ) 182: END DO 183: ELSE 184: DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 185: TMP = TMP + ABS( AB( KE-I+J, I ) / C( J ) ) 186: END DO 187: END IF 188: WORK( 2*N+I ) = TMP 189: END DO 190: END IF 191: * 192: * Estimate the norm of inv(op(A)). 193: * 194: AINVNM = 0.0D+0 195: 196: KASE = 0 197: 10 CONTINUE 198: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 199: IF( KASE.NE.0 ) THEN 200: IF( KASE.EQ.2 ) THEN 201: * 202: * Multiply by R. 203: * 204: DO I = 1, N 205: WORK( I ) = WORK( I ) * WORK( 2*N+I ) 206: END DO 207: 208: IF ( NOTRANS ) THEN 209: CALL DGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 210: $ IPIV, WORK, N, INFO ) 211: ELSE 212: CALL DGBTRS( 'Transpose', N, KL, KU, 1, AFB, LDAFB, IPIV, 213: $ WORK, N, INFO ) 214: END IF 215: * 216: * Multiply by inv(C). 217: * 218: IF ( CMODE .EQ. 1 ) THEN 219: DO I = 1, N 220: WORK( I ) = WORK( I ) / C( I ) 221: END DO 222: ELSE IF ( CMODE .EQ. -1 ) THEN 223: DO I = 1, N 224: WORK( I ) = WORK( I ) * C( I ) 225: END DO 226: END IF 227: ELSE 228: * 229: * Multiply by inv(C'). 230: * 231: IF ( CMODE .EQ. 1 ) THEN 232: DO I = 1, N 233: WORK( I ) = WORK( I ) / C( I ) 234: END DO 235: ELSE IF ( CMODE .EQ. -1 ) THEN 236: DO I = 1, N 237: WORK( I ) = WORK( I ) * C( I ) 238: END DO 239: END IF 240: 241: IF ( NOTRANS ) THEN 242: CALL DGBTRS( 'Transpose', N, KL, KU, 1, AFB, LDAFB, IPIV, 243: $ WORK, N, INFO ) 244: ELSE 245: CALL DGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 246: $ IPIV, WORK, N, INFO ) 247: END IF 248: * 249: * Multiply by R. 250: * 251: DO I = 1, N 252: WORK( I ) = WORK( I ) * WORK( 2*N+I ) 253: END DO 254: END IF 255: GO TO 10 256: END IF 257: * 258: * Compute the estimate of the reciprocal condition number. 259: * 260: IF( AINVNM .NE. 0.0D+0 ) 261: $ DLA_GBRCOND = ( 1.0D+0 / AINVNM ) 262: * 263: RETURN 264: * 265: END