![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DPFTRF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DPFTRF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpftrf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpftrf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpftrf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DPFTRF( TRANSR, UPLO, N, A, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER TRANSR, UPLO 00025 * INTEGER N, INFO 00026 * .. 00027 * .. Array Arguments .. 00028 * DOUBLE PRECISION A( 0: * ) 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> DPFTRF computes the Cholesky factorization of a real symmetric 00037 *> positive definite matrix A. 00038 *> 00039 *> The factorization has the form 00040 *> A = U**T * U, if UPLO = 'U', or 00041 *> A = L * L**T, if UPLO = 'L', 00042 *> where U is an upper triangular matrix and L is lower triangular. 00043 *> 00044 *> This is the block version of the algorithm, calling Level 3 BLAS. 00045 *> \endverbatim 00046 * 00047 * Arguments: 00048 * ========== 00049 * 00050 *> \param[in] TRANSR 00051 *> \verbatim 00052 *> TRANSR is CHARACTER*1 00053 *> = 'N': The Normal TRANSR of RFP A is stored; 00054 *> = 'T': The Transpose TRANSR of RFP A is stored. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] UPLO 00058 *> \verbatim 00059 *> UPLO is CHARACTER*1 00060 *> = 'U': Upper triangle of RFP A is stored; 00061 *> = 'L': Lower triangle of RFP A is stored. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] N 00065 *> \verbatim 00066 *> N is INTEGER 00067 *> The order of the matrix A. N >= 0. 00068 *> \endverbatim 00069 *> 00070 *> \param[in,out] A 00071 *> \verbatim 00072 *> A is DOUBLE PRECISION array, dimension ( N*(N+1)/2 ); 00073 *> On entry, the symmetric matrix A in RFP format. RFP format is 00074 *> described by TRANSR, UPLO, and N as follows: If TRANSR = 'N' 00075 *> then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is 00076 *> (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is 00077 *> the transpose of RFP A as defined when 00078 *> TRANSR = 'N'. The contents of RFP A are defined by UPLO as 00079 *> follows: If UPLO = 'U' the RFP A contains the NT elements of 00080 *> upper packed A. If UPLO = 'L' the RFP A contains the elements 00081 *> of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR = 00082 *> 'T'. When TRANSR is 'N' the LDA is N+1 when N is even and N 00083 *> is odd. See the Note below for more details. 00084 *> 00085 *> On exit, if INFO = 0, the factor U or L from the Cholesky 00086 *> factorization RFP A = U**T*U or RFP A = L*L**T. 00087 *> \endverbatim 00088 *> 00089 *> \param[out] INFO 00090 *> \verbatim 00091 *> INFO is INTEGER 00092 *> = 0: successful exit 00093 *> < 0: if INFO = -i, the i-th argument had an illegal value 00094 *> > 0: if INFO = i, the leading minor of order i is not 00095 *> positive definite, and the factorization could not be 00096 *> completed. 00097 *> \endverbatim 00098 * 00099 * Authors: 00100 * ======== 00101 * 00102 *> \author Univ. of Tennessee 00103 *> \author Univ. of California Berkeley 00104 *> \author Univ. of Colorado Denver 00105 *> \author NAG Ltd. 00106 * 00107 *> \date November 2011 00108 * 00109 *> \ingroup doubleOTHERcomputational 00110 * 00111 *> \par Further Details: 00112 * ===================== 00113 *> 00114 *> \verbatim 00115 *> 00116 *> We first consider Rectangular Full Packed (RFP) Format when N is 00117 *> even. We give an example where N = 6. 00118 *> 00119 *> AP is Upper AP is Lower 00120 *> 00121 *> 00 01 02 03 04 05 00 00122 *> 11 12 13 14 15 10 11 00123 *> 22 23 24 25 20 21 22 00124 *> 33 34 35 30 31 32 33 00125 *> 44 45 40 41 42 43 44 00126 *> 55 50 51 52 53 54 55 00127 *> 00128 *> 00129 *> Let TRANSR = 'N'. RFP holds AP as follows: 00130 *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 00131 *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of 00132 *> the transpose of the first three columns of AP upper. 00133 *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 00134 *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of 00135 *> the transpose of the last three columns of AP lower. 00136 *> This covers the case N even and TRANSR = 'N'. 00137 *> 00138 *> RFP A RFP A 00139 *> 00140 *> 03 04 05 33 43 53 00141 *> 13 14 15 00 44 54 00142 *> 23 24 25 10 11 55 00143 *> 33 34 35 20 21 22 00144 *> 00 44 45 30 31 32 00145 *> 01 11 55 40 41 42 00146 *> 02 12 22 50 51 52 00147 *> 00148 *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 00149 *> transpose of RFP A above. One therefore gets: 00150 *> 00151 *> 00152 *> RFP A RFP A 00153 *> 00154 *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50 00155 *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51 00156 *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52 00157 *> 00158 *> 00159 *> We then consider Rectangular Full Packed (RFP) Format when N is 00160 *> odd. We give an example where N = 5. 00161 *> 00162 *> AP is Upper AP is Lower 00163 *> 00164 *> 00 01 02 03 04 00 00165 *> 11 12 13 14 10 11 00166 *> 22 23 24 20 21 22 00167 *> 33 34 30 31 32 33 00168 *> 44 40 41 42 43 44 00169 *> 00170 *> 00171 *> Let TRANSR = 'N'. RFP holds AP as follows: 00172 *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 00173 *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of 00174 *> the transpose of the first two columns of AP upper. 00175 *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 00176 *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of 00177 *> the transpose of the last two columns of AP lower. 00178 *> This covers the case N odd and TRANSR = 'N'. 00179 *> 00180 *> RFP A RFP A 00181 *> 00182 *> 02 03 04 00 33 43 00183 *> 12 13 14 10 11 44 00184 *> 22 23 24 20 21 22 00185 *> 00 33 34 30 31 32 00186 *> 01 11 44 40 41 42 00187 *> 00188 *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 00189 *> transpose of RFP A above. One therefore gets: 00190 *> 00191 *> RFP A RFP A 00192 *> 00193 *> 02 12 22 00 01 00 10 20 30 40 50 00194 *> 03 13 23 33 11 33 11 21 31 41 51 00195 *> 04 14 24 34 44 43 44 22 32 42 52 00196 *> \endverbatim 00197 *> 00198 * ===================================================================== 00199 SUBROUTINE DPFTRF( TRANSR, UPLO, N, A, INFO ) 00200 * 00201 * -- LAPACK computational routine (version 3.4.0) -- 00202 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00203 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00204 * November 2011 00205 * 00206 * .. Scalar Arguments .. 00207 CHARACTER TRANSR, UPLO 00208 INTEGER N, INFO 00209 * .. 00210 * .. Array Arguments .. 00211 DOUBLE PRECISION A( 0: * ) 00212 * 00213 * ===================================================================== 00214 * 00215 * .. Parameters .. 00216 DOUBLE PRECISION ONE 00217 PARAMETER ( ONE = 1.0D+0 ) 00218 * .. 00219 * .. Local Scalars .. 00220 LOGICAL LOWER, NISODD, NORMALTRANSR 00221 INTEGER N1, N2, K 00222 * .. 00223 * .. External Functions .. 00224 LOGICAL LSAME 00225 EXTERNAL LSAME 00226 * .. 00227 * .. External Subroutines .. 00228 EXTERNAL XERBLA, DSYRK, DPOTRF, DTRSM 00229 * .. 00230 * .. Intrinsic Functions .. 00231 INTRINSIC MOD 00232 * .. 00233 * .. Executable Statements .. 00234 * 00235 * Test the input parameters. 00236 * 00237 INFO = 0 00238 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00239 LOWER = LSAME( UPLO, 'L' ) 00240 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN 00241 INFO = -1 00242 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00243 INFO = -2 00244 ELSE IF( N.LT.0 ) THEN 00245 INFO = -3 00246 END IF 00247 IF( INFO.NE.0 ) THEN 00248 CALL XERBLA( 'DPFTRF', -INFO ) 00249 RETURN 00250 END IF 00251 * 00252 * Quick return if possible 00253 * 00254 IF( N.EQ.0 ) 00255 $ RETURN 00256 * 00257 * If N is odd, set NISODD = .TRUE. 00258 * If N is even, set K = N/2 and NISODD = .FALSE. 00259 * 00260 IF( MOD( N, 2 ).EQ.0 ) THEN 00261 K = N / 2 00262 NISODD = .FALSE. 00263 ELSE 00264 NISODD = .TRUE. 00265 END IF 00266 * 00267 * Set N1 and N2 depending on LOWER 00268 * 00269 IF( LOWER ) THEN 00270 N2 = N / 2 00271 N1 = N - N2 00272 ELSE 00273 N1 = N / 2 00274 N2 = N - N1 00275 END IF 00276 * 00277 * start execution: there are eight cases 00278 * 00279 IF( NISODD ) THEN 00280 * 00281 * N is odd 00282 * 00283 IF( NORMALTRANSR ) THEN 00284 * 00285 * N is odd and TRANSR = 'N' 00286 * 00287 IF( LOWER ) THEN 00288 * 00289 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) ) 00290 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0) 00291 * T1 -> a(0), T2 -> a(n), S -> a(n1) 00292 * 00293 CALL DPOTRF( 'L', N1, A( 0 ), N, INFO ) 00294 IF( INFO.GT.0 ) 00295 $ RETURN 00296 CALL DTRSM( 'R', 'L', 'T', 'N', N2, N1, ONE, A( 0 ), N, 00297 $ A( N1 ), N ) 00298 CALL DSYRK( 'U', 'N', N2, N1, -ONE, A( N1 ), N, ONE, 00299 $ A( N ), N ) 00300 CALL DPOTRF( 'U', N2, A( N ), N, INFO ) 00301 IF( INFO.GT.0 ) 00302 $ INFO = INFO + N1 00303 * 00304 ELSE 00305 * 00306 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1) 00307 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0) 00308 * T1 -> a(n2), T2 -> a(n1), S -> a(0) 00309 * 00310 CALL DPOTRF( 'L', N1, A( N2 ), N, INFO ) 00311 IF( INFO.GT.0 ) 00312 $ RETURN 00313 CALL DTRSM( 'L', 'L', 'N', 'N', N1, N2, ONE, A( N2 ), N, 00314 $ A( 0 ), N ) 00315 CALL DSYRK( 'U', 'T', N2, N1, -ONE, A( 0 ), N, ONE, 00316 $ A( N1 ), N ) 00317 CALL DPOTRF( 'U', N2, A( N1 ), N, INFO ) 00318 IF( INFO.GT.0 ) 00319 $ INFO = INFO + N1 00320 * 00321 END IF 00322 * 00323 ELSE 00324 * 00325 * N is odd and TRANSR = 'T' 00326 * 00327 IF( LOWER ) THEN 00328 * 00329 * SRPA for LOWER, TRANSPOSE and N is odd 00330 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1) 00331 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1 00332 * 00333 CALL DPOTRF( 'U', N1, A( 0 ), N1, INFO ) 00334 IF( INFO.GT.0 ) 00335 $ RETURN 00336 CALL DTRSM( 'L', 'U', 'T', 'N', N1, N2, ONE, A( 0 ), N1, 00337 $ A( N1*N1 ), N1 ) 00338 CALL DSYRK( 'L', 'T', N2, N1, -ONE, A( N1*N1 ), N1, ONE, 00339 $ A( 1 ), N1 ) 00340 CALL DPOTRF( 'L', N2, A( 1 ), N1, INFO ) 00341 IF( INFO.GT.0 ) 00342 $ INFO = INFO + N1 00343 * 00344 ELSE 00345 * 00346 * SRPA for UPPER, TRANSPOSE and N is odd 00347 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0) 00348 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2 00349 * 00350 CALL DPOTRF( 'U', N1, A( N2*N2 ), N2, INFO ) 00351 IF( INFO.GT.0 ) 00352 $ RETURN 00353 CALL DTRSM( 'R', 'U', 'N', 'N', N2, N1, ONE, A( N2*N2 ), 00354 $ N2, A( 0 ), N2 ) 00355 CALL DSYRK( 'L', 'N', N2, N1, -ONE, A( 0 ), N2, ONE, 00356 $ A( N1*N2 ), N2 ) 00357 CALL DPOTRF( 'L', N2, A( N1*N2 ), N2, INFO ) 00358 IF( INFO.GT.0 ) 00359 $ INFO = INFO + N1 00360 * 00361 END IF 00362 * 00363 END IF 00364 * 00365 ELSE 00366 * 00367 * N is even 00368 * 00369 IF( NORMALTRANSR ) THEN 00370 * 00371 * N is even and TRANSR = 'N' 00372 * 00373 IF( LOWER ) THEN 00374 * 00375 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) ) 00376 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0) 00377 * T1 -> a(1), T2 -> a(0), S -> a(k+1) 00378 * 00379 CALL DPOTRF( 'L', K, A( 1 ), N+1, INFO ) 00380 IF( INFO.GT.0 ) 00381 $ RETURN 00382 CALL DTRSM( 'R', 'L', 'T', 'N', K, K, ONE, A( 1 ), N+1, 00383 $ A( K+1 ), N+1 ) 00384 CALL DSYRK( 'U', 'N', K, K, -ONE, A( K+1 ), N+1, ONE, 00385 $ A( 0 ), N+1 ) 00386 CALL DPOTRF( 'U', K, A( 0 ), N+1, INFO ) 00387 IF( INFO.GT.0 ) 00388 $ INFO = INFO + K 00389 * 00390 ELSE 00391 * 00392 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) ) 00393 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0) 00394 * T1 -> a(k+1), T2 -> a(k), S -> a(0) 00395 * 00396 CALL DPOTRF( 'L', K, A( K+1 ), N+1, INFO ) 00397 IF( INFO.GT.0 ) 00398 $ RETURN 00399 CALL DTRSM( 'L', 'L', 'N', 'N', K, K, ONE, A( K+1 ), 00400 $ N+1, A( 0 ), N+1 ) 00401 CALL DSYRK( 'U', 'T', K, K, -ONE, A( 0 ), N+1, ONE, 00402 $ A( K ), N+1 ) 00403 CALL DPOTRF( 'U', K, A( K ), N+1, INFO ) 00404 IF( INFO.GT.0 ) 00405 $ INFO = INFO + K 00406 * 00407 END IF 00408 * 00409 ELSE 00410 * 00411 * N is even and TRANSR = 'T' 00412 * 00413 IF( LOWER ) THEN 00414 * 00415 * SRPA for LOWER, TRANSPOSE and N is even (see paper) 00416 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1) 00417 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k 00418 * 00419 CALL DPOTRF( 'U', K, A( 0+K ), K, INFO ) 00420 IF( INFO.GT.0 ) 00421 $ RETURN 00422 CALL DTRSM( 'L', 'U', 'T', 'N', K, K, ONE, A( K ), N1, 00423 $ A( K*( K+1 ) ), K ) 00424 CALL DSYRK( 'L', 'T', K, K, -ONE, A( K*( K+1 ) ), K, ONE, 00425 $ A( 0 ), K ) 00426 CALL DPOTRF( 'L', K, A( 0 ), K, INFO ) 00427 IF( INFO.GT.0 ) 00428 $ INFO = INFO + K 00429 * 00430 ELSE 00431 * 00432 * SRPA for UPPER, TRANSPOSE and N is even (see paper) 00433 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0) 00434 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k 00435 * 00436 CALL DPOTRF( 'U', K, A( K*( K+1 ) ), K, INFO ) 00437 IF( INFO.GT.0 ) 00438 $ RETURN 00439 CALL DTRSM( 'R', 'U', 'N', 'N', K, K, ONE, 00440 $ A( K*( K+1 ) ), K, A( 0 ), K ) 00441 CALL DSYRK( 'L', 'N', K, K, -ONE, A( 0 ), K, ONE, 00442 $ A( K*K ), K ) 00443 CALL DPOTRF( 'L', K, A( K*K ), K, INFO ) 00444 IF( INFO.GT.0 ) 00445 $ INFO = INFO + K 00446 * 00447 END IF 00448 * 00449 END IF 00450 * 00451 END IF 00452 * 00453 RETURN 00454 * 00455 * End of DPFTRF 00456 * 00457 END