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