![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: DOUBLE PRECISION FUNCTION ZLA_GERCOND_X( TRANS, N, A, LDA, AF, 2: $ LDAF, IPIV, X, INFO, 3: $ 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, LDA, LDAF, INFO 18: * .. 19: * .. Array Arguments .. 20: INTEGER IPIV( * ) 21: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * ) 22: DOUBLE PRECISION RWORK( * ) 23: * .. 24: * 25: * Purpose 26: * ======= 27: * 28: * ZLA_GERCOND_X computes the infinity norm condition number of 29: * op(A) * diag(X) where X is a COMPLEX*16 vector. 30: * 31: * Arguments 32: * ========= 33: * 34: * TRANS (input) CHARACTER*1 35: * Specifies the form of the system of equations: 36: * = 'N': A * X = B (No transpose) 37: * = 'T': A**T * X = B (Transpose) 38: * = 'C': A**H * X = B (Conjugate Transpose = Transpose) 39: * 40: * N (input) INTEGER 41: * The number of linear equations, i.e., the order of the 42: * matrix A. N >= 0. 43: * 44: * A (input) COMPLEX*16 array, dimension (LDA,N) 45: * On entry, the N-by-N matrix A. 46: * 47: * LDA (input) INTEGER 48: * The leading dimension of the array A. LDA >= max(1,N). 49: * 50: * AF (input) COMPLEX*16 array, dimension (LDAF,N) 51: * The factors L and U from the factorization 52: * A = P*L*U as computed by ZGETRF. 53: * 54: * LDAF (input) INTEGER 55: * The leading dimension of the array AF. LDAF >= max(1,N). 56: * 57: * IPIV (input) INTEGER array, dimension (N) 58: * The pivot indices from the factorization A = P*L*U 59: * as computed by ZGETRF; row i of the matrix was interchanged 60: * with row IPIV(i). 61: * 62: * X (input) COMPLEX*16 array, dimension (N) 63: * The vector X in the formula op(A) * diag(X). 64: * 65: * INFO (output) INTEGER 66: * = 0: Successful exit. 67: * i > 0: The ith argument is invalid. 68: * 69: * WORK (input) COMPLEX*16 array, dimension (2*N). 70: * Workspace. 71: * 72: * RWORK (input) DOUBLE PRECISION array, dimension (N). 73: * Workspace. 74: * 75: * ===================================================================== 76: * 77: * .. Local Scalars .. 78: LOGICAL NOTRANS 79: INTEGER KASE 80: DOUBLE PRECISION AINVNM, ANORM, TMP 81: INTEGER I, J 82: COMPLEX*16 ZDUM 83: * .. 84: * .. Local Arrays .. 85: INTEGER ISAVE( 3 ) 86: * .. 87: * .. External Functions .. 88: LOGICAL LSAME 89: EXTERNAL LSAME 90: * .. 91: * .. External Subroutines .. 92: EXTERNAL ZLACN2, ZGETRS, XERBLA 93: * .. 94: * .. Intrinsic Functions .. 95: INTRINSIC ABS, MAX, REAL, DIMAG 96: * .. 97: * .. Statement Functions .. 98: DOUBLE PRECISION CABS1 99: * .. 100: * .. Statement Function Definitions .. 101: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 102: * .. 103: * .. Executable Statements .. 104: * 105: ZLA_GERCOND_X = 0.0D+0 106: * 107: INFO = 0 108: NOTRANS = LSAME( TRANS, 'N' ) 109: IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT. 110: $ LSAME( TRANS, 'C' ) ) THEN 111: INFO = -1 112: ELSE IF( N.LT.0 ) THEN 113: INFO = -2 114: END IF 115: IF( INFO.NE.0 ) THEN 116: CALL XERBLA( 'ZLA_GERCOND_X', -INFO ) 117: RETURN 118: END IF 119: * 120: * Compute norm of op(A)*op2(C). 121: * 122: ANORM = 0.0D+0 123: IF ( NOTRANS ) THEN 124: DO I = 1, N 125: TMP = 0.0D+0 126: DO J = 1, N 127: TMP = TMP + CABS1( A( I, J ) * X( J ) ) 128: END DO 129: RWORK( I ) = TMP 130: ANORM = MAX( ANORM, TMP ) 131: END DO 132: ELSE 133: DO I = 1, N 134: TMP = 0.0D+0 135: DO J = 1, N 136: TMP = TMP + CABS1( A( J, I ) * X( J ) ) 137: END DO 138: RWORK( I ) = TMP 139: ANORM = MAX( ANORM, TMP ) 140: END DO 141: END IF 142: * 143: * Quick return if possible. 144: * 145: IF( N.EQ.0 ) THEN 146: ZLA_GERCOND_X = 1.0D+0 147: RETURN 148: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN 149: RETURN 150: END IF 151: * 152: * Estimate the norm of inv(op(A)). 153: * 154: AINVNM = 0.0D+0 155: * 156: KASE = 0 157: 10 CONTINUE 158: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 159: IF( KASE.NE.0 ) THEN 160: IF( KASE.EQ.2 ) THEN 161: * Multiply by R. 162: DO I = 1, N 163: WORK( I ) = WORK( I ) * RWORK( I ) 164: END DO 165: * 166: IF ( NOTRANS ) THEN 167: CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, 168: $ WORK, N, INFO ) 169: ELSE 170: CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV, 171: $ WORK, N, INFO ) 172: ENDIF 173: * 174: * Multiply by inv(X). 175: * 176: DO I = 1, N 177: WORK( I ) = WORK( I ) / X( I ) 178: END DO 179: ELSE 180: * 181: * Multiply by inv(X'). 182: * 183: DO I = 1, N 184: WORK( I ) = WORK( I ) / X( I ) 185: END DO 186: * 187: IF ( NOTRANS ) THEN 188: CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV, 189: $ WORK, N, INFO ) 190: ELSE 191: CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV, 192: $ WORK, N, INFO ) 193: END IF 194: * 195: * Multiply by R. 196: * 197: DO I = 1, N 198: WORK( I ) = WORK( I ) * RWORK( I ) 199: END DO 200: END IF 201: GO TO 10 202: END IF 203: * 204: * Compute the estimate of the reciprocal condition number. 205: * 206: IF( AINVNM .NE. 0.0D+0 ) 207: $ ZLA_GERCOND_X = 1.0D+0 / AINVNM 208: * 209: RETURN 210: * 211: END