1*> \brief \b CPFTRF 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CPFTRF + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpftrf.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpftrf.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpftrf.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CPFTRF( TRANSR, UPLO, N, A, INFO ) 22* 23* .. Scalar Arguments .. 24* CHARACTER TRANSR, UPLO 25* INTEGER N, INFO 26* .. 27* .. Array Arguments .. 28* COMPLEX A( 0: * ) 29* 30* 31*> \par Purpose: 32* ============= 33*> 34*> \verbatim 35*> 36*> CPFTRF computes the Cholesky factorization of a complex Hermitian 37*> positive definite matrix A. 38*> 39*> The factorization has the form 40*> A = U**H * U, if UPLO = 'U', or 41*> A = L * L**H, if UPLO = 'L', 42*> where U is an upper triangular matrix and L is lower triangular. 43*> 44*> This is the block version of the algorithm, calling Level 3 BLAS. 45*> \endverbatim 46* 47* Arguments: 48* ========== 49* 50*> \param[in] TRANSR 51*> \verbatim 52*> TRANSR is CHARACTER*1 53*> = 'N': The Normal TRANSR of RFP A is stored; 54*> = 'C': The Conjugate-transpose TRANSR of RFP A is stored. 55*> \endverbatim 56*> 57*> \param[in] UPLO 58*> \verbatim 59*> UPLO is CHARACTER*1 60*> = 'U': Upper triangle of RFP A is stored; 61*> = 'L': Lower triangle of RFP A is stored. 62*> \endverbatim 63*> 64*> \param[in] N 65*> \verbatim 66*> N is INTEGER 67*> The order of the matrix A. N >= 0. 68*> \endverbatim 69*> 70*> \param[in,out] A 71*> \verbatim 72*> A is COMPLEX array, dimension ( N*(N+1)/2 ); 73*> On entry, the Hermitian matrix A in RFP format. RFP format is 74*> described by TRANSR, UPLO, and N as follows: If TRANSR = 'N' 75*> then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is 76*> (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is 77*> the Conjugate-transpose of RFP A as defined when 78*> TRANSR = 'N'. The contents of RFP A are defined by UPLO as 79*> follows: If UPLO = 'U' the RFP A contains the nt elements of 80*> upper packed A. If UPLO = 'L' the RFP A contains the elements 81*> of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR = 82*> 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N 83*> is odd. See the Note below for more details. 84*> 85*> On exit, if INFO = 0, the factor U or L from the Cholesky 86*> factorization RFP A = U**H*U or RFP A = L*L**H. 87*> \endverbatim 88*> 89*> \param[out] INFO 90*> \verbatim 91*> INFO is INTEGER 92*> = 0: successful exit 93*> < 0: if INFO = -i, the i-th argument had an illegal value 94*> > 0: if INFO = i, the leading minor of order i is not 95*> positive definite, and the factorization could not be 96*> completed. 97*> 98*> Further Notes on RFP Format: 99*> ============================ 100*> 101*> We first consider Standard Packed Format when N is even. 102*> We give an example where N = 6. 103*> 104*> AP is Upper AP is Lower 105*> 106*> 00 01 02 03 04 05 00 107*> 11 12 13 14 15 10 11 108*> 22 23 24 25 20 21 22 109*> 33 34 35 30 31 32 33 110*> 44 45 40 41 42 43 44 111*> 55 50 51 52 53 54 55 112*> 113*> Let TRANSR = 'N'. RFP holds AP as follows: 114*> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 115*> three columns of AP upper. The lower triangle A(4:6,0:2) consists of 116*> conjugate-transpose of the first three columns of AP upper. 117*> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 118*> three columns of AP lower. The upper triangle A(0:2,0:2) consists of 119*> conjugate-transpose of the last three columns of AP lower. 120*> To denote conjugate we place -- above the element. This covers the 121*> case N even and TRANSR = 'N'. 122*> 123*> RFP A RFP A 124*> 125*> -- -- -- 126*> 03 04 05 33 43 53 127*> -- -- 128*> 13 14 15 00 44 54 129*> -- 130*> 23 24 25 10 11 55 131*> 132*> 33 34 35 20 21 22 133*> -- 134*> 00 44 45 30 31 32 135*> -- -- 136*> 01 11 55 40 41 42 137*> -- -- -- 138*> 02 12 22 50 51 52 139*> 140*> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 141*> transpose of RFP A above. One therefore gets: 142*> 143*> RFP A RFP A 144*> 145*> -- -- -- -- -- -- -- -- -- -- 146*> 03 13 23 33 00 01 02 33 00 10 20 30 40 50 147*> -- -- -- -- -- -- -- -- -- -- 148*> 04 14 24 34 44 11 12 43 44 11 21 31 41 51 149*> -- -- -- -- -- -- -- -- -- -- 150*> 05 15 25 35 45 55 22 53 54 55 22 32 42 52 151*> 152*> We next consider Standard Packed Format when N is odd. 153*> We give an example where N = 5. 154*> 155*> AP is Upper AP is Lower 156*> 157*> 00 01 02 03 04 00 158*> 11 12 13 14 10 11 159*> 22 23 24 20 21 22 160*> 33 34 30 31 32 33 161*> 44 40 41 42 43 44 162*> 163*> Let TRANSR = 'N'. RFP holds AP as follows: 164*> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 165*> three columns of AP upper. The lower triangle A(3:4,0:1) consists of 166*> conjugate-transpose of the first two columns of AP upper. 167*> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 168*> three columns of AP lower. The upper triangle A(0:1,1:2) consists of 169*> conjugate-transpose of the last two columns of AP lower. 170*> To denote conjugate we place -- above the element. This covers the 171*> case N odd and TRANSR = 'N'. 172*> 173*> RFP A RFP A 174*> 175*> -- -- 176*> 02 03 04 00 33 43 177*> -- 178*> 12 13 14 10 11 44 179*> 180*> 22 23 24 20 21 22 181*> -- 182*> 00 33 34 30 31 32 183*> -- -- 184*> 01 11 44 40 41 42 185*> 186*> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 187*> transpose of RFP A above. One therefore gets: 188*> 189*> RFP A RFP A 190*> 191*> -- -- -- -- -- -- -- -- -- 192*> 02 12 22 00 01 00 10 20 30 40 50 193*> -- -- -- -- -- -- -- -- -- 194*> 03 13 23 33 11 33 11 21 31 41 51 195*> -- -- -- -- -- -- -- -- -- 196*> 04 14 24 34 44 43 44 22 32 42 52 197*> \endverbatim 198* 199* Authors: 200* ======== 201* 202*> \author Univ. of Tennessee 203*> \author Univ. of California Berkeley 204*> \author Univ. of Colorado Denver 205*> \author NAG Ltd. 206* 207*> \date December 2016 208* 209*> \ingroup complexOTHERcomputational 210* 211* ===================================================================== 212 SUBROUTINE CPFTRF( TRANSR, UPLO, N, A, INFO ) 213* 214* -- LAPACK computational routine (version 3.7.0) -- 215* -- LAPACK is a software package provided by Univ. of Tennessee, -- 216* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 217* December 2016 218* 219* .. Scalar Arguments .. 220 CHARACTER TRANSR, UPLO 221 INTEGER N, INFO 222* .. 223* .. Array Arguments .. 224 COMPLEX A( 0: * ) 225* 226* ===================================================================== 227* 228* .. Parameters .. 229 REAL ONE 230 COMPLEX CONE 231 PARAMETER ( ONE = 1.0E+0, CONE = ( 1.0E+0, 0.0E+0 ) ) 232* .. 233* .. Local Scalars .. 234 LOGICAL LOWER, NISODD, NORMALTRANSR 235 INTEGER N1, N2, K 236* .. 237* .. External Functions .. 238 LOGICAL LSAME 239 EXTERNAL LSAME 240* .. 241* .. External Subroutines .. 242 EXTERNAL XERBLA, CHERK, CPOTRF, CTRSM 243* .. 244* .. Intrinsic Functions .. 245 INTRINSIC MOD 246* .. 247* .. Executable Statements .. 248* 249* Test the input parameters. 250* 251 INFO = 0 252 NORMALTRANSR = LSAME( TRANSR, 'N' ) 253 LOWER = LSAME( UPLO, 'L' ) 254 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 255 INFO = -1 256 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 257 INFO = -2 258 ELSE IF( N.LT.0 ) THEN 259 INFO = -3 260 END IF 261 IF( INFO.NE.0 ) THEN 262 CALL XERBLA( 'CPFTRF', -INFO ) 263 RETURN 264 END IF 265* 266* Quick return if possible 267* 268 IF( N.EQ.0 ) 269 $ RETURN 270* 271* If N is odd, set NISODD = .TRUE. 272* If N is even, set K = N/2 and NISODD = .FALSE. 273* 274 IF( MOD( N, 2 ).EQ.0 ) THEN 275 K = N / 2 276 NISODD = .FALSE. 277 ELSE 278 NISODD = .TRUE. 279 END IF 280* 281* Set N1 and N2 depending on LOWER 282* 283 IF( LOWER ) THEN 284 N2 = N / 2 285 N1 = N - N2 286 ELSE 287 N1 = N / 2 288 N2 = N - N1 289 END IF 290* 291* start execution: there are eight cases 292* 293 IF( NISODD ) THEN 294* 295* N is odd 296* 297 IF( NORMALTRANSR ) THEN 298* 299* N is odd and TRANSR = 'N' 300* 301 IF( LOWER ) THEN 302* 303* SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) ) 304* T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0) 305* T1 -> a(0), T2 -> a(n), S -> a(n1) 306* 307 CALL CPOTRF( 'L', N1, A( 0 ), N, INFO ) 308 IF( INFO.GT.0 ) 309 $ RETURN 310 CALL CTRSM( 'R', 'L', 'C', 'N', N2, N1, CONE, A( 0 ), N, 311 $ A( N1 ), N ) 312 CALL CHERK( 'U', 'N', N2, N1, -ONE, A( N1 ), N, ONE, 313 $ A( N ), N ) 314 CALL CPOTRF( 'U', N2, A( N ), N, INFO ) 315 IF( INFO.GT.0 ) 316 $ INFO = INFO + N1 317* 318 ELSE 319* 320* SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1) 321* T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0) 322* T1 -> a(n2), T2 -> a(n1), S -> a(0) 323* 324 CALL CPOTRF( 'L', N1, A( N2 ), N, INFO ) 325 IF( INFO.GT.0 ) 326 $ RETURN 327 CALL CTRSM( 'L', 'L', 'N', 'N', N1, N2, CONE, A( N2 ), N, 328 $ A( 0 ), N ) 329 CALL CHERK( 'U', 'C', N2, N1, -ONE, A( 0 ), N, ONE, 330 $ A( N1 ), N ) 331 CALL CPOTRF( 'U', N2, A( N1 ), N, INFO ) 332 IF( INFO.GT.0 ) 333 $ INFO = INFO + N1 334* 335 END IF 336* 337 ELSE 338* 339* N is odd and TRANSR = 'C' 340* 341 IF( LOWER ) THEN 342* 343* SRPA for LOWER, TRANSPOSE and N is odd 344* T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1) 345* T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1 346* 347 CALL CPOTRF( 'U', N1, A( 0 ), N1, INFO ) 348 IF( INFO.GT.0 ) 349 $ RETURN 350 CALL CTRSM( 'L', 'U', 'C', 'N', N1, N2, CONE, A( 0 ), N1, 351 $ A( N1*N1 ), N1 ) 352 CALL CHERK( 'L', 'C', N2, N1, -ONE, A( N1*N1 ), N1, ONE, 353 $ A( 1 ), N1 ) 354 CALL CPOTRF( 'L', N2, A( 1 ), N1, INFO ) 355 IF( INFO.GT.0 ) 356 $ INFO = INFO + N1 357* 358 ELSE 359* 360* SRPA for UPPER, TRANSPOSE and N is odd 361* T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0) 362* T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2 363* 364 CALL CPOTRF( 'U', N1, A( N2*N2 ), N2, INFO ) 365 IF( INFO.GT.0 ) 366 $ RETURN 367 CALL CTRSM( 'R', 'U', 'N', 'N', N2, N1, CONE, A( N2*N2 ), 368 $ N2, A( 0 ), N2 ) 369 CALL CHERK( 'L', 'N', N2, N1, -ONE, A( 0 ), N2, ONE, 370 $ A( N1*N2 ), N2 ) 371 CALL CPOTRF( 'L', N2, A( N1*N2 ), N2, INFO ) 372 IF( INFO.GT.0 ) 373 $ INFO = INFO + N1 374* 375 END IF 376* 377 END IF 378* 379 ELSE 380* 381* N is even 382* 383 IF( NORMALTRANSR ) THEN 384* 385* N is even and TRANSR = 'N' 386* 387 IF( LOWER ) THEN 388* 389* SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) ) 390* T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0) 391* T1 -> a(1), T2 -> a(0), S -> a(k+1) 392* 393 CALL CPOTRF( 'L', K, A( 1 ), N+1, INFO ) 394 IF( INFO.GT.0 ) 395 $ RETURN 396 CALL CTRSM( 'R', 'L', 'C', 'N', K, K, CONE, A( 1 ), N+1, 397 $ A( K+1 ), N+1 ) 398 CALL CHERK( 'U', 'N', K, K, -ONE, A( K+1 ), N+1, ONE, 399 $ A( 0 ), N+1 ) 400 CALL CPOTRF( 'U', K, A( 0 ), N+1, INFO ) 401 IF( INFO.GT.0 ) 402 $ INFO = INFO + K 403* 404 ELSE 405* 406* SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) ) 407* T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0) 408* T1 -> a(k+1), T2 -> a(k), S -> a(0) 409* 410 CALL CPOTRF( 'L', K, A( K+1 ), N+1, INFO ) 411 IF( INFO.GT.0 ) 412 $ RETURN 413 CALL CTRSM( 'L', 'L', 'N', 'N', K, K, CONE, A( K+1 ), 414 $ N+1, A( 0 ), N+1 ) 415 CALL CHERK( 'U', 'C', K, K, -ONE, A( 0 ), N+1, ONE, 416 $ A( K ), N+1 ) 417 CALL CPOTRF( 'U', K, A( K ), N+1, INFO ) 418 IF( INFO.GT.0 ) 419 $ INFO = INFO + K 420* 421 END IF 422* 423 ELSE 424* 425* N is even and TRANSR = 'C' 426* 427 IF( LOWER ) THEN 428* 429* SRPA for LOWER, TRANSPOSE and N is even (see paper) 430* T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1) 431* T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k 432* 433 CALL CPOTRF( 'U', K, A( 0+K ), K, INFO ) 434 IF( INFO.GT.0 ) 435 $ RETURN 436 CALL CTRSM( 'L', 'U', 'C', 'N', K, K, CONE, A( K ), N1, 437 $ A( K*( K+1 ) ), K ) 438 CALL CHERK( 'L', 'C', K, K, -ONE, A( K*( K+1 ) ), K, ONE, 439 $ A( 0 ), K ) 440 CALL CPOTRF( 'L', K, A( 0 ), K, INFO ) 441 IF( INFO.GT.0 ) 442 $ INFO = INFO + K 443* 444 ELSE 445* 446* SRPA for UPPER, TRANSPOSE and N is even (see paper) 447* T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0) 448* T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k 449* 450 CALL CPOTRF( 'U', K, A( K*( K+1 ) ), K, INFO ) 451 IF( INFO.GT.0 ) 452 $ RETURN 453 CALL CTRSM( 'R', 'U', 'N', 'N', K, K, CONE, 454 $ A( K*( K+1 ) ), K, A( 0 ), K ) 455 CALL CHERK( 'L', 'N', K, K, -ONE, A( 0 ), K, ONE, 456 $ A( K*K ), K ) 457 CALL CPOTRF( 'L', K, A( K*K ), K, INFO ) 458 IF( INFO.GT.0 ) 459 $ INFO = INFO + K 460* 461 END IF 462* 463 END IF 464* 465 END IF 466* 467 RETURN 468* 469* End of CPFTRF 470* 471 END 472