1*> \brief \b CHETF2_ROOK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method (unblocked algorithm). 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CHETF2_ROOK + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetf2_rook.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetf2_rook.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetf2_rook.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CHETF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) 22* 23* .. Scalar Arguments .. 24* CHARACTER UPLO 25* INTEGER INFO, LDA, N 26* .. 27* .. Array Arguments .. 28* INTEGER IPIV( * ) 29* COMPLEX A( LDA, * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> CHETF2_ROOK computes the factorization of a complex Hermitian matrix A 39*> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method: 40*> 41*> A = U*D*U**H or A = L*D*L**H 42*> 43*> where U (or L) is a product of permutation and unit upper (lower) 44*> triangular matrices, U**H is the conjugate transpose of U, and D is 45*> Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. 46*> 47*> This is the unblocked version of the algorithm, calling Level 2 BLAS. 48*> \endverbatim 49* 50* Arguments: 51* ========== 52* 53*> \param[in] UPLO 54*> \verbatim 55*> UPLO is CHARACTER*1 56*> Specifies whether the upper or lower triangular part of the 57*> Hermitian matrix A is stored: 58*> = 'U': Upper triangular 59*> = 'L': Lower triangular 60*> \endverbatim 61*> 62*> \param[in] N 63*> \verbatim 64*> N is INTEGER 65*> The order of the matrix A. N >= 0. 66*> \endverbatim 67*> 68*> \param[in,out] A 69*> \verbatim 70*> A is COMPLEX array, dimension (LDA,N) 71*> On entry, the Hermitian matrix A. If UPLO = 'U', the leading 72*> n-by-n upper triangular part of A contains the upper 73*> triangular part of the matrix A, and the strictly lower 74*> triangular part of A is not referenced. If UPLO = 'L', the 75*> leading n-by-n lower triangular part of A contains the lower 76*> triangular part of the matrix A, and the strictly upper 77*> triangular part of A is not referenced. 78*> 79*> On exit, the block diagonal matrix D and the multipliers used 80*> to obtain the factor U or L (see below for further details). 81*> \endverbatim 82*> 83*> \param[in] LDA 84*> \verbatim 85*> LDA is INTEGER 86*> The leading dimension of the array A. LDA >= max(1,N). 87*> \endverbatim 88*> 89*> \param[out] IPIV 90*> \verbatim 91*> IPIV is INTEGER array, dimension (N) 92*> Details of the interchanges and the block structure of D. 93*> 94*> If UPLO = 'U': 95*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were 96*> interchanged and D(k,k) is a 1-by-1 diagonal block. 97*> 98*> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and 99*> columns k and -IPIV(k) were interchanged and rows and 100*> columns k-1 and -IPIV(k-1) were inerchaged, 101*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. 102*> 103*> If UPLO = 'L': 104*> If IPIV(k) > 0, then rows and columns k and IPIV(k) 105*> were interchanged and D(k,k) is a 1-by-1 diagonal block. 106*> 107*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and 108*> columns k and -IPIV(k) were interchanged and rows and 109*> columns k+1 and -IPIV(k+1) were inerchaged, 110*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 111*> \endverbatim 112*> 113*> \param[out] INFO 114*> \verbatim 115*> INFO is INTEGER 116*> = 0: successful exit 117*> < 0: if INFO = -k, the k-th argument had an illegal value 118*> > 0: if INFO = k, D(k,k) is exactly zero. The factorization 119*> has been completed, but the block diagonal matrix D is 120*> exactly singular, and division by zero will occur if it 121*> is used to solve a system of equations. 122*> \endverbatim 123* 124* Authors: 125* ======== 126* 127*> \author Univ. of Tennessee 128*> \author Univ. of California Berkeley 129*> \author Univ. of Colorado Denver 130*> \author NAG Ltd. 131* 132*> \ingroup complexHEcomputational 133* 134*> \par Further Details: 135* ===================== 136*> 137*> \verbatim 138*> 139*> If UPLO = 'U', then A = U*D*U**H, where 140*> U = P(n)*U(n)* ... *P(k)U(k)* ..., 141*> i.e., U is a product of terms P(k)*U(k), where k decreases from n to 142*> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 143*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 144*> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such 145*> that if the diagonal block D(k) is of order s (s = 1 or 2), then 146*> 147*> ( I v 0 ) k-s 148*> U(k) = ( 0 I 0 ) s 149*> ( 0 0 I ) n-k 150*> k-s s n-k 151*> 152*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). 153*> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), 154*> and A(k,k), and v overwrites A(1:k-2,k-1:k). 155*> 156*> If UPLO = 'L', then A = L*D*L**H, where 157*> L = P(1)*L(1)* ... *P(k)*L(k)* ..., 158*> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to 159*> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 160*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 161*> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such 162*> that if the diagonal block D(k) is of order s (s = 1 or 2), then 163*> 164*> ( I 0 0 ) k-1 165*> L(k) = ( 0 I 0 ) s 166*> ( 0 v I ) n-k-s+1 167*> k-1 s n-k-s+1 168*> 169*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). 170*> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), 171*> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). 172*> \endverbatim 173* 174*> \par Contributors: 175* ================== 176*> 177*> \verbatim 178*> 179*> November 2013, Igor Kozachenko, 180*> Computer Science Division, 181*> University of California, Berkeley 182*> 183*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, 184*> School of Mathematics, 185*> University of Manchester 186*> 187*> 01-01-96 - Based on modifications by 188*> J. Lewis, Boeing Computer Services Company 189*> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA 190*> \endverbatim 191* 192* ===================================================================== 193 SUBROUTINE CHETF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) 194* 195* -- LAPACK computational routine -- 196* -- LAPACK is a software package provided by Univ. of Tennessee, -- 197* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 198* 199* .. Scalar Arguments .. 200 CHARACTER UPLO 201 INTEGER INFO, LDA, N 202* .. 203* .. Array Arguments .. 204 INTEGER IPIV( * ) 205 COMPLEX A( LDA, * ) 206* .. 207* 208* ====================================================================== 209* 210* .. Parameters .. 211 REAL ZERO, ONE 212 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 213 REAL EIGHT, SEVTEN 214 PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 ) 215* .. 216* .. Local Scalars .. 217 LOGICAL DONE, UPPER 218 INTEGER I, II, IMAX, ITEMP, J, JMAX, K, KK, KP, KSTEP, 219 $ P 220 REAL ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, STEMP, 221 $ ROWMAX, TT, SFMIN 222 COMPLEX D12, D21, T, WK, WKM1, WKP1, Z 223* .. 224* .. External Functions .. 225* 226 LOGICAL LSAME 227 INTEGER ICAMAX 228 REAL SLAMCH, SLAPY2 229 EXTERNAL LSAME, ICAMAX, SLAMCH, SLAPY2 230* .. 231* .. External Subroutines .. 232 EXTERNAL XERBLA, CSSCAL, CHER, CSWAP 233* .. 234* .. Intrinsic Functions .. 235 INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, REAL, SQRT 236* .. 237* .. Statement Functions .. 238 REAL CABS1 239* .. 240* .. Statement Function definitions .. 241 CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) ) 242* .. 243* .. Executable Statements .. 244* 245* Test the input parameters. 246* 247 INFO = 0 248 UPPER = LSAME( UPLO, 'U' ) 249 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 250 INFO = -1 251 ELSE IF( N.LT.0 ) THEN 252 INFO = -2 253 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 254 INFO = -4 255 END IF 256 IF( INFO.NE.0 ) THEN 257 CALL XERBLA( 'CHETF2_ROOK', -INFO ) 258 RETURN 259 END IF 260* 261* Initialize ALPHA for use in choosing pivot block size. 262* 263 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 264* 265* Compute machine safe minimum 266* 267 SFMIN = SLAMCH( 'S' ) 268* 269 IF( UPPER ) THEN 270* 271* Factorize A as U*D*U**H using the upper triangle of A 272* 273* K is the main loop index, decreasing from N to 1 in steps of 274* 1 or 2 275* 276 K = N 277 10 CONTINUE 278* 279* If K < 1, exit from loop 280* 281 IF( K.LT.1 ) 282 $ GO TO 70 283 KSTEP = 1 284 P = K 285* 286* Determine rows and columns to be interchanged and whether 287* a 1-by-1 or 2-by-2 pivot block will be used 288* 289 ABSAKK = ABS( REAL( A( K, K ) ) ) 290* 291* IMAX is the row-index of the largest off-diagonal element in 292* column K, and COLMAX is its absolute value. 293* Determine both COLMAX and IMAX. 294* 295 IF( K.GT.1 ) THEN 296 IMAX = ICAMAX( K-1, A( 1, K ), 1 ) 297 COLMAX = CABS1( A( IMAX, K ) ) 298 ELSE 299 COLMAX = ZERO 300 END IF 301* 302 IF( ( MAX( ABSAKK, COLMAX ).EQ.ZERO ) ) THEN 303* 304* Column K is zero or underflow: set INFO and continue 305* 306 IF( INFO.EQ.0 ) 307 $ INFO = K 308 KP = K 309 A( K, K ) = REAL( A( K, K ) ) 310 ELSE 311* 312* ============================================================ 313* 314* BEGIN pivot search 315* 316* Case(1) 317* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX 318* (used to handle NaN and Inf) 319* 320 IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN 321* 322* no interchange, use 1-by-1 pivot block 323* 324 KP = K 325* 326 ELSE 327* 328 DONE = .FALSE. 329* 330* Loop until pivot found 331* 332 12 CONTINUE 333* 334* BEGIN pivot search loop body 335* 336* 337* JMAX is the column-index of the largest off-diagonal 338* element in row IMAX, and ROWMAX is its absolute value. 339* Determine both ROWMAX and JMAX. 340* 341 IF( IMAX.NE.K ) THEN 342 JMAX = IMAX + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ), 343 $ LDA ) 344 ROWMAX = CABS1( A( IMAX, JMAX ) ) 345 ELSE 346 ROWMAX = ZERO 347 END IF 348* 349 IF( IMAX.GT.1 ) THEN 350 ITEMP = ICAMAX( IMAX-1, A( 1, IMAX ), 1 ) 351 STEMP = CABS1( A( ITEMP, IMAX ) ) 352 IF( STEMP.GT.ROWMAX ) THEN 353 ROWMAX = STEMP 354 JMAX = ITEMP 355 END IF 356 END IF 357* 358* Case(2) 359* Equivalent to testing for 360* ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX 361* (used to handle NaN and Inf) 362* 363 IF( .NOT.( ABS( REAL( A( IMAX, IMAX ) ) ) 364 $ .LT.ALPHA*ROWMAX ) ) THEN 365* 366* interchange rows and columns K and IMAX, 367* use 1-by-1 pivot block 368* 369 KP = IMAX 370 DONE = .TRUE. 371* 372* Case(3) 373* Equivalent to testing for ROWMAX.EQ.COLMAX, 374* (used to handle NaN and Inf) 375* 376 ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) ) 377 $ THEN 378* 379* interchange rows and columns K-1 and IMAX, 380* use 2-by-2 pivot block 381* 382 KP = IMAX 383 KSTEP = 2 384 DONE = .TRUE. 385* 386* Case(4) 387 ELSE 388* 389* Pivot not found: set params and repeat 390* 391 P = IMAX 392 COLMAX = ROWMAX 393 IMAX = JMAX 394 END IF 395* 396* END pivot search loop body 397* 398 IF( .NOT.DONE ) GOTO 12 399* 400 END IF 401* 402* END pivot search 403* 404* ============================================================ 405* 406* KK is the column of A where pivoting step stopped 407* 408 KK = K - KSTEP + 1 409* 410* For only a 2x2 pivot, interchange rows and columns K and P 411* in the leading submatrix A(1:k,1:k) 412* 413 IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN 414* (1) Swap columnar parts 415 IF( P.GT.1 ) 416 $ CALL CSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 ) 417* (2) Swap and conjugate middle parts 418 DO 14 J = P + 1, K - 1 419 T = CONJG( A( J, K ) ) 420 A( J, K ) = CONJG( A( P, J ) ) 421 A( P, J ) = T 422 14 CONTINUE 423* (3) Swap and conjugate corner elements at row-col interserction 424 A( P, K ) = CONJG( A( P, K ) ) 425* (4) Swap diagonal elements at row-col intersection 426 R1 = REAL( A( K, K ) ) 427 A( K, K ) = REAL( A( P, P ) ) 428 A( P, P ) = R1 429 END IF 430* 431* For both 1x1 and 2x2 pivots, interchange rows and 432* columns KK and KP in the leading submatrix A(1:k,1:k) 433* 434 IF( KP.NE.KK ) THEN 435* (1) Swap columnar parts 436 IF( KP.GT.1 ) 437 $ CALL CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 ) 438* (2) Swap and conjugate middle parts 439 DO 15 J = KP + 1, KK - 1 440 T = CONJG( A( J, KK ) ) 441 A( J, KK ) = CONJG( A( KP, J ) ) 442 A( KP, J ) = T 443 15 CONTINUE 444* (3) Swap and conjugate corner elements at row-col interserction 445 A( KP, KK ) = CONJG( A( KP, KK ) ) 446* (4) Swap diagonal elements at row-col intersection 447 R1 = REAL( A( KK, KK ) ) 448 A( KK, KK ) = REAL( A( KP, KP ) ) 449 A( KP, KP ) = R1 450* 451 IF( KSTEP.EQ.2 ) THEN 452* (*) Make sure that diagonal element of pivot is real 453 A( K, K ) = REAL( A( K, K ) ) 454* (5) Swap row elements 455 T = A( K-1, K ) 456 A( K-1, K ) = A( KP, K ) 457 A( KP, K ) = T 458 END IF 459 ELSE 460* (*) Make sure that diagonal element of pivot is real 461 A( K, K ) = REAL( A( K, K ) ) 462 IF( KSTEP.EQ.2 ) 463 $ A( K-1, K-1 ) = REAL( A( K-1, K-1 ) ) 464 END IF 465* 466* Update the leading submatrix 467* 468 IF( KSTEP.EQ.1 ) THEN 469* 470* 1-by-1 pivot block D(k): column k now holds 471* 472* W(k) = U(k)*D(k) 473* 474* where U(k) is the k-th column of U 475* 476 IF( K.GT.1 ) THEN 477* 478* Perform a rank-1 update of A(1:k-1,1:k-1) and 479* store U(k) in column k 480* 481 IF( ABS( REAL( A( K, K ) ) ).GE.SFMIN ) THEN 482* 483* Perform a rank-1 update of A(1:k-1,1:k-1) as 484* A := A - U(k)*D(k)*U(k)**T 485* = A - W(k)*1/D(k)*W(k)**T 486* 487 D11 = ONE / REAL( A( K, K ) ) 488 CALL CHER( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA ) 489* 490* Store U(k) in column k 491* 492 CALL CSSCAL( K-1, D11, A( 1, K ), 1 ) 493 ELSE 494* 495* Store L(k) in column K 496* 497 D11 = REAL( A( K, K ) ) 498 DO 16 II = 1, K - 1 499 A( II, K ) = A( II, K ) / D11 500 16 CONTINUE 501* 502* Perform a rank-1 update of A(k+1:n,k+1:n) as 503* A := A - U(k)*D(k)*U(k)**T 504* = A - W(k)*(1/D(k))*W(k)**T 505* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T 506* 507 CALL CHER( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA ) 508 END IF 509 END IF 510* 511 ELSE 512* 513* 2-by-2 pivot block D(k): columns k and k-1 now hold 514* 515* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 516* 517* where U(k) and U(k-1) are the k-th and (k-1)-th columns 518* of U 519* 520* Perform a rank-2 update of A(1:k-2,1:k-2) as 521* 522* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T 523* = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T 524* 525* and store L(k) and L(k+1) in columns k and k+1 526* 527 IF( K.GT.2 ) THEN 528* D = |A12| 529 D = SLAPY2( REAL( A( K-1, K ) ), 530 $ AIMAG( A( K-1, K ) ) ) 531 D11 = REAL( A( K, K ) / D ) 532 D22 = REAL( A( K-1, K-1 ) / D ) 533 D12 = A( K-1, K ) / D 534 TT = ONE / ( D11*D22-ONE ) 535* 536 DO 30 J = K - 2, 1, -1 537* 538* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J 539* 540 WKM1 = TT*( D11*A( J, K-1 )-CONJG( D12 )* 541 $ A( J, K ) ) 542 WK = TT*( D22*A( J, K )-D12*A( J, K-1 ) ) 543* 544* Perform a rank-2 update of A(1:k-2,1:k-2) 545* 546 DO 20 I = J, 1, -1 547 A( I, J ) = A( I, J ) - 548 $ ( A( I, K ) / D )*CONJG( WK ) - 549 $ ( A( I, K-1 ) / D )*CONJG( WKM1 ) 550 20 CONTINUE 551* 552* Store U(k) and U(k-1) in cols k and k-1 for row J 553* 554 A( J, K ) = WK / D 555 A( J, K-1 ) = WKM1 / D 556* (*) Make sure that diagonal element of pivot is real 557 A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO ) 558* 559 30 CONTINUE 560* 561 END IF 562* 563 END IF 564* 565 END IF 566* 567* Store details of the interchanges in IPIV 568* 569 IF( KSTEP.EQ.1 ) THEN 570 IPIV( K ) = KP 571 ELSE 572 IPIV( K ) = -P 573 IPIV( K-1 ) = -KP 574 END IF 575* 576* Decrease K and return to the start of the main loop 577* 578 K = K - KSTEP 579 GO TO 10 580* 581 ELSE 582* 583* Factorize A as L*D*L**H using the lower triangle of A 584* 585* K is the main loop index, increasing from 1 to N in steps of 586* 1 or 2 587* 588 K = 1 589 40 CONTINUE 590* 591* If K > N, exit from loop 592* 593 IF( K.GT.N ) 594 $ GO TO 70 595 KSTEP = 1 596 P = K 597* 598* Determine rows and columns to be interchanged and whether 599* a 1-by-1 or 2-by-2 pivot block will be used 600* 601 ABSAKK = ABS( REAL( A( K, K ) ) ) 602* 603* IMAX is the row-index of the largest off-diagonal element in 604* column K, and COLMAX is its absolute value. 605* Determine both COLMAX and IMAX. 606* 607 IF( K.LT.N ) THEN 608 IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 ) 609 COLMAX = CABS1( A( IMAX, K ) ) 610 ELSE 611 COLMAX = ZERO 612 END IF 613* 614 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 615* 616* Column K is zero or underflow: set INFO and continue 617* 618 IF( INFO.EQ.0 ) 619 $ INFO = K 620 KP = K 621 A( K, K ) = REAL( A( K, K ) ) 622 ELSE 623* 624* ============================================================ 625* 626* BEGIN pivot search 627* 628* Case(1) 629* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX 630* (used to handle NaN and Inf) 631* 632 IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN 633* 634* no interchange, use 1-by-1 pivot block 635* 636 KP = K 637* 638 ELSE 639* 640 DONE = .FALSE. 641* 642* Loop until pivot found 643* 644 42 CONTINUE 645* 646* BEGIN pivot search loop body 647* 648* 649* JMAX is the column-index of the largest off-diagonal 650* element in row IMAX, and ROWMAX is its absolute value. 651* Determine both ROWMAX and JMAX. 652* 653 IF( IMAX.NE.K ) THEN 654 JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA ) 655 ROWMAX = CABS1( A( IMAX, JMAX ) ) 656 ELSE 657 ROWMAX = ZERO 658 END IF 659* 660 IF( IMAX.LT.N ) THEN 661 ITEMP = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ), 662 $ 1 ) 663 STEMP = CABS1( A( ITEMP, IMAX ) ) 664 IF( STEMP.GT.ROWMAX ) THEN 665 ROWMAX = STEMP 666 JMAX = ITEMP 667 END IF 668 END IF 669* 670* Case(2) 671* Equivalent to testing for 672* ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX 673* (used to handle NaN and Inf) 674* 675 IF( .NOT.( ABS( REAL( A( IMAX, IMAX ) ) ) 676 $ .LT.ALPHA*ROWMAX ) ) THEN 677* 678* interchange rows and columns K and IMAX, 679* use 1-by-1 pivot block 680* 681 KP = IMAX 682 DONE = .TRUE. 683* 684* Case(3) 685* Equivalent to testing for ROWMAX.EQ.COLMAX, 686* (used to handle NaN and Inf) 687* 688 ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) ) 689 $ THEN 690* 691* interchange rows and columns K+1 and IMAX, 692* use 2-by-2 pivot block 693* 694 KP = IMAX 695 KSTEP = 2 696 DONE = .TRUE. 697* 698* Case(4) 699 ELSE 700* 701* Pivot not found: set params and repeat 702* 703 P = IMAX 704 COLMAX = ROWMAX 705 IMAX = JMAX 706 END IF 707* 708* 709* END pivot search loop body 710* 711 IF( .NOT.DONE ) GOTO 42 712* 713 END IF 714* 715* END pivot search 716* 717* ============================================================ 718* 719* KK is the column of A where pivoting step stopped 720* 721 KK = K + KSTEP - 1 722* 723* For only a 2x2 pivot, interchange rows and columns K and P 724* in the trailing submatrix A(k:n,k:n) 725* 726 IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN 727* (1) Swap columnar parts 728 IF( P.LT.N ) 729 $ CALL CSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 ) 730* (2) Swap and conjugate middle parts 731 DO 44 J = K + 1, P - 1 732 T = CONJG( A( J, K ) ) 733 A( J, K ) = CONJG( A( P, J ) ) 734 A( P, J ) = T 735 44 CONTINUE 736* (3) Swap and conjugate corner elements at row-col interserction 737 A( P, K ) = CONJG( A( P, K ) ) 738* (4) Swap diagonal elements at row-col intersection 739 R1 = REAL( A( K, K ) ) 740 A( K, K ) = REAL( A( P, P ) ) 741 A( P, P ) = R1 742 END IF 743* 744* For both 1x1 and 2x2 pivots, interchange rows and 745* columns KK and KP in the trailing submatrix A(k:n,k:n) 746* 747 IF( KP.NE.KK ) THEN 748* (1) Swap columnar parts 749 IF( KP.LT.N ) 750 $ CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 ) 751* (2) Swap and conjugate middle parts 752 DO 45 J = KK + 1, KP - 1 753 T = CONJG( A( J, KK ) ) 754 A( J, KK ) = CONJG( A( KP, J ) ) 755 A( KP, J ) = T 756 45 CONTINUE 757* (3) Swap and conjugate corner elements at row-col interserction 758 A( KP, KK ) = CONJG( A( KP, KK ) ) 759* (4) Swap diagonal elements at row-col intersection 760 R1 = REAL( A( KK, KK ) ) 761 A( KK, KK ) = REAL( A( KP, KP ) ) 762 A( KP, KP ) = R1 763* 764 IF( KSTEP.EQ.2 ) THEN 765* (*) Make sure that diagonal element of pivot is real 766 A( K, K ) = REAL( A( K, K ) ) 767* (5) Swap row elements 768 T = A( K+1, K ) 769 A( K+1, K ) = A( KP, K ) 770 A( KP, K ) = T 771 END IF 772 ELSE 773* (*) Make sure that diagonal element of pivot is real 774 A( K, K ) = REAL( A( K, K ) ) 775 IF( KSTEP.EQ.2 ) 776 $ A( K+1, K+1 ) = REAL( A( K+1, K+1 ) ) 777 END IF 778* 779* Update the trailing submatrix 780* 781 IF( KSTEP.EQ.1 ) THEN 782* 783* 1-by-1 pivot block D(k): column k of A now holds 784* 785* W(k) = L(k)*D(k), 786* 787* where L(k) is the k-th column of L 788* 789 IF( K.LT.N ) THEN 790* 791* Perform a rank-1 update of A(k+1:n,k+1:n) and 792* store L(k) in column k 793* 794* Handle division by a small number 795* 796 IF( ABS( REAL( A( K, K ) ) ).GE.SFMIN ) THEN 797* 798* Perform a rank-1 update of A(k+1:n,k+1:n) as 799* A := A - L(k)*D(k)*L(k)**T 800* = A - W(k)*(1/D(k))*W(k)**T 801* 802 D11 = ONE / REAL( A( K, K ) ) 803 CALL CHER( UPLO, N-K, -D11, A( K+1, K ), 1, 804 $ A( K+1, K+1 ), LDA ) 805* 806* Store L(k) in column k 807* 808 CALL CSSCAL( N-K, D11, A( K+1, K ), 1 ) 809 ELSE 810* 811* Store L(k) in column k 812* 813 D11 = REAL( A( K, K ) ) 814 DO 46 II = K + 1, N 815 A( II, K ) = A( II, K ) / D11 816 46 CONTINUE 817* 818* Perform a rank-1 update of A(k+1:n,k+1:n) as 819* A := A - L(k)*D(k)*L(k)**T 820* = A - W(k)*(1/D(k))*W(k)**T 821* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T 822* 823 CALL CHER( UPLO, N-K, -D11, A( K+1, K ), 1, 824 $ A( K+1, K+1 ), LDA ) 825 END IF 826 END IF 827* 828 ELSE 829* 830* 2-by-2 pivot block D(k): columns k and k+1 now hold 831* 832* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) 833* 834* where L(k) and L(k+1) are the k-th and (k+1)-th columns 835* of L 836* 837* 838* Perform a rank-2 update of A(k+2:n,k+2:n) as 839* 840* A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T 841* = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T 842* 843* and store L(k) and L(k+1) in columns k and k+1 844* 845 IF( K.LT.N-1 ) THEN 846* D = |A21| 847 D = SLAPY2( REAL( A( K+1, K ) ), 848 $ AIMAG( A( K+1, K ) ) ) 849 D11 = REAL( A( K+1, K+1 ) ) / D 850 D22 = REAL( A( K, K ) ) / D 851 D21 = A( K+1, K ) / D 852 TT = ONE / ( D11*D22-ONE ) 853* 854 DO 60 J = K + 2, N 855* 856* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J 857* 858 WK = TT*( D11*A( J, K )-D21*A( J, K+1 ) ) 859 WKP1 = TT*( D22*A( J, K+1 )-CONJG( D21 )* 860 $ A( J, K ) ) 861* 862* Perform a rank-2 update of A(k+2:n,k+2:n) 863* 864 DO 50 I = J, N 865 A( I, J ) = A( I, J ) - 866 $ ( A( I, K ) / D )*CONJG( WK ) - 867 $ ( A( I, K+1 ) / D )*CONJG( WKP1 ) 868 50 CONTINUE 869* 870* Store L(k) and L(k+1) in cols k and k+1 for row J 871* 872 A( J, K ) = WK / D 873 A( J, K+1 ) = WKP1 / D 874* (*) Make sure that diagonal element of pivot is real 875 A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO ) 876* 877 60 CONTINUE 878* 879 END IF 880* 881 END IF 882* 883 END IF 884* 885* Store details of the interchanges in IPIV 886* 887 IF( KSTEP.EQ.1 ) THEN 888 IPIV( K ) = KP 889 ELSE 890 IPIV( K ) = -P 891 IPIV( K+1 ) = -KP 892 END IF 893* 894* Increase K and return to the start of the main loop 895* 896 K = K + KSTEP 897 GO TO 40 898* 899 END IF 900* 901 70 CONTINUE 902* 903 RETURN 904* 905* End of CHETF2_ROOK 906* 907 END 908