![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, 2: $ INFO ) 3: * 4: * -- LAPACK routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH. 10: * 11: * .. Scalar Arguments .. 12: CHARACTER UPLO 13: INTEGER INFO, LDA, N 14: DOUBLE PRECISION ANORM, RCOND 15: * .. 16: * .. Array Arguments .. 17: INTEGER IWORK( * ) 18: DOUBLE PRECISION A( LDA, * ), WORK( * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * DPOCON estimates the reciprocal of the condition number (in the 25: * 1-norm) of a real symmetric positive definite matrix using the 26: * Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF. 27: * 28: * An estimate is obtained for norm(inv(A)), and the reciprocal of the 29: * condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). 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 order of the matrix A. N >= 0. 40: * 41: * A (input) DOUBLE PRECISION array, dimension (LDA,N) 42: * The triangular factor U or L from the Cholesky factorization 43: * A = U**T*U or A = L*L**T, as computed by DPOTRF. 44: * 45: * LDA (input) INTEGER 46: * The leading dimension of the array A. LDA >= max(1,N). 47: * 48: * ANORM (input) DOUBLE PRECISION 49: * The 1-norm (or infinity-norm) of the symmetric matrix A. 50: * 51: * RCOND (output) DOUBLE PRECISION 52: * The reciprocal of the condition number of the matrix A, 53: * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an 54: * estimate of the 1-norm of inv(A) computed in this routine. 55: * 56: * WORK (workspace) DOUBLE PRECISION array, dimension (3*N) 57: * 58: * IWORK (workspace) INTEGER array, dimension (N) 59: * 60: * INFO (output) INTEGER 61: * = 0: successful exit 62: * < 0: if INFO = -i, the i-th argument had an illegal value 63: * 64: * ===================================================================== 65: * 66: * .. Parameters .. 67: DOUBLE PRECISION ONE, ZERO 68: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 69: * .. 70: * .. Local Scalars .. 71: LOGICAL UPPER 72: CHARACTER NORMIN 73: INTEGER IX, KASE 74: DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM 75: * .. 76: * .. Local Arrays .. 77: INTEGER ISAVE( 3 ) 78: * .. 79: * .. External Functions .. 80: LOGICAL LSAME 81: INTEGER IDAMAX 82: DOUBLE PRECISION DLAMCH 83: EXTERNAL LSAME, IDAMAX, DLAMCH 84: * .. 85: * .. External Subroutines .. 86: EXTERNAL DLACN2, DLATRS, DRSCL, XERBLA 87: * .. 88: * .. Intrinsic Functions .. 89: INTRINSIC ABS, MAX 90: * .. 91: * .. Executable Statements .. 92: * 93: * Test the input parameters. 94: * 95: INFO = 0 96: UPPER = LSAME( UPLO, 'U' ) 97: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 98: INFO = -1 99: ELSE IF( N.LT.0 ) THEN 100: INFO = -2 101: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 102: INFO = -4 103: ELSE IF( ANORM.LT.ZERO ) THEN 104: INFO = -5 105: END IF 106: IF( INFO.NE.0 ) THEN 107: CALL XERBLA( 'DPOCON', -INFO ) 108: RETURN 109: END IF 110: * 111: * Quick return if possible 112: * 113: RCOND = ZERO 114: IF( N.EQ.0 ) THEN 115: RCOND = ONE 116: RETURN 117: ELSE IF( ANORM.EQ.ZERO ) THEN 118: RETURN 119: END IF 120: * 121: SMLNUM = DLAMCH( 'Safe minimum' ) 122: * 123: * Estimate the 1-norm of inv(A). 124: * 125: KASE = 0 126: NORMIN = 'N' 127: 10 CONTINUE 128: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 129: IF( KASE.NE.0 ) THEN 130: IF( UPPER ) THEN 131: * 132: * Multiply by inv(U'). 133: * 134: CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A, 135: $ LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO ) 136: NORMIN = 'Y' 137: * 138: * Multiply by inv(U). 139: * 140: CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 141: $ A, LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO ) 142: ELSE 143: * 144: * Multiply by inv(L). 145: * 146: CALL DLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, 147: $ A, LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO ) 148: NORMIN = 'Y' 149: * 150: * Multiply by inv(L'). 151: * 152: CALL DLATRS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N, A, 153: $ LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO ) 154: END IF 155: * 156: * Multiply by 1/SCALE if doing so will not cause overflow. 157: * 158: SCALE = SCALEL*SCALEU 159: IF( SCALE.NE.ONE ) THEN 160: IX = IDAMAX( N, WORK, 1 ) 161: IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 162: $ GO TO 20 163: CALL DRSCL( N, SCALE, WORK, 1 ) 164: END IF 165: GO TO 10 166: END IF 167: * 168: * Compute the estimate of the reciprocal condition number. 169: * 170: IF( AINVNM.NE.ZERO ) 171: $ RCOND = ( ONE / AINVNM ) / ANORM 172: * 173: 20 CONTINUE 174: RETURN 175: * 176: * End of DPOCON 177: * 178: END