![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: DOUBLE PRECISION FUNCTION ZLA_PORCOND_C( UPLO, N, A, LDA, AF, 2: $ LDAF, C, CAPPLY, 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: LOGICAL CAPPLY 18: INTEGER N, LDA, LDAF, INFO 19: * .. 20: * .. Array Arguments .. 21: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), WORK( * ) 22: DOUBLE PRECISION C( * ), RWORK( * ) 23: * .. 24: * 25: * Purpose 26: * ======= 27: * 28: * ZLA_PORCOND_C Computes the infinity norm condition number of 29: * op(A) * inv(diag(C)) where C is a DOUBLE PRECISION 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 triangular factor U or L from the Cholesky factorization 50: * A = U**T*U or A = L*L**T, as computed by ZPOTRF. 51: * 52: * LDAF (input) INTEGER 53: * The leading dimension of the array AF. LDAF >= max(1,N). 54: * 55: * C (input) DOUBLE PRECISION array, dimension (N) 56: * The vector C in the formula op(A) * inv(diag(C)). 57: * 58: * CAPPLY (input) LOGICAL 59: * If .TRUE. then access the vector C in the formula above. 60: * 61: * INFO (output) INTEGER 62: * = 0: Successful exit. 63: * i > 0: The ith argument is invalid. 64: * 65: * WORK (input) COMPLEX*16 array, dimension (2*N). 66: * Workspace. 67: * 68: * RWORK (input) DOUBLE PRECISION array, dimension (N). 69: * Workspace. 70: * 71: * ===================================================================== 72: * 73: * .. Local Scalars .. 74: INTEGER KASE 75: DOUBLE PRECISION AINVNM, ANORM, TMP 76: INTEGER I, J 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, ZPOTRS, XERBLA 89: * .. 90: * .. Intrinsic Functions .. 91: INTRINSIC ABS, MAX, REAL, DIMAG 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_PORCOND_C = 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_PORCOND_C', -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: IF ( CAPPLY ) THEN 121: DO J = 1, I 122: TMP = TMP + CABS1( A( J, I ) ) / C( J ) 123: END DO 124: DO J = I+1, N 125: TMP = TMP + CABS1( A( I, J ) ) / C( J ) 126: END DO 127: ELSE 128: DO J = 1, I 129: TMP = TMP + CABS1( A( J, I ) ) 130: END DO 131: DO J = I+1, N 132: TMP = TMP + CABS1( A( I, J ) ) 133: END DO 134: END IF 135: RWORK( I ) = TMP 136: ANORM = MAX( ANORM, TMP ) 137: END DO 138: ELSE 139: DO I = 1, N 140: TMP = 0.0D+0 141: IF ( CAPPLY ) THEN 142: DO J = 1, I 143: TMP = TMP + CABS1( A( I, J ) ) / C( J ) 144: END DO 145: DO J = I+1, N 146: TMP = TMP + CABS1( A( J, I ) ) / C( J ) 147: END DO 148: ELSE 149: DO J = 1, I 150: TMP = TMP + CABS1( A( I, J ) ) 151: END DO 152: DO J = I+1, N 153: TMP = TMP + CABS1( A( J, I ) ) 154: END DO 155: END IF 156: RWORK( I ) = TMP 157: ANORM = MAX( ANORM, TMP ) 158: END DO 159: END IF 160: * 161: * Quick return if possible. 162: * 163: IF( N.EQ.0 ) THEN 164: ZLA_PORCOND_C = 1.0D+0 165: RETURN 166: ELSE IF( ANORM .EQ. 0.0D+0 ) THEN 167: RETURN 168: END IF 169: * 170: * Estimate the norm of inv(op(A)). 171: * 172: AINVNM = 0.0D+0 173: * 174: KASE = 0 175: 10 CONTINUE 176: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 177: IF( KASE.NE.0 ) THEN 178: IF( KASE.EQ.2 ) THEN 179: * 180: * Multiply by R. 181: * 182: DO I = 1, N 183: WORK( I ) = WORK( I ) * RWORK( I ) 184: END DO 185: * 186: IF ( UP ) THEN 187: CALL ZPOTRS( 'U', N, 1, AF, LDAF, 188: $ WORK, N, INFO ) 189: ELSE 190: CALL ZPOTRS( 'L', N, 1, AF, LDAF, 191: $ WORK, N, INFO ) 192: ENDIF 193: * 194: * Multiply by inv(C). 195: * 196: IF ( CAPPLY ) THEN 197: DO I = 1, N 198: WORK( I ) = WORK( I ) * C( I ) 199: END DO 200: END IF 201: ELSE 202: * 203: * Multiply by inv(C'). 204: * 205: IF ( CAPPLY ) THEN 206: DO I = 1, N 207: WORK( I ) = WORK( I ) * C( I ) 208: END DO 209: END IF 210: * 211: IF ( UP ) THEN 212: CALL ZPOTRS( 'U', N, 1, AF, LDAF, 213: $ WORK, N, INFO ) 214: ELSE 215: CALL ZPOTRS( 'L', N, 1, AF, LDAF, 216: $ WORK, N, INFO ) 217: END IF 218: * 219: * Multiply by R. 220: * 221: DO I = 1, N 222: WORK( I ) = WORK( I ) * RWORK( I ) 223: END DO 224: END IF 225: GO TO 10 226: END IF 227: * 228: * Compute the estimate of the reciprocal condition number. 229: * 230: IF( AINVNM .NE. 0.0D+0 ) 231: $ ZLA_PORCOND_C = 1.0D+0 / AINVNM 232: * 233: RETURN 234: * 235: END