![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZPPTRF( UPLO, N, AP, 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, N 11: * .. 12: * .. Array Arguments .. 13: COMPLEX*16 AP( * ) 14: * .. 15: * 16: * Purpose 17: * ======= 18: * 19: * ZPPTRF computes the Cholesky factorization of a complex Hermitian 20: * positive definite matrix A stored in packed format. 21: * 22: * The factorization has the form 23: * A = U**H * U, if UPLO = 'U', or 24: * A = L * L**H, if UPLO = 'L', 25: * where U is an upper triangular matrix and L is lower triangular. 26: * 27: * Arguments 28: * ========= 29: * 30: * UPLO (input) CHARACTER*1 31: * = 'U': Upper triangle of A is stored; 32: * = 'L': Lower triangle of A is stored. 33: * 34: * N (input) INTEGER 35: * The order of the matrix A. N >= 0. 36: * 37: * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) 38: * On entry, the upper or lower triangle of the Hermitian matrix 39: * A, packed columnwise in a linear array. The j-th column of A 40: * is stored in the array AP as follows: 41: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 42: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 43: * See below for further details. 44: * 45: * On exit, if INFO = 0, the triangular factor U or L from the 46: * Cholesky factorization A = U**H*U or A = L*L**H, in the same 47: * storage format as A. 48: * 49: * INFO (output) INTEGER 50: * = 0: successful exit 51: * < 0: if INFO = -i, the i-th argument had an illegal value 52: * > 0: if INFO = i, the leading minor of order i is not 53: * positive definite, and the factorization could not be 54: * completed. 55: * 56: * Further Details 57: * =============== 58: * 59: * The packed storage scheme is illustrated by the following example 60: * when N = 4, UPLO = 'U': 61: * 62: * Two-dimensional storage of the Hermitian matrix A: 63: * 64: * a11 a12 a13 a14 65: * a22 a23 a24 66: * a33 a34 (aij = conjg(aji)) 67: * a44 68: * 69: * Packed storage of the upper triangle of A: 70: * 71: * AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] 72: * 73: * ===================================================================== 74: * 75: * .. Parameters .. 76: DOUBLE PRECISION ZERO, ONE 77: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 78: * .. 79: * .. Local Scalars .. 80: LOGICAL UPPER 81: INTEGER J, JC, JJ 82: DOUBLE PRECISION AJJ 83: * .. 84: * .. External Functions .. 85: LOGICAL LSAME 86: COMPLEX*16 ZDOTC 87: EXTERNAL LSAME, ZDOTC 88: * .. 89: * .. External Subroutines .. 90: EXTERNAL XERBLA, ZDSCAL, ZHPR, ZTPSV 91: * .. 92: * .. Intrinsic Functions .. 93: INTRINSIC DBLE, SQRT 94: * .. 95: * .. Executable Statements .. 96: * 97: * Test the input parameters. 98: * 99: INFO = 0 100: UPPER = LSAME( UPLO, 'U' ) 101: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 102: INFO = -1 103: ELSE IF( N.LT.0 ) THEN 104: INFO = -2 105: END IF 106: IF( INFO.NE.0 ) THEN 107: CALL XERBLA( 'ZPPTRF', -INFO ) 108: RETURN 109: END IF 110: * 111: * Quick return if possible 112: * 113: IF( N.EQ.0 ) 114: $ RETURN 115: * 116: IF( UPPER ) THEN 117: * 118: * Compute the Cholesky factorization A = U'*U. 119: * 120: JJ = 0 121: DO 10 J = 1, N 122: JC = JJ + 1 123: JJ = JJ + J 124: * 125: * Compute elements 1:J-1 of column J. 126: * 127: IF( J.GT.1 ) 128: $ CALL ZTPSV( 'Upper', 'Conjugate transpose', 'Non-unit', 129: $ J-1, AP, AP( JC ), 1 ) 130: * 131: * Compute U(J,J) and test for non-positive-definiteness. 132: * 133: AJJ = DBLE( AP( JJ ) ) - ZDOTC( J-1, AP( JC ), 1, AP( JC ), 134: $ 1 ) 135: IF( AJJ.LE.ZERO ) THEN 136: AP( JJ ) = AJJ 137: GO TO 30 138: END IF 139: AP( JJ ) = SQRT( AJJ ) 140: 10 CONTINUE 141: ELSE 142: * 143: * Compute the Cholesky factorization A = L*L'. 144: * 145: JJ = 1 146: DO 20 J = 1, N 147: * 148: * Compute L(J,J) and test for non-positive-definiteness. 149: * 150: AJJ = DBLE( AP( JJ ) ) 151: IF( AJJ.LE.ZERO ) THEN 152: AP( JJ ) = AJJ 153: GO TO 30 154: END IF 155: AJJ = SQRT( AJJ ) 156: AP( JJ ) = AJJ 157: * 158: * Compute elements J+1:N of column J and update the trailing 159: * submatrix. 160: * 161: IF( J.LT.N ) THEN 162: CALL ZDSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 ) 163: CALL ZHPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1, 164: $ AP( JJ+N-J+1 ) ) 165: JJ = JJ + N - J + 1 166: END IF 167: 20 CONTINUE 168: END IF 169: GO TO 40 170: * 171: 30 CONTINUE 172: INFO = J 173: * 174: 40 CONTINUE 175: RETURN 176: * 177: * End of ZPPTRF 178: * 179: END