![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: DOUBLE PRECISION FUNCTION ZLA_GBRCOND_X( TRANS, N, KL, KU, AB, 2: $ LDAB, AFB, LDAFB, IPIV, 3: $ X, INFO, WORK, RWORK ) 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, KL, KU, KD, KE, LDAB, LDAFB, INFO 18: * .. 19: * .. Array Arguments .. 20: INTEGER IPIV( * ) 21: COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ), 22: $ X( * ) 23: DOUBLE PRECISION RWORK( * ) 24: * 25: * 26: * Purpose 27: * ======= 28: * 29: * ZLA_GBRCOND_X Computes the infinity norm condition number of 30: * op(A) * diag(X) where X is a COMPLEX*16 vector. 31: * 32: * Arguments 33: * ========= 34: * 35: * TRANS (input) CHARACTER*1 36: * Specifies the form of the system of equations: 37: * = 'N': A * X = B (No transpose) 38: * = 'T': A**T * X = B (Transpose) 39: * = 'C': A**H * X = B (Conjugate Transpose = Transpose) 40: * 41: * N (input) INTEGER 42: * The number of linear equations, i.e., the order of the 43: * matrix A. N >= 0. 44: * 45: * KL (input) INTEGER 46: * The number of subdiagonals within the band of A. KL >= 0. 47: * 48: * KU (input) INTEGER 49: * The number of superdiagonals within the band of A. KU >= 0. 50: * 51: * AB (input) COMPLEX*16 array, dimension (LDAB,N) 52: * On entry, the matrix A in band storage, in rows 1 to KL+KU+1. 53: * The j-th column of A is stored in the j-th column of the 54: * array AB as follows: 55: * AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) 56: * 57: * LDAB (input) INTEGER 58: * The leading dimension of the array AB. LDAB >= KL+KU+1. 59: * 60: * AFB (input) COMPLEX*16 array, dimension (LDAFB,N) 61: * Details of the LU factorization of the band matrix A, as 62: * computed by ZGBTRF. U is stored as an upper triangular 63: * band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, 64: * and the multipliers used during the factorization are stored 65: * in rows KL+KU+2 to 2*KL+KU+1. 66: * 67: * LDAFB (input) INTEGER 68: * The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1. 69: * 70: * IPIV (input) INTEGER array, dimension (N) 71: * The pivot indices from the factorization A = P*L*U 72: * as computed by ZGBTRF; row i of the matrix was interchanged 73: * with row IPIV(i). 74: * 75: * X (input) COMPLEX*16 array, dimension (N) 76: * The vector X in the formula op(A) * diag(X). 77: * 78: * INFO (output) INTEGER 79: * = 0: Successful exit. 80: * i > 0: The ith argument is invalid. 81: * 82: * WORK (input) COMPLEX*16 array, dimension (2*N). 83: * Workspace. 84: * 85: * RWORK (input) DOUBLE PRECISION array, dimension (N). 86: * Workspace. 87: * 88: * ===================================================================== 89: * 90: * .. Local Scalars .. 91: LOGICAL NOTRANS 92: INTEGER KASE, I, J 93: DOUBLE PRECISION AINVNM, ANORM, TMP 94: COMPLEX*16 ZDUM 95: * .. 96: * .. Local Arrays .. 97: INTEGER ISAVE( 3 ) 98: * .. 99: * .. External Functions .. 100: LOGICAL LSAME 101: EXTERNAL LSAME 102: * .. 103: * .. External Subroutines .. 104: EXTERNAL ZLACN2, ZGBTRS, XERBLA 105: * .. 106: * .. Intrinsic Functions .. 107: INTRINSIC ABS, MAX 108: * .. 109: * .. Statement Functions .. 110: DOUBLE PRECISION CABS1 111: * .. 112: * .. Statement Function Definitions .. 113: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 114: * .. 115: * .. Executable Statements .. 116: * 117: ZLA_GBRCOND_X = 0.0D+0 118: * 119: INFO = 0 120: NOTRANS = LSAME( TRANS, 'N' ) 121: IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T') .AND. .NOT. 122: $ LSAME( TRANS, 'C' ) ) THEN 123: INFO = -1 124: ELSE IF( N.LT.0 ) THEN 125: INFO = -2 126: ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN 127: INFO = -3 128: ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 129: INFO = -4 130: ELSE IF( LDAB.LT.KL+KU+1 ) THEN 131: INFO = -6 132: ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 133: INFO = -8 134: END IF 135: IF( INFO.NE.0 ) THEN 136: CALL XERBLA( 'ZLA_GBRCOND_X', -INFO ) 137: RETURN 138: END IF 139: * 140: * Compute norm of op(A)*op2(C). 141: * 142: KD = KU + 1 143: KE = KL + 1 144: ANORM = 0.0D+0 145: IF ( NOTRANS ) THEN 146: DO I = 1, N 147: TMP = 0.0D+0 148: DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 149: TMP = TMP + CABS1( AB( KD+I-J, J) * X( J ) ) 150: END DO 151: RWORK( I ) = TMP 152: ANORM = MAX( ANORM, TMP ) 153: END DO 154: ELSE 155: DO I = 1, N 156: TMP = 0.0D+0 157: DO J = MAX( I-KL, 1 ), MIN( I+KU, N ) 158: TMP = TMP + CABS1( AB( KE-I+J, I ) * X( J ) ) 159: END DO 160: RWORK( I ) = TMP 161: ANORM = MAX( ANORM, TMP ) 162: END DO 163: END IF 164: * 165: * Quick return if possible. 166: * 167: IF( N.EQ.0 ) THEN 168: ZLA_GBRCOND_X = 1.0D+0 169: RETURN 170: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN 171: RETURN 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 ZLACN2( N, WORK( N+1 ), WORK, 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 ) * RWORK( I ) 188: END DO 189: * 190: IF ( NOTRANS ) THEN 191: CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 192: $ IPIV, WORK, N, INFO ) 193: ELSE 194: CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB, 195: $ LDAFB, IPIV, WORK, N, INFO ) 196: ENDIF 197: * 198: * Multiply by inv(X). 199: * 200: DO I = 1, N 201: WORK( I ) = WORK( I ) / X( I ) 202: END DO 203: ELSE 204: * 205: * Multiply by inv(X'). 206: * 207: DO I = 1, N 208: WORK( I ) = WORK( I ) / X( I ) 209: END DO 210: * 211: IF ( NOTRANS ) THEN 212: CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB, 213: $ LDAFB, IPIV, WORK, N, INFO ) 214: ELSE 215: CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB, 216: $ IPIV, WORK, N, INFO ) 217: END IF 218: * 219: * Multiply by R. 220: * 221: DO I = 1, N 222: WORK( I ) = WORK( I ) * RWORK( I ) 223: END DO 224: END IF 225: GO TO 10 226: END IF 227: * 228: * Compute the estimate of the reciprocal condition number. 229: * 230: IF( AINVNM .NE. 0.0D+0 ) 231: $ ZLA_GBRCOND_X = 1.0D+0 / AINVNM 232: * 233: RETURN 234: * 235: END