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