![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DPOTRF( 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: * DPOTRF computes the Cholesky factorization of a real symmetric 20: * positive definite matrix A. 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: * This is the block version of the algorithm, calling Level 3 BLAS. 28: * 29: * Arguments 30: * ========= 31: * 32: * UPLO (input) CHARACTER*1 33: * = 'U': Upper triangle of A is stored; 34: * = 'L': Lower triangle of A is stored. 35: * 36: * N (input) INTEGER 37: * The order of the matrix A. N >= 0. 38: * 39: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 40: * On entry, the symmetric matrix A. If UPLO = 'U', the leading 41: * N-by-N upper triangular part of A contains the upper 42: * triangular part of the matrix A, and the strictly lower 43: * triangular part of A is not referenced. If UPLO = 'L', the 44: * leading N-by-N lower triangular part of A contains the lower 45: * triangular part of the matrix A, and the strictly upper 46: * triangular part of A is not referenced. 47: * 48: * On exit, if INFO = 0, the factor U or L from the Cholesky 49: * factorization A = U**T*U or A = L*L**T. 50: * 51: * LDA (input) INTEGER 52: * The leading dimension of the array A. LDA >= max(1,N). 53: * 54: * INFO (output) INTEGER 55: * = 0: successful exit 56: * < 0: if INFO = -i, the i-th argument had an illegal value 57: * > 0: if INFO = i, the leading minor of order i is not 58: * positive definite, and the factorization could not be 59: * completed. 60: * 61: * ===================================================================== 62: * 63: * .. Parameters .. 64: DOUBLE PRECISION ONE 65: PARAMETER ( ONE = 1.0D+0 ) 66: * .. 67: * .. Local Scalars .. 68: LOGICAL UPPER 69: INTEGER J, JB, NB 70: * .. 71: * .. External Functions .. 72: LOGICAL LSAME 73: INTEGER ILAENV 74: EXTERNAL LSAME, ILAENV 75: * .. 76: * .. External Subroutines .. 77: EXTERNAL DGEMM, DPOTF2, DSYRK, DTRSM, XERBLA 78: * .. 79: * .. Intrinsic Functions .. 80: INTRINSIC MAX, MIN 81: * .. 82: * .. Executable Statements .. 83: * 84: * Test the input parameters. 85: * 86: INFO = 0 87: UPPER = LSAME( UPLO, 'U' ) 88: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 89: INFO = -1 90: ELSE IF( N.LT.0 ) THEN 91: INFO = -2 92: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 93: INFO = -4 94: END IF 95: IF( INFO.NE.0 ) THEN 96: CALL XERBLA( 'DPOTRF', -INFO ) 97: RETURN 98: END IF 99: * 100: * Quick return if possible 101: * 102: IF( N.EQ.0 ) 103: $ RETURN 104: * 105: * Determine the block size for this environment. 106: * 107: NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 ) 108: IF( NB.LE.1 .OR. NB.GE.N ) THEN 109: * 110: * Use unblocked code. 111: * 112: CALL DPOTF2( UPLO, N, A, LDA, INFO ) 113: ELSE 114: * 115: * Use blocked code. 116: * 117: IF( UPPER ) THEN 118: * 119: * Compute the Cholesky factorization A = U'*U. 120: * 121: DO 10 J = 1, N, NB 122: * 123: * Update and factorize the current diagonal block and test 124: * for non-positive-definiteness. 125: * 126: JB = MIN( NB, N-J+1 ) 127: CALL DSYRK( 'Upper', 'Transpose', JB, J-1, -ONE, 128: $ A( 1, J ), LDA, ONE, A( J, J ), LDA ) 129: CALL DPOTF2( 'Upper', JB, A( J, J ), LDA, INFO ) 130: IF( INFO.NE.0 ) 131: $ GO TO 30 132: IF( J+JB.LE.N ) THEN 133: * 134: * Compute the current block row. 135: * 136: CALL DGEMM( 'Transpose', 'No transpose', JB, N-J-JB+1, 137: $ J-1, -ONE, A( 1, J ), LDA, A( 1, J+JB ), 138: $ LDA, ONE, A( J, J+JB ), LDA ) 139: CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', 140: $ JB, N-J-JB+1, ONE, A( J, J ), LDA, 141: $ A( J, J+JB ), LDA ) 142: END IF 143: 10 CONTINUE 144: * 145: ELSE 146: * 147: * Compute the Cholesky factorization A = L*L'. 148: * 149: DO 20 J = 1, N, NB 150: * 151: * Update and factorize the current diagonal block and test 152: * for non-positive-definiteness. 153: * 154: JB = MIN( NB, N-J+1 ) 155: CALL DSYRK( 'Lower', 'No transpose', JB, J-1, -ONE, 156: $ A( J, 1 ), LDA, ONE, A( J, J ), LDA ) 157: CALL DPOTF2( 'Lower', JB, A( J, J ), LDA, INFO ) 158: IF( INFO.NE.0 ) 159: $ GO TO 30 160: IF( J+JB.LE.N ) THEN 161: * 162: * Compute the current block column. 163: * 164: CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB, 165: $ J-1, -ONE, A( J+JB, 1 ), LDA, A( J, 1 ), 166: $ LDA, ONE, A( J+JB, J ), LDA ) 167: CALL DTRSM( 'Right', 'Lower', 'Transpose', 'Non-unit', 168: $ N-J-JB+1, JB, ONE, A( J, J ), LDA, 169: $ A( J+JB, J ), LDA ) 170: END IF 171: 20 CONTINUE 172: END IF 173: END IF 174: GO TO 40 175: * 176: 30 CONTINUE 177: INFO = INFO + J - 1 178: * 179: 40 CONTINUE 180: RETURN 181: * 182: * End of DPOTRF 183: * 184: END