![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DPPTRF( 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: DOUBLE PRECISION AP( * ) 14: * .. 15: * 16: * Purpose 17: * ======= 18: * 19: * DPPTRF computes the Cholesky factorization of a real symmetric 20: * positive definite matrix A stored in packed format. 21: * 22: * The factorization has the form 23: * A = U**T * U, if UPLO = 'U', or 24: * A = L * L**T, 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) DOUBLE PRECISION array, dimension (N*(N+1)/2) 38: * On entry, the upper or lower triangle of the symmetric 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**T*U or A = L*L**T, 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 symmetric matrix A: 63: * 64: * a11 a12 a13 a14 65: * a22 a23 a24 66: * a33 a34 (aij = 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 ONE, ZERO 77: PARAMETER ( ONE = 1.0D+0, ZERO = 0.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: DOUBLE PRECISION DDOT 87: EXTERNAL LSAME, DDOT 88: * .. 89: * .. External Subroutines .. 90: EXTERNAL DSCAL, DSPR, DTPSV, XERBLA 91: * .. 92: * .. Intrinsic Functions .. 93: INTRINSIC 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( 'DPPTRF', -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 DTPSV( 'Upper', 'Transpose', 'Non-unit', J-1, AP, 129: $ AP( JC ), 1 ) 130: * 131: * Compute U(J,J) and test for non-positive-definiteness. 132: * 133: AJJ = AP( JJ ) - DDOT( J-1, AP( JC ), 1, AP( JC ), 1 ) 134: IF( AJJ.LE.ZERO ) THEN 135: AP( JJ ) = AJJ 136: GO TO 30 137: END IF 138: AP( JJ ) = SQRT( AJJ ) 139: 10 CONTINUE 140: ELSE 141: * 142: * Compute the Cholesky factorization A = L*L'. 143: * 144: JJ = 1 145: DO 20 J = 1, N 146: * 147: * Compute L(J,J) and test for non-positive-definiteness. 148: * 149: AJJ = AP( JJ ) 150: IF( AJJ.LE.ZERO ) THEN 151: AP( JJ ) = AJJ 152: GO TO 30 153: END IF 154: AJJ = SQRT( AJJ ) 155: AP( JJ ) = AJJ 156: * 157: * Compute elements J+1:N of column J and update the trailing 158: * submatrix. 159: * 160: IF( J.LT.N ) THEN 161: CALL DSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 ) 162: CALL DSPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1, 163: $ AP( JJ+N-J+1 ) ) 164: JJ = JJ + N - J + 1 165: END IF 166: 20 CONTINUE 167: END IF 168: GO TO 40 169: * 170: 30 CONTINUE 171: INFO = J 172: * 173: 40 CONTINUE 174: RETURN 175: * 176: * End of DPPTRF 177: * 178: END