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