![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: DOUBLE PRECISION FUNCTION ZLA_HERCOND_C( UPLO, 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 Arguments .. 16: CHARACTER UPLO 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_HERCOND_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: * UPLO (input) CHARACTER*1 36: * = 'U': Upper triangle of A is stored; 37: * = 'L': Lower triangle of A is stored. 38: * 39: * N (input) INTEGER 40: * The number of linear equations, i.e., the order of the 41: * matrix A. N >= 0. 42: * 43: * A (input) COMPLEX*16 array, dimension (LDA,N) 44: * On entry, the N-by-N matrix A 45: * 46: * LDA (input) INTEGER 47: * The leading dimension of the array A. LDA >= max(1,N). 48: * 49: * AF (input) COMPLEX*16 array, dimension (LDAF,N) 50: * The block diagonal matrix D and the multipliers used to 51: * obtain the factor U or L as computed by ZHETRF. 52: * 53: * LDAF (input) INTEGER 54: * The leading dimension of the array AF. LDAF >= max(1,N). 55: * 56: * IPIV (input) INTEGER array, dimension (N) 57: * Details of the interchanges and the block structure of D 58: * as determined by CHETRF. 59: * 60: * C (input) DOUBLE PRECISION array, dimension (N) 61: * The vector C in the formula op(A) * inv(diag(C)). 62: * 63: * CAPPLY (input) LOGICAL 64: * If .TRUE. then access the vector C in the formula above. 65: * 66: * INFO (output) INTEGER 67: * = 0: Successful exit. 68: * i > 0: The ith argument is invalid. 69: * 70: * WORK (input) COMPLEX*16 array, dimension (2*N). 71: * Workspace. 72: * 73: * RWORK (input) DOUBLE PRECISION array, dimension (N). 74: * Workspace. 75: * 76: * ===================================================================== 77: * 78: * .. Local Scalars .. 79: INTEGER KASE, I, J 80: DOUBLE PRECISION AINVNM, ANORM, TMP 81: LOGICAL UP 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, ZHETRS, XERBLA 93: * .. 94: * .. Intrinsic Functions .. 95: INTRINSIC ABS, MAX 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_HERCOND_C = 0.0D+0 106: * 107: INFO = 0 108: IF( N.LT.0 ) THEN 109: INFO = -2 110: END IF 111: IF( INFO.NE.0 ) THEN 112: CALL XERBLA( 'ZLA_HERCOND_C', -INFO ) 113: RETURN 114: END IF 115: UP = .FALSE. 116: IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE. 117: * 118: * Compute norm of op(A)*op2(C). 119: * 120: ANORM = 0.0D+0 121: IF ( UP ) THEN 122: DO I = 1, N 123: TMP = 0.0D+0 124: IF ( CAPPLY ) THEN 125: DO J = 1, I 126: TMP = TMP + CABS1( A( J, I ) ) / C( J ) 127: END DO 128: DO J = I+1, N 129: TMP = TMP + CABS1( A( I, J ) ) / C( J ) 130: END DO 131: ELSE 132: DO J = 1, I 133: TMP = TMP + CABS1( A( J, I ) ) 134: END DO 135: DO J = I+1, N 136: TMP = TMP + CABS1( A( I, J ) ) 137: END DO 138: END IF 139: RWORK( I ) = TMP 140: ANORM = MAX( ANORM, TMP ) 141: END DO 142: ELSE 143: DO I = 1, N 144: TMP = 0.0D+0 145: IF ( CAPPLY ) THEN 146: DO J = 1, I 147: TMP = TMP + CABS1( A( I, J ) ) / C( J ) 148: END DO 149: DO J = I+1, N 150: TMP = TMP + CABS1( A( J, I ) ) / C( J ) 151: END DO 152: ELSE 153: DO J = 1, I 154: TMP = TMP + CABS1( A( I, J ) ) 155: END DO 156: DO J = I+1, N 157: TMP = TMP + CABS1( A( J, I ) ) 158: END DO 159: END IF 160: RWORK( I ) = TMP 161: ANORM = MAX( ANORM, TMP ) 162: END DO 163: END IF 164: * 165: * Quick return if possible. 166: * 167: IF( N.EQ.0 ) THEN 168: ZLA_HERCOND_C = 1.0D+0 169: RETURN 170: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN 171: RETURN 172: END IF 173: * 174: * Estimate the norm of inv(op(A)). 175: * 176: AINVNM = 0.0D+0 177: * 178: KASE = 0 179: 10 CONTINUE 180: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 181: IF( KASE.NE.0 ) THEN 182: IF( KASE.EQ.2 ) THEN 183: * 184: * Multiply by R. 185: * 186: DO I = 1, N 187: WORK( I ) = WORK( I ) * RWORK( I ) 188: END DO 189: * 190: IF ( UP ) THEN 191: CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV, 192: $ WORK, N, INFO ) 193: ELSE 194: CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV, 195: $ WORK, N, INFO ) 196: ENDIF 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: ELSE 206: * 207: * Multiply by inv(C'). 208: * 209: IF ( CAPPLY ) THEN 210: DO I = 1, N 211: WORK( I ) = WORK( I ) * C( I ) 212: END DO 213: END IF 214: * 215: IF ( UP ) THEN 216: CALL ZHETRS( 'U', N, 1, AF, LDAF, IPIV, 217: $ WORK, N, INFO ) 218: ELSE 219: CALL ZHETRS( 'L', N, 1, AF, LDAF, IPIV, 220: $ WORK, N, INFO ) 221: END IF 222: * 223: * Multiply by R. 224: * 225: DO I = 1, N 226: WORK( I ) = WORK( I ) * RWORK( I ) 227: END DO 228: END IF 229: GO TO 10 230: END IF 231: * 232: * Compute the estimate of the reciprocal condition number. 233: * 234: IF( AINVNM .NE. 0.0D+0 ) 235: $ ZLA_HERCOND_C = 1.0D+0 / AINVNM 236: * 237: RETURN 238: * 239: END