LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DPOTRF ( UPLO, N, A, LDA, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.1) -- 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * March 2008 00006 * 00007 * .. Scalar Arguments .. 00008 CHARACTER UPLO 00009 INTEGER INFO, LDA, N 00010 * .. 00011 * .. Array Arguments .. 00012 DOUBLE PRECISION A( LDA, * ) 00013 * .. 00014 * 00015 * Purpose 00016 * ======= 00017 * 00018 * DPOTRF computes the Cholesky factorization of a real symmetric 00019 * positive definite matrix A. 00020 * 00021 * The factorization has the form 00022 * A = U**T * U, if UPLO = 'U', or 00023 * A = L * L**T, if UPLO = 'L', 00024 * where U is an upper triangular matrix and L is lower triangular. 00025 * 00026 * This is the right looking block version of the algorithm, calling Level 3 BLAS. 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * UPLO (input) CHARACTER*1 00032 * = 'U': Upper triangle of A is stored; 00033 * = 'L': Lower triangle of A is stored. 00034 * 00035 * N (input) INTEGER 00036 * The order of the matrix A. N >= 0. 00037 * 00038 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00039 * On entry, the symmetric matrix A. If UPLO = 'U', the leading 00040 * N-by-N upper triangular part of A contains the upper 00041 * triangular part of the matrix A, and the strictly lower 00042 * triangular part of A is not referenced. If UPLO = 'L', the 00043 * leading N-by-N lower triangular part of A contains the lower 00044 * triangular part of the matrix A, and the strictly upper 00045 * triangular part of A is not referenced. 00046 * 00047 * On exit, if INFO = 0, the factor U or L from the Cholesky 00048 * factorization A = U**T*U or A = L*L**T. 00049 * 00050 * LDA (input) INTEGER 00051 * The leading dimension of the array A. LDA >= max(1,N). 00052 * 00053 * INFO (output) INTEGER 00054 * = 0: successful exit 00055 * < 0: if INFO = -i, the i-th argument had an illegal value 00056 * > 0: if INFO = i, the leading minor of order i is not 00057 * positive definite, and the factorization could not be 00058 * completed. 00059 * 00060 * ===================================================================== 00061 * 00062 * .. Parameters .. 00063 DOUBLE PRECISION ONE 00064 PARAMETER ( ONE = 1.0D+0 ) 00065 * .. 00066 * .. Local Scalars .. 00067 LOGICAL UPPER 00068 INTEGER J, JB, NB 00069 * .. 00070 * .. External Functions .. 00071 LOGICAL LSAME 00072 INTEGER ILAENV 00073 EXTERNAL LSAME, ILAENV 00074 * .. 00075 * .. External Subroutines .. 00076 EXTERNAL DGEMM, DPOTF2, DSYRK, DTRSM, XERBLA 00077 * .. 00078 * .. Intrinsic Functions .. 00079 INTRINSIC MAX, MIN 00080 * .. 00081 * .. Executable Statements .. 00082 * 00083 * Test the input parameters. 00084 * 00085 INFO = 0 00086 UPPER = LSAME( UPLO, 'U' ) 00087 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00088 INFO = -1 00089 ELSE IF( N.LT.0 ) THEN 00090 INFO = -2 00091 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00092 INFO = -4 00093 END IF 00094 IF( INFO.NE.0 ) THEN 00095 CALL XERBLA( 'DPOTRF', -INFO ) 00096 RETURN 00097 END IF 00098 * 00099 * Quick return if possible 00100 * 00101 IF( N.EQ.0 ) 00102 $ RETURN 00103 * 00104 * Determine the block size for this environment. 00105 * 00106 NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 ) 00107 IF( NB.LE.1 .OR. NB.GE.N ) THEN 00108 * 00109 * Use unblocked code. 00110 * 00111 CALL DPOTF2( UPLO, N, A, LDA, INFO ) 00112 ELSE 00113 * 00114 * Use blocked code. 00115 * 00116 IF( UPPER ) THEN 00117 * 00118 * Compute the Cholesky factorization A = U'*U. 00119 * 00120 DO 10 J = 1, N, NB 00121 * 00122 * Update and factorize the current diagonal block and test 00123 * for non-positive-definiteness. 00124 * 00125 JB = MIN( NB, N-J+1 ) 00126 00127 CALL DPOTF2( 'Upper', JB, A( J, J ), LDA, INFO ) 00128 00129 IF( INFO.NE.0 ) 00130 $ GO TO 30 00131 00132 IF( J+JB.LE.N ) THEN 00133 * 00134 * Updating the trailing submatrix. 00135 * 00136 CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', 00137 $ JB, N-J-JB+1, ONE, A( J, J ), LDA, 00138 $ A( J, J+JB ), LDA ) 00139 CALL DSYRK( 'Upper', 'Transpose', N-J-JB+1, JB, -ONE, 00140 $ A( J, J+JB ), LDA, 00141 $ ONE, A( J+JB, J+JB ), LDA ) 00142 END IF 00143 10 CONTINUE 00144 * 00145 ELSE 00146 * 00147 * Compute the Cholesky factorization A = L*L'. 00148 * 00149 DO 20 J = 1, N, NB 00150 * 00151 * Update and factorize the current diagonal block and test 00152 * for non-positive-definiteness. 00153 * 00154 JB = MIN( NB, N-J+1 ) 00155 00156 CALL DPOTF2( 'Lower', JB, A( J, J ), LDA, INFO ) 00157 00158 IF( INFO.NE.0 ) 00159 $ GO TO 30 00160 00161 IF( J+JB.LE.N ) THEN 00162 * 00163 * Updating the trailing submatrix. 00164 * 00165 CALL DTRSM( 'Right', 'Lower', 'Transpose', 'Non-unit', 00166 $ N-J-JB+1, JB, ONE, A( J, J ), LDA, 00167 $ A( J+JB, J ), LDA ) 00168 00169 CALL DSYRK( 'Lower', 'No Transpose', N-J-JB+1, JB, 00170 $ -ONE, A( J+JB, J ), LDA, 00171 $ ONE, A( J+JB, J+JB ), LDA ) 00172 END IF 00173 20 CONTINUE 00174 END IF 00175 END IF 00176 GO TO 40 00177 * 00178 30 CONTINUE 00179 INFO = INFO + J - 1 00180 * 00181 40 CONTINUE 00182 RETURN 00183 * 00184 * End of DPOTRF 00185 * 00186 END