![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZPOTF2( UPLO, N, A, LDA, 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: * .. Scalar Arguments .. 9: CHARACTER UPLO 10: INTEGER INFO, LDA, N 11: * .. 12: * .. Array Arguments .. 13: COMPLEX*16 A( LDA, * ) 14: * .. 15: * 16: * Purpose 17: * ======= 18: * 19: * ZPOTF2 computes the Cholesky factorization of a complex Hermitian 20: * positive definite matrix A. 21: * 22: * The factorization has the form 23: * A = U' * U , if UPLO = 'U', or 24: * A = L * L', if UPLO = 'L', 25: * where U is an upper triangular matrix and L is lower triangular. 26: * 27: * This is the unblocked version of the algorithm, calling Level 2 BLAS. 28: * 29: * Arguments 30: * ========= 31: * 32: * UPLO (input) CHARACTER*1 33: * Specifies whether the upper or lower triangular part of the 34: * Hermitian matrix A is stored. 35: * = 'U': Upper triangular 36: * = 'L': Lower triangular 37: * 38: * N (input) INTEGER 39: * The order of the matrix A. N >= 0. 40: * 41: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 42: * On entry, the Hermitian matrix A. If UPLO = 'U', the leading 43: * n by n upper triangular part of A contains the upper 44: * triangular part of the matrix A, and the strictly lower 45: * triangular part of A is not referenced. If UPLO = 'L', the 46: * leading n by n lower triangular part of A contains the lower 47: * triangular part of the matrix A, and the strictly upper 48: * triangular part of A is not referenced. 49: * 50: * On exit, if INFO = 0, the factor U or L from the Cholesky 51: * factorization A = U'*U or A = L*L'. 52: * 53: * LDA (input) INTEGER 54: * The leading dimension of the array A. LDA >= max(1,N). 55: * 56: * INFO (output) INTEGER 57: * = 0: successful exit 58: * < 0: if INFO = -k, the k-th argument had an illegal value 59: * > 0: if INFO = k, the leading minor of order k is not 60: * positive definite, and the factorization could not be 61: * completed. 62: * 63: * ===================================================================== 64: * 65: * .. Parameters .. 66: DOUBLE PRECISION ONE, ZERO 67: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 68: COMPLEX*16 CONE 69: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 70: * .. 71: * .. Local Scalars .. 72: LOGICAL UPPER 73: INTEGER J 74: DOUBLE PRECISION AJJ 75: * .. 76: * .. External Functions .. 77: LOGICAL LSAME, DISNAN 78: COMPLEX*16 ZDOTC 79: EXTERNAL LSAME, ZDOTC, DISNAN 80: * .. 81: * .. External Subroutines .. 82: EXTERNAL XERBLA, ZDSCAL, ZGEMV, ZLACGV 83: * .. 84: * .. Intrinsic Functions .. 85: INTRINSIC DBLE, MAX, SQRT 86: * .. 87: * .. Executable Statements .. 88: * 89: * Test the input parameters. 90: * 91: INFO = 0 92: UPPER = LSAME( UPLO, 'U' ) 93: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 94: INFO = -1 95: ELSE IF( N.LT.0 ) THEN 96: INFO = -2 97: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 98: INFO = -4 99: END IF 100: IF( INFO.NE.0 ) THEN 101: CALL XERBLA( 'ZPOTF2', -INFO ) 102: RETURN 103: END IF 104: * 105: * Quick return if possible 106: * 107: IF( N.EQ.0 ) 108: $ RETURN 109: * 110: IF( UPPER ) THEN 111: * 112: * Compute the Cholesky factorization A = U'*U. 113: * 114: DO 10 J = 1, N 115: * 116: * Compute U(J,J) and test for non-positive-definiteness. 117: * 118: AJJ = DBLE( A( J, J ) ) - ZDOTC( J-1, A( 1, J ), 1, 119: $ A( 1, J ), 1 ) 120: IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN 121: A( J, J ) = AJJ 122: GO TO 30 123: END IF 124: AJJ = SQRT( AJJ ) 125: A( J, J ) = AJJ 126: * 127: * Compute elements J+1:N of row J. 128: * 129: IF( J.LT.N ) THEN 130: CALL ZLACGV( J-1, A( 1, J ), 1 ) 131: CALL ZGEMV( 'Transpose', J-1, N-J, -CONE, A( 1, J+1 ), 132: $ LDA, A( 1, J ), 1, CONE, A( J, J+1 ), LDA ) 133: CALL ZLACGV( J-1, A( 1, J ), 1 ) 134: CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) 135: END IF 136: 10 CONTINUE 137: ELSE 138: * 139: * Compute the Cholesky factorization A = L*L'. 140: * 141: DO 20 J = 1, N 142: * 143: * Compute L(J,J) and test for non-positive-definiteness. 144: * 145: AJJ = DBLE( A( J, J ) ) - ZDOTC( J-1, A( J, 1 ), LDA, 146: $ A( J, 1 ), LDA ) 147: IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN 148: A( J, J ) = AJJ 149: GO TO 30 150: END IF 151: AJJ = SQRT( AJJ ) 152: A( J, J ) = AJJ 153: * 154: * Compute elements J+1:N of column J. 155: * 156: IF( J.LT.N ) THEN 157: CALL ZLACGV( J-1, A( J, 1 ), LDA ) 158: CALL ZGEMV( 'No transpose', N-J, J-1, -CONE, A( J+1, 1 ), 159: $ LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 ) 160: CALL ZLACGV( J-1, A( J, 1 ), LDA ) 161: CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) 162: END IF 163: 20 CONTINUE 164: END IF 165: GO TO 40 166: * 167: 30 CONTINUE 168: INFO = J 169: * 170: 40 CONTINUE 171: RETURN 172: * 173: * End of ZPOTF2 174: * 175: END