![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.2.2.
1: SUBROUTINE DSPCON( UPLO, N, AP, IPIV, ANORM, RCOND, WORK, IWORK, 2: $ INFO ) 3: * 4: * -- LAPACK routine (version 3.2.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: * June 2010 8: * 9: * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH. 10: * 11: * .. Scalar Arguments .. 12: CHARACTER UPLO 13: INTEGER INFO, N 14: DOUBLE PRECISION ANORM, RCOND 15: * .. 16: * .. Array Arguments .. 17: INTEGER IPIV( * ), IWORK( * ) 18: DOUBLE PRECISION AP( * ), WORK( * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * DSPCON estimates the reciprocal of the condition number (in the 25: * 1-norm) of a real symmetric packed matrix A using the factorization 26: * A = U*D*U**T or A = L*D*L**T computed by DSPTRF. 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: * Specifies whether the details of the factorization are stored 36: * as an upper or lower triangular matrix. 37: * = 'U': Upper triangular, form is A = U*D*U**T; 38: * = 'L': Lower triangular, form is A = L*D*L**T. 39: * 40: * N (input) INTEGER 41: * The order of the matrix A. N >= 0. 42: * 43: * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2) 44: * The block diagonal matrix D and the multipliers used to 45: * obtain the factor U or L as computed by DSPTRF, stored as a 46: * packed triangular matrix. 47: * 48: * IPIV (input) INTEGER array, dimension (N) 49: * Details of the interchanges and the block structure of D 50: * as determined by DSPTRF. 51: * 52: * ANORM (input) DOUBLE PRECISION 53: * The 1-norm of the original matrix A. 54: * 55: * RCOND (output) DOUBLE PRECISION 56: * The reciprocal of the condition number of the matrix A, 57: * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an 58: * estimate of the 1-norm of inv(A) computed in this routine. 59: * 60: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N) 61: * 62: * IWORK (workspace) INTEGER array, dimension (N) 63: * 64: * INFO (output) INTEGER 65: * = 0: successful exit 66: * < 0: if INFO = -i, the i-th argument had an illegal value 67: * 68: * ===================================================================== 69: * 70: * .. Parameters .. 71: DOUBLE PRECISION ONE, ZERO 72: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 73: * .. 74: * .. Local Scalars .. 75: LOGICAL UPPER 76: INTEGER I, IP, KASE 77: DOUBLE PRECISION AINVNM 78: * .. 79: * .. Local Arrays .. 80: INTEGER ISAVE( 3 ) 81: * .. 82: * .. External Functions .. 83: LOGICAL LSAME 84: EXTERNAL LSAME 85: * .. 86: * .. External Subroutines .. 87: EXTERNAL DLACN2, DSPTRS, XERBLA 88: * .. 89: * .. Executable Statements .. 90: * 91: * Test the input parameters. 92: * 93: INFO = 0 94: UPPER = LSAME( UPLO, 'U' ) 95: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 96: INFO = -1 97: ELSE IF( N.LT.0 ) THEN 98: INFO = -2 99: ELSE IF( ANORM.LT.ZERO ) THEN 100: INFO = -5 101: END IF 102: IF( INFO.NE.0 ) THEN 103: CALL XERBLA( 'DSPCON', -INFO ) 104: RETURN 105: END IF 106: * 107: * Quick return if possible 108: * 109: RCOND = ZERO 110: IF( N.EQ.0 ) THEN 111: RCOND = ONE 112: RETURN 113: ELSE IF( ANORM.LE.ZERO ) THEN 114: RETURN 115: END IF 116: * 117: * Check that the diagonal matrix D is nonsingular. 118: * 119: IF( UPPER ) THEN 120: * 121: * Upper triangular storage: examine D from bottom to top 122: * 123: IP = N*( N+1 ) / 2 124: DO 10 I = N, 1, -1 125: IF( IPIV( I ).GT.0 .AND. AP( IP ).EQ.ZERO ) 126: $ RETURN 127: IP = IP - I 128: 10 CONTINUE 129: ELSE 130: * 131: * Lower triangular storage: examine D from top to bottom. 132: * 133: IP = 1 134: DO 20 I = 1, N 135: IF( IPIV( I ).GT.0 .AND. AP( IP ).EQ.ZERO ) 136: $ RETURN 137: IP = IP + N - I + 1 138: 20 CONTINUE 139: END IF 140: * 141: * Estimate the 1-norm of the inverse. 142: * 143: KASE = 0 144: 30 CONTINUE 145: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 146: IF( KASE.NE.0 ) THEN 147: * 148: * Multiply by inv(L*D*L') or inv(U*D*U'). 149: * 150: CALL DSPTRS( UPLO, N, 1, AP, IPIV, WORK, N, INFO ) 151: GO TO 30 152: END IF 153: * 154: * Compute the estimate of the reciprocal condition number. 155: * 156: IF( AINVNM.NE.ZERO ) 157: $ RCOND = ( ONE / AINVNM ) / ANORM 158: * 159: RETURN 160: * 161: * End of DSPCON 162: * 163: END