![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DPOTF2( 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: DOUBLE PRECISION A( LDA, * ) 14: * .. 15: * 16: * Purpose 17: * ======= 18: * 19: * DPOTF2 computes the Cholesky factorization of a real symmetric 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: * symmetric 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) DOUBLE PRECISION array, dimension (LDA,N) 42: * On entry, the symmetric 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: * .. 69: * .. Local Scalars .. 70: LOGICAL UPPER 71: INTEGER J 72: DOUBLE PRECISION AJJ 73: * .. 74: * .. External Functions .. 75: LOGICAL LSAME, DISNAN 76: DOUBLE PRECISION DDOT 77: EXTERNAL LSAME, DDOT, DISNAN 78: * .. 79: * .. External Subroutines .. 80: EXTERNAL DGEMV, DSCAL, XERBLA 81: * .. 82: * .. Intrinsic Functions .. 83: INTRINSIC MAX, SQRT 84: * .. 85: * .. Executable Statements .. 86: * 87: * Test the input parameters. 88: * 89: INFO = 0 90: UPPER = LSAME( UPLO, 'U' ) 91: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 92: INFO = -1 93: ELSE IF( N.LT.0 ) THEN 94: INFO = -2 95: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 96: INFO = -4 97: END IF 98: IF( INFO.NE.0 ) THEN 99: CALL XERBLA( 'DPOTF2', -INFO ) 100: RETURN 101: END IF 102: * 103: * Quick return if possible 104: * 105: IF( N.EQ.0 ) 106: $ RETURN 107: * 108: IF( UPPER ) THEN 109: * 110: * Compute the Cholesky factorization A = U'*U. 111: * 112: DO 10 J = 1, N 113: * 114: * Compute U(J,J) and test for non-positive-definiteness. 115: * 116: AJJ = A( J, J ) - DDOT( J-1, A( 1, J ), 1, A( 1, J ), 1 ) 117: IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN 118: A( J, J ) = AJJ 119: GO TO 30 120: END IF 121: AJJ = SQRT( AJJ ) 122: A( J, J ) = AJJ 123: * 124: * Compute elements J+1:N of row J. 125: * 126: IF( J.LT.N ) THEN 127: CALL DGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ), 128: $ LDA, A( 1, J ), 1, ONE, A( J, J+1 ), LDA ) 129: CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) 130: END IF 131: 10 CONTINUE 132: ELSE 133: * 134: * Compute the Cholesky factorization A = L*L'. 135: * 136: DO 20 J = 1, N 137: * 138: * Compute L(J,J) and test for non-positive-definiteness. 139: * 140: AJJ = A( J, J ) - DDOT( J-1, A( J, 1 ), LDA, A( J, 1 ), 141: $ LDA ) 142: IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN 143: A( J, J ) = AJJ 144: GO TO 30 145: END IF 146: AJJ = SQRT( AJJ ) 147: A( J, J ) = AJJ 148: * 149: * Compute elements J+1:N of column J. 150: * 151: IF( J.LT.N ) THEN 152: CALL DGEMV( 'No transpose', N-J, J-1, -ONE, A( J+1, 1 ), 153: $ LDA, A( J, 1 ), LDA, ONE, A( J+1, J ), 1 ) 154: CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) 155: END IF 156: 20 CONTINUE 157: END IF 158: GO TO 40 159: * 160: 30 CONTINUE 161: INFO = J 162: * 163: 40 CONTINUE 164: RETURN 165: * 166: * End of DPOTF2 167: * 168: END