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