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