![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: DOUBLE PRECISION FUNCTION ZLA_HERCOND_X( UPLO, 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 UPLO 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_HERCOND_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: * UPLO (input) CHARACTER*1 35: * = 'U': Upper triangle of A is stored; 36: * = 'L': Lower triangle of A is stored. 37: * 38: * N (input) INTEGER 39: * The number of linear equations, i.e., the order of the 40: * matrix A. N >= 0. 41: * 42: * A (input) COMPLEX*16 array, dimension (LDA,N) 43: * On entry, the N-by-N matrix A. 44: * 45: * LDA (input) INTEGER 46: * The leading dimension of the array A. LDA >= max(1,N). 47: * 48: * AF (input) COMPLEX*16 array, dimension (LDAF,N) 49: * The block diagonal matrix D and the multipliers used to 50: * obtain the factor U or L as computed by ZHETRF. 51: * 52: * LDAF (input) INTEGER 53: * The leading dimension of the array AF. LDAF >= max(1,N). 54: * 55: * IPIV (input) INTEGER array, dimension (N) 56: * Details of the interchanges and the block structure of D 57: * as determined by CHETRF. 58: * 59: * X (input) COMPLEX*16 array, dimension (N) 60: * The vector X in the formula op(A) * diag(X). 61: * 62: * INFO (output) INTEGER 63: * = 0: Successful exit. 64: * i > 0: The ith argument is invalid. 65: * 66: * WORK (input) COMPLEX*16 array, dimension (2*N). 67: * Workspace. 68: * 69: * RWORK (input) DOUBLE PRECISION array, dimension (N). 70: * Workspace. 71: * 72: * ===================================================================== 73: * 74: * .. Local Scalars .. 75: INTEGER KASE, I, J 76: DOUBLE PRECISION AINVNM, ANORM, TMP 77: LOGICAL UP 78: COMPLEX*16 ZDUM 79: * .. 80: * .. Local Arrays .. 81: INTEGER ISAVE( 3 ) 82: * .. 83: * .. External Functions .. 84: LOGICAL LSAME 85: EXTERNAL LSAME 86: * .. 87: * .. External Subroutines .. 88: EXTERNAL ZLACN2, ZHETRS, XERBLA 89: * .. 90: * .. Intrinsic Functions .. 91: INTRINSIC ABS, MAX 92: * .. 93: * .. Statement Functions .. 94: DOUBLE PRECISION CABS1 95: * .. 96: * .. Statement Function Definitions .. 97: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 98: * .. 99: * .. Executable Statements .. 100: * 101: ZLA_HERCOND_X = 0.0D+0 102: * 103: INFO = 0 104: IF( N.LT.0 ) THEN 105: INFO = -2 106: END IF 107: IF( INFO.NE.0 ) THEN 108: CALL XERBLA( 'ZLA_HERCOND_X', -INFO ) 109: RETURN 110: END IF 111: UP = .FALSE. 112: IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE. 113: * 114: * Compute norm of op(A)*op2(C). 115: * 116: ANORM = 0.0D+0 117: IF ( UP ) THEN 118: DO I = 1, N 119: TMP = 0.0D+0 120: DO J = 1, I 121: TMP = TMP + CABS1( A( J, I ) * X( J ) ) 122: END DO 123: DO J = I+1, N 124: TMP = TMP + CABS1( A( I, J ) * X( J ) ) 125: END DO 126: RWORK( I ) = TMP 127: ANORM = MAX( ANORM, TMP ) 128: END DO 129: ELSE 130: DO I = 1, N 131: TMP = 0.0D+0 132: DO J = 1, I 133: TMP = TMP + CABS1( A( I, J ) * X( J ) ) 134: END DO 135: DO J = I+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_HERCOND_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: * 162: * Multiply by R. 163: * 164: DO I = 1, N 165: WORK( I ) = WORK( I ) * RWORK( I ) 166: END DO 167: * 168: IF ( UP ) THEN 169: CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV, 170: $ WORK, N, INFO ) 171: ELSE 172: CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV, 173: $ WORK, N, INFO ) 174: ENDIF 175: * 176: * Multiply by inv(X). 177: * 178: DO I = 1, N 179: WORK( I ) = WORK( I ) / X( I ) 180: END DO 181: ELSE 182: * 183: * Multiply by inv(X'). 184: * 185: DO I = 1, N 186: WORK( I ) = WORK( I ) / X( I ) 187: END DO 188: * 189: IF ( UP ) THEN 190: CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV, 191: $ WORK, N, INFO ) 192: ELSE 193: CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV, 194: $ WORK, N, INFO ) 195: END IF 196: * 197: * Multiply by R. 198: * 199: DO I = 1, N 200: WORK( I ) = WORK( I ) * RWORK( I ) 201: END DO 202: END IF 203: GO TO 10 204: END IF 205: * 206: * Compute the estimate of the reciprocal condition number. 207: * 208: IF( AINVNM .NE. 0.0D+0 ) 209: $ ZLA_HERCOND_X = 1.0D+0 / AINVNM 210: * 211: RETURN 212: * 213: END