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