![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZPPCON( UPLO, N, AP, ANORM, RCOND, WORK, RWORK, INFO ) 2: * 3: * -- LAPACK routine (version 3.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * November 2006 7: * 8: * Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH. 9: * 10: * .. Scalar Arguments .. 11: CHARACTER UPLO 12: INTEGER INFO, N 13: DOUBLE PRECISION ANORM, RCOND 14: * .. 15: * .. Array Arguments .. 16: DOUBLE PRECISION RWORK( * ) 17: COMPLEX*16 AP( * ), WORK( * ) 18: * .. 19: * 20: * Purpose 21: * ======= 22: * 23: * ZPPCON estimates the reciprocal of the condition number (in the 24: * 1-norm) of a complex Hermitian positive definite packed matrix using 25: * the Cholesky factorization A = U**H*U or A = L*L**H computed by 26: * ZPPTRF. 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 triangle of A is stored; 36: * = 'L': Lower triangle of A is stored. 37: * 38: * N (input) INTEGER 39: * The order of the matrix A. N >= 0. 40: * 41: * AP (input) COMPLEX*16 array, dimension (N*(N+1)/2) 42: * The triangular factor U or L from the Cholesky factorization 43: * A = U**H*U or A = L*L**H, packed columnwise in a linear 44: * array. The j-th column of U or L is stored in the array AP 45: * as follows: 46: * if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j; 47: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n. 48: * 49: * ANORM (input) DOUBLE PRECISION 50: * The 1-norm (or infinity-norm) of the Hermitian matrix A. 51: * 52: * RCOND (output) DOUBLE PRECISION 53: * The reciprocal of the condition number of the matrix A, 54: * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an 55: * estimate of the 1-norm of inv(A) computed in this routine. 56: * 57: * WORK (workspace) COMPLEX*16 array, dimension (2*N) 58: * 59: * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 60: * 61: * INFO (output) INTEGER 62: * = 0: successful exit 63: * < 0: if INFO = -i, the i-th argument had an illegal value 64: * 65: * ===================================================================== 66: * 67: * .. Parameters .. 68: DOUBLE PRECISION ONE, ZERO 69: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 70: * .. 71: * .. Local Scalars .. 72: LOGICAL UPPER 73: CHARACTER NORMIN 74: INTEGER IX, KASE 75: DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM 76: COMPLEX*16 ZDUM 77: * .. 78: * .. Local Arrays .. 79: INTEGER ISAVE( 3 ) 80: * .. 81: * .. External Functions .. 82: LOGICAL LSAME 83: INTEGER IZAMAX 84: DOUBLE PRECISION DLAMCH 85: EXTERNAL LSAME, IZAMAX, DLAMCH 86: * .. 87: * .. External Subroutines .. 88: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLATPS 89: * .. 90: * .. Intrinsic Functions .. 91: INTRINSIC ABS, DBLE, DIMAG 92: * .. 93: * .. Statement Functions .. 94: DOUBLE PRECISION CABS1 95: * .. 96: * .. Statement Function definitions .. 97: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 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( ANORM.LT.ZERO ) THEN 110: INFO = -4 111: END IF 112: IF( INFO.NE.0 ) THEN 113: CALL XERBLA( 'ZPPCON', -INFO ) 114: RETURN 115: END IF 116: * 117: * Quick return if possible 118: * 119: RCOND = ZERO 120: IF( N.EQ.0 ) THEN 121: RCOND = ONE 122: RETURN 123: ELSE IF( ANORM.EQ.ZERO ) THEN 124: RETURN 125: END IF 126: * 127: SMLNUM = DLAMCH( 'Safe minimum' ) 128: * 129: * Estimate the 1-norm of the inverse. 130: * 131: KASE = 0 132: NORMIN = 'N' 133: 10 CONTINUE 134: CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 135: IF( KASE.NE.0 ) THEN 136: IF( UPPER ) THEN 137: * 138: * Multiply by inv(U'). 139: * 140: CALL ZLATPS( 'Upper', 'Conjugate transpose', 'Non-unit', 141: $ NORMIN, N, AP, WORK, SCALEL, RWORK, INFO ) 142: NORMIN = 'Y' 143: * 144: * Multiply by inv(U). 145: * 146: CALL ZLATPS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 147: $ AP, WORK, SCALEU, RWORK, INFO ) 148: ELSE 149: * 150: * Multiply by inv(L). 151: * 152: CALL ZLATPS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, 153: $ AP, WORK, SCALEL, RWORK, INFO ) 154: NORMIN = 'Y' 155: * 156: * Multiply by inv(L'). 157: * 158: CALL ZLATPS( 'Lower', 'Conjugate transpose', 'Non-unit', 159: $ NORMIN, N, AP, WORK, SCALEU, RWORK, INFO ) 160: END IF 161: * 162: * Multiply by 1/SCALE if doing so will not cause overflow. 163: * 164: SCALE = SCALEL*SCALEU 165: IF( SCALE.NE.ONE ) THEN 166: IX = IZAMAX( N, WORK, 1 ) 167: IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 168: $ GO TO 20 169: CALL ZDRSCL( N, SCALE, WORK, 1 ) 170: END IF 171: GO TO 10 172: END IF 173: * 174: * Compute the estimate of the reciprocal condition number. 175: * 176: IF( AINVNM.NE.ZERO ) 177: $ RCOND = ( ONE / AINVNM ) / ANORM 178: * 179: 20 CONTINUE 180: RETURN 181: * 182: * End of ZPPCON 183: * 184: END