![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZTBCON( NORM, UPLO, DIAG, N, KD, AB, LDAB, RCOND, WORK, 2: $ RWORK, 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 ZLACN2 in place of ZLACON, 10 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: DOUBLE PRECISION RWORK( * ) 18: COMPLEX*16 AB( LDAB, * ), WORK( * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * ZTBCON 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) COMPLEX*16 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) COMPLEX*16 array, dimension (2*N) 73: * 74: * RWORK (workspace) DOUBLE PRECISION 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: COMPLEX*16 ZDUM 92: * .. 93: * .. Local Arrays .. 94: INTEGER ISAVE( 3 ) 95: * .. 96: * .. External Functions .. 97: LOGICAL LSAME 98: INTEGER IZAMAX 99: DOUBLE PRECISION DLAMCH, ZLANTB 100: EXTERNAL LSAME, IZAMAX, DLAMCH, ZLANTB 101: * .. 102: * .. External Subroutines .. 103: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLATBS 104: * .. 105: * .. Intrinsic Functions .. 106: INTRINSIC ABS, DBLE, DIMAG, MAX 107: * .. 108: * .. Statement Functions .. 109: DOUBLE PRECISION CABS1 110: * .. 111: * .. Statement Function definitions .. 112: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 113: * .. 114: * .. Executable Statements .. 115: * 116: * Test the input parameters. 117: * 118: INFO = 0 119: UPPER = LSAME( UPLO, 'U' ) 120: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' ) 121: NOUNIT = LSAME( DIAG, 'N' ) 122: * 123: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN 124: INFO = -1 125: ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 126: INFO = -2 127: ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 128: INFO = -3 129: ELSE IF( N.LT.0 ) THEN 130: INFO = -4 131: ELSE IF( KD.LT.0 ) THEN 132: INFO = -5 133: ELSE IF( LDAB.LT.KD+1 ) THEN 134: INFO = -7 135: END IF 136: IF( INFO.NE.0 ) THEN 137: CALL XERBLA( 'ZTBCON', -INFO ) 138: RETURN 139: END IF 140: * 141: * Quick return if possible 142: * 143: IF( N.EQ.0 ) THEN 144: RCOND = ONE 145: RETURN 146: END IF 147: * 148: RCOND = ZERO 149: SMLNUM = DLAMCH( 'Safe minimum' )*DBLE( MAX( N, 1 ) ) 150: * 151: * Compute the 1-norm of the triangular matrix A or A'. 152: * 153: ANORM = ZLANTB( NORM, UPLO, DIAG, N, KD, AB, LDAB, RWORK ) 154: * 155: * Continue only if ANORM > 0. 156: * 157: IF( ANORM.GT.ZERO ) THEN 158: * 159: * Estimate the 1-norm of the inverse of A. 160: * 161: AINVNM = ZERO 162: NORMIN = 'N' 163: IF( ONENRM ) THEN 164: KASE1 = 1 165: ELSE 166: KASE1 = 2 167: END IF 168: KASE = 0 169: 10 CONTINUE 170: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 171: IF( KASE.NE.0 ) THEN 172: IF( KASE.EQ.KASE1 ) THEN 173: * 174: * Multiply by inv(A). 175: * 176: CALL ZLATBS( UPLO, 'No transpose', DIAG, NORMIN, N, KD, 177: $ AB, LDAB, WORK, SCALE, RWORK, INFO ) 178: ELSE 179: * 180: * Multiply by inv(A'). 181: * 182: CALL ZLATBS( UPLO, 'Conjugate transpose', DIAG, NORMIN, 183: $ N, KD, AB, LDAB, WORK, SCALE, RWORK, INFO ) 184: END IF 185: NORMIN = 'Y' 186: * 187: * Multiply by 1/SCALE if doing so will not cause overflow. 188: * 189: IF( SCALE.NE.ONE ) THEN 190: IX = IZAMAX( N, WORK, 1 ) 191: XNORM = CABS1( WORK( IX ) ) 192: IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO ) 193: $ GO TO 20 194: CALL ZDRSCL( N, SCALE, WORK, 1 ) 195: END IF 196: GO TO 10 197: END IF 198: * 199: * Compute the estimate of the reciprocal condition number. 200: * 201: IF( AINVNM.NE.ZERO ) 202: $ RCOND = ( ONE / ANORM ) / AINVNM 203: END IF 204: * 205: 20 CONTINUE 206: RETURN 207: * 208: * End of ZTBCON 209: * 210: END