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