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