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