![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZTPCON( NORM, UPLO, DIAG, N, AP, RCOND, WORK, RWORK, 2: $ 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, N 14: DOUBLE PRECISION RCOND 15: * .. 16: * .. Array Arguments .. 17: DOUBLE PRECISION RWORK( * ) 18: COMPLEX*16 AP( * ), WORK( * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * ZTPCON estimates the reciprocal of the condition number of a packed 25: * triangular 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: * AP (input) COMPLEX*16 array, dimension (N*(N+1)/2) 53: * The upper or lower triangular matrix A, packed columnwise in 54: * a linear array. The j-th column of A is stored in the array 55: * AP as follows: 56: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 57: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 58: * If DIAG = 'U', the diagonal elements of A are not referenced 59: * and are assumed to be 1. 60: * 61: * RCOND (output) DOUBLE PRECISION 62: * The reciprocal of the condition number of the matrix A, 63: * computed as RCOND = 1/(norm(A) * norm(inv(A))). 64: * 65: * WORK (workspace) COMPLEX*16 array, dimension (2*N) 66: * 67: * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 68: * 69: * INFO (output) INTEGER 70: * = 0: successful exit 71: * < 0: if INFO = -i, the i-th argument had an illegal value 72: * 73: * ===================================================================== 74: * 75: * .. Parameters .. 76: DOUBLE PRECISION ONE, ZERO 77: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 78: * .. 79: * .. Local Scalars .. 80: LOGICAL NOUNIT, ONENRM, UPPER 81: CHARACTER NORMIN 82: INTEGER IX, KASE, KASE1 83: DOUBLE PRECISION AINVNM, ANORM, SCALE, SMLNUM, XNORM 84: COMPLEX*16 ZDUM 85: * .. 86: * .. Local Arrays .. 87: INTEGER ISAVE( 3 ) 88: * .. 89: * .. External Functions .. 90: LOGICAL LSAME 91: INTEGER IZAMAX 92: DOUBLE PRECISION DLAMCH, ZLANTP 93: EXTERNAL LSAME, IZAMAX, DLAMCH, ZLANTP 94: * .. 95: * .. External Subroutines .. 96: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLATPS 97: * .. 98: * .. Intrinsic Functions .. 99: INTRINSIC ABS, DBLE, DIMAG, MAX 100: * .. 101: * .. Statement Functions .. 102: DOUBLE PRECISION CABS1 103: * .. 104: * .. Statement Function definitions .. 105: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 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: END IF 125: IF( INFO.NE.0 ) THEN 126: CALL XERBLA( 'ZTPCON', -INFO ) 127: RETURN 128: END IF 129: * 130: * Quick return if possible 131: * 132: IF( N.EQ.0 ) THEN 133: RCOND = ONE 134: RETURN 135: END IF 136: * 137: RCOND = ZERO 138: SMLNUM = DLAMCH( 'Safe minimum' )*DBLE( MAX( 1, N ) ) 139: * 140: * Compute the norm of the triangular matrix A. 141: * 142: ANORM = ZLANTP( NORM, UPLO, DIAG, N, AP, RWORK ) 143: * 144: * Continue only if ANORM > 0. 145: * 146: IF( ANORM.GT.ZERO ) THEN 147: * 148: * Estimate the norm of the inverse of A. 149: * 150: AINVNM = ZERO 151: NORMIN = 'N' 152: IF( ONENRM ) THEN 153: KASE1 = 1 154: ELSE 155: KASE1 = 2 156: END IF 157: KASE = 0 158: 10 CONTINUE 159: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 160: IF( KASE.NE.0 ) THEN 161: IF( KASE.EQ.KASE1 ) THEN 162: * 163: * Multiply by inv(A). 164: * 165: CALL ZLATPS( UPLO, 'No transpose', DIAG, NORMIN, N, AP, 166: $ WORK, SCALE, RWORK, INFO ) 167: ELSE 168: * 169: * Multiply by inv(A'). 170: * 171: CALL ZLATPS( UPLO, 'Conjugate transpose', DIAG, NORMIN, 172: $ N, AP, WORK, SCALE, RWORK, INFO ) 173: END IF 174: NORMIN = 'Y' 175: * 176: * Multiply by 1/SCALE if doing so will not cause overflow. 177: * 178: IF( SCALE.NE.ONE ) THEN 179: IX = IZAMAX( N, WORK, 1 ) 180: XNORM = CABS1( WORK( IX ) ) 181: IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO ) 182: $ GO TO 20 183: CALL ZDRSCL( N, SCALE, WORK, 1 ) 184: END IF 185: GO TO 10 186: END IF 187: * 188: * Compute the estimate of the reciprocal condition number. 189: * 190: IF( AINVNM.NE.ZERO ) 191: $ RCOND = ( ONE / ANORM ) / AINVNM 192: END IF 193: * 194: 20 CONTINUE 195: RETURN 196: * 197: * End of ZTPCON 198: * 199: END