![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: DOUBLE PRECISION FUNCTION DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF, 2: $ CMODE, C, INFO, WORK, 3: $ IWORK ) 4: * 5: * -- LAPACK routine (version 3.2.2) -- 6: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- 7: * -- Jason Riedy of Univ. of California Berkeley. -- 8: * -- June 2010 -- 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, CMODE 18: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ), 19: $ C( * ) 20: * .. 21: * .. Array Arguments .. 22: INTEGER IWORK( * ) 23: * .. 24: * 25: * Purpose 26: * ======= 27: * 28: * DLA_PORCOND Estimates the Skeel condition number of op(A) * op2(C) 29: * where op2 is determined by CMODE as follows 30: * CMODE = 1 op2(C) = C 31: * CMODE = 0 op2(C) = I 32: * CMODE = -1 op2(C) = inv(C) 33: * The Skeel condition number cond(A) = norminf( |inv(A)||A| ) 34: * is computed by computing scaling factors R such that 35: * diag(R)*A*op2(C) is row equilibrated and computing the standard 36: * infinity-norm condition number. 37: * 38: * Arguments 39: * ========== 40: * 41: * UPLO (input) CHARACTER*1 42: * = 'U': Upper triangle of A is stored; 43: * = 'L': Lower triangle of A is stored. 44: * 45: * N (input) INTEGER 46: * The number of linear equations, i.e., the order of the 47: * matrix A. N >= 0. 48: * 49: * A (input) DOUBLE PRECISION array, dimension (LDA,N) 50: * On entry, the N-by-N matrix A. 51: * 52: * LDA (input) INTEGER 53: * The leading dimension of the array A. LDA >= max(1,N). 54: * 55: * AF (input) DOUBLE PRECISION array, dimension (LDAF,N) 56: * The triangular factor U or L from the Cholesky factorization 57: * A = U**T*U or A = L*L**T, as computed by DPOTRF. 58: * 59: * LDAF (input) INTEGER 60: * The leading dimension of the array AF. LDAF >= max(1,N). 61: * 62: * CMODE (input) INTEGER 63: * Determines op2(C) in the formula op(A) * op2(C) as follows: 64: * CMODE = 1 op2(C) = C 65: * CMODE = 0 op2(C) = I 66: * CMODE = -1 op2(C) = inv(C) 67: * 68: * C (input) DOUBLE PRECISION array, dimension (N) 69: * The vector C in the formula op(A) * op2(C). 70: * 71: * INFO (output) INTEGER 72: * = 0: Successful exit. 73: * i > 0: The ith argument is invalid. 74: * 75: * WORK (input) DOUBLE PRECISION array, dimension (3*N). 76: * Workspace. 77: * 78: * IWORK (input) INTEGER array, dimension (N). 79: * Workspace. 80: * 81: * ===================================================================== 82: * 83: * .. Local Scalars .. 84: INTEGER KASE, I, J 85: DOUBLE PRECISION AINVNM, TMP 86: LOGICAL UP 87: * .. 88: * .. Array Arguments .. 89: INTEGER ISAVE( 3 ) 90: * .. 91: * .. External Functions .. 92: LOGICAL LSAME 93: INTEGER IDAMAX 94: EXTERNAL LSAME, IDAMAX 95: * .. 96: * .. External Subroutines .. 97: EXTERNAL DLACN2, DPOTRS, XERBLA 98: * .. 99: * .. Intrinsic Functions .. 100: INTRINSIC ABS, MAX 101: * .. 102: * .. Executable Statements .. 103: * 104: DLA_PORCOND = 0.0D+0 105: * 106: INFO = 0 107: IF( N.LT.0 ) THEN 108: INFO = -2 109: END IF 110: IF( INFO.NE.0 ) THEN 111: CALL XERBLA( 'DLA_PORCOND', -INFO ) 112: RETURN 113: END IF 114: 115: IF( N.EQ.0 ) THEN 116: DLA_PORCOND = 1.0D+0 117: RETURN 118: END IF 119: UP = .FALSE. 120: IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE. 121: * 122: * Compute the equilibration matrix R such that 123: * inv(R)*A*C has unit 1-norm. 124: * 125: IF ( UP ) THEN 126: DO I = 1, N 127: TMP = 0.0D+0 128: IF ( CMODE .EQ. 1 ) THEN 129: DO J = 1, I 130: TMP = TMP + ABS( A( J, I ) * C( J ) ) 131: END DO 132: DO J = I+1, N 133: TMP = TMP + ABS( A( I, J ) * C( J ) ) 134: END DO 135: ELSE IF ( CMODE .EQ. 0 ) THEN 136: DO J = 1, I 137: TMP = TMP + ABS( A( J, I ) ) 138: END DO 139: DO J = I+1, N 140: TMP = TMP + ABS( A( I, J ) ) 141: END DO 142: ELSE 143: DO J = 1, I 144: TMP = TMP + ABS( A( J ,I ) / C( J ) ) 145: END DO 146: DO J = I+1, N 147: TMP = TMP + ABS( A( I, J ) / C( J ) ) 148: END DO 149: END IF 150: WORK( 2*N+I ) = TMP 151: END DO 152: ELSE 153: DO I = 1, N 154: TMP = 0.0D+0 155: IF ( CMODE .EQ. 1 ) THEN 156: DO J = 1, I 157: TMP = TMP + ABS( A( I, J ) * C( J ) ) 158: END DO 159: DO J = I+1, N 160: TMP = TMP + ABS( A( J, I ) * C( J ) ) 161: END DO 162: ELSE IF ( CMODE .EQ. 0 ) THEN 163: DO J = 1, I 164: TMP = TMP + ABS( A( I, J ) ) 165: END DO 166: DO J = I+1, N 167: TMP = TMP + ABS( A( J, I ) ) 168: END DO 169: ELSE 170: DO J = 1, I 171: TMP = TMP + ABS( A( I, J ) / C( J ) ) 172: END DO 173: DO J = I+1, N 174: TMP = TMP + ABS( A( J, I ) / C( J ) ) 175: END DO 176: END IF 177: WORK( 2*N+I ) = TMP 178: END DO 179: ENDIF 180: * 181: * Estimate the norm of inv(op(A)). 182: * 183: AINVNM = 0.0D+0 184: 185: KASE = 0 186: 10 CONTINUE 187: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 188: IF( KASE.NE.0 ) THEN 189: IF( KASE.EQ.2 ) THEN 190: * 191: * Multiply by R. 192: * 193: DO I = 1, N 194: WORK( I ) = WORK( I ) * WORK( 2*N+I ) 195: END DO 196: 197: IF (UP) THEN 198: CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO ) 199: ELSE 200: CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO ) 201: ENDIF 202: * 203: * Multiply by inv(C). 204: * 205: IF ( CMODE .EQ. 1 ) THEN 206: DO I = 1, N 207: WORK( I ) = WORK( I ) / C( I ) 208: END DO 209: ELSE IF ( CMODE .EQ. -1 ) THEN 210: DO I = 1, N 211: WORK( I ) = WORK( I ) * C( I ) 212: END DO 213: END IF 214: ELSE 215: * 216: * Multiply by inv(C'). 217: * 218: IF ( CMODE .EQ. 1 ) THEN 219: DO I = 1, N 220: WORK( I ) = WORK( I ) / C( I ) 221: END DO 222: ELSE IF ( CMODE .EQ. -1 ) THEN 223: DO I = 1, N 224: WORK( I ) = WORK( I ) * C( I ) 225: END DO 226: END IF 227: 228: IF ( UP ) THEN 229: CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO ) 230: ELSE 231: CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO ) 232: ENDIF 233: * 234: * Multiply by R. 235: * 236: DO I = 1, N 237: WORK( I ) = WORK( I ) * WORK( 2*N+I ) 238: END DO 239: END IF 240: GO TO 10 241: END IF 242: * 243: * Compute the estimate of the reciprocal condition number. 244: * 245: IF( AINVNM .NE. 0.0D+0 ) 246: $ DLA_PORCOND = ( 1.0D+0 / AINVNM ) 247: * 248: RETURN 249: * 250: END