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