1*> \brief \b DSYTF2_RK computes the factorization of a real symmetric 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 DSYTF2_RK + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytf2_rk.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytf2_rk.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytf2_rk.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DSYTF2_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* DOUBLE PRECISION A( LDA, * ), E ( * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> DSYTF2_RK computes the factorization of a real symmetric matrix A 38*> using the bounded Bunch-Kaufman (rook) diagonal pivoting method: 39*> 40*> A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T), 41*> 42*> where U (or L) is unit upper (or lower) triangular matrix, 43*> U**T (or L**T) is the transpose of U (or L), P is a permutation 44*> matrix, P**T is the transpose of P, and D is symmetric 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*> symmetric 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 DOUBLE PRECISION array, dimension (LDA,N) 72*> On entry, the symmetric 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 symmetric 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 DOUBLE PRECISION array, dimension (N) 101*> On exit, contains the superdiagonal (or subdiagonal) 102*> elements of the symmetric 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 symmetric 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 doubleSYcomputational 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 DSYTF2_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 DOUBLE PRECISION 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* .. 263* .. Local Scalars .. 264 LOGICAL UPPER, DONE 265 INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP, 266 $ P, II 267 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, 268 $ ROWMAX, DTEMP, T, WK, WKM1, WKP1, SFMIN 269* .. 270* .. External Functions .. 271 LOGICAL LSAME 272 INTEGER IDAMAX 273 DOUBLE PRECISION DLAMCH 274 EXTERNAL LSAME, IDAMAX, DLAMCH 275* .. 276* .. External Subroutines .. 277 EXTERNAL DSCAL, DSWAP, DSYR, XERBLA 278* .. 279* .. Intrinsic Functions .. 280 INTRINSIC ABS, MAX, SQRT 281* .. 282* .. Executable Statements .. 283* 284* Test the input parameters. 285* 286 INFO = 0 287 UPPER = LSAME( UPLO, 'U' ) 288 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 289 INFO = -1 290 ELSE IF( N.LT.0 ) THEN 291 INFO = -2 292 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 293 INFO = -4 294 END IF 295 IF( INFO.NE.0 ) THEN 296 CALL XERBLA( 'DSYTF2_RK', -INFO ) 297 RETURN 298 END IF 299* 300* Initialize ALPHA for use in choosing pivot block size. 301* 302 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 303* 304* Compute machine safe minimum 305* 306 SFMIN = DLAMCH( 'S' ) 307* 308 IF( UPPER ) THEN 309* 310* Factorize A as U*D*U**T using the upper triangle of A 311* 312* Initialize the first entry of array E, where superdiagonal 313* elements of D are stored 314* 315 E( 1 ) = ZERO 316* 317* K is the main loop index, decreasing from N to 1 in steps of 318* 1 or 2 319* 320 K = N 321 10 CONTINUE 322* 323* If K < 1, exit from loop 324* 325 IF( K.LT.1 ) 326 $ GO TO 34 327 KSTEP = 1 328 P = K 329* 330* Determine rows and columns to be interchanged and whether 331* a 1-by-1 or 2-by-2 pivot block will be used 332* 333 ABSAKK = ABS( A( K, K ) ) 334* 335* IMAX is the row-index of the largest off-diagonal element in 336* column K, and COLMAX is its absolute value. 337* Determine both COLMAX and IMAX. 338* 339 IF( K.GT.1 ) THEN 340 IMAX = IDAMAX( K-1, A( 1, K ), 1 ) 341 COLMAX = ABS( A( IMAX, K ) ) 342 ELSE 343 COLMAX = ZERO 344 END IF 345* 346 IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) ) THEN 347* 348* Column K is zero or underflow: set INFO and continue 349* 350 IF( INFO.EQ.0 ) 351 $ INFO = K 352 KP = K 353* 354* Set E( K ) to zero 355* 356 IF( K.GT.1 ) 357 $ E( K ) = ZERO 358* 359 ELSE 360* 361* Test for interchange 362* 363* Equivalent to testing for (used to handle NaN and Inf) 364* ABSAKK.GE.ALPHA*COLMAX 365* 366 IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN 367* 368* no interchange, 369* use 1-by-1 pivot block 370* 371 KP = K 372 ELSE 373* 374 DONE = .FALSE. 375* 376* Loop until pivot found 377* 378 12 CONTINUE 379* 380* Begin pivot search loop body 381* 382* JMAX is the column-index of the largest off-diagonal 383* element in row IMAX, and ROWMAX is its absolute value. 384* Determine both ROWMAX and JMAX. 385* 386 IF( IMAX.NE.K ) THEN 387 JMAX = IMAX + IDAMAX( K-IMAX, A( IMAX, IMAX+1 ), 388 $ LDA ) 389 ROWMAX = ABS( A( IMAX, JMAX ) ) 390 ELSE 391 ROWMAX = ZERO 392 END IF 393* 394 IF( IMAX.GT.1 ) THEN 395 ITEMP = IDAMAX( IMAX-1, A( 1, IMAX ), 1 ) 396 DTEMP = ABS( A( ITEMP, IMAX ) ) 397 IF( DTEMP.GT.ROWMAX ) THEN 398 ROWMAX = DTEMP 399 JMAX = ITEMP 400 END IF 401 END IF 402* 403* Equivalent to testing for (used to handle NaN and Inf) 404* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX 405* 406 IF( .NOT.( ABS( A( IMAX, IMAX ) ).LT.ALPHA*ROWMAX ) ) 407 $ THEN 408* 409* interchange rows and columns K and IMAX, 410* use 1-by-1 pivot block 411* 412 KP = IMAX 413 DONE = .TRUE. 414* 415* Equivalent to testing for ROWMAX .EQ. COLMAX, 416* used to handle NaN and Inf 417* 418 ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN 419* 420* interchange rows and columns K+1 and IMAX, 421* use 2-by-2 pivot block 422* 423 KP = IMAX 424 KSTEP = 2 425 DONE = .TRUE. 426 ELSE 427* 428* Pivot NOT found, set variables and repeat 429* 430 P = IMAX 431 COLMAX = ROWMAX 432 IMAX = JMAX 433 END IF 434* 435* End pivot search loop body 436* 437 IF( .NOT. DONE ) GOTO 12 438* 439 END IF 440* 441* Swap TWO rows and TWO columns 442* 443* First swap 444* 445 IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN 446* 447* Interchange rows and column K and P in the leading 448* submatrix A(1:k,1:k) if we have a 2-by-2 pivot 449* 450 IF( P.GT.1 ) 451 $ CALL DSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 ) 452 IF( P.LT.(K-1) ) 453 $ CALL DSWAP( K-P-1, A( P+1, K ), 1, A( P, P+1 ), 454 $ LDA ) 455 T = A( K, K ) 456 A( K, K ) = A( P, P ) 457 A( P, P ) = T 458* 459* Convert upper triangle of A into U form by applying 460* the interchanges in columns k+1:N. 461* 462 IF( K.LT.N ) 463 $ CALL DSWAP( N-K, A( K, K+1 ), LDA, A( P, K+1 ), LDA ) 464* 465 END IF 466* 467* Second swap 468* 469 KK = K - KSTEP + 1 470 IF( KP.NE.KK ) THEN 471* 472* Interchange rows and columns KK and KP in the leading 473* submatrix A(1:k,1:k) 474* 475 IF( KP.GT.1 ) 476 $ CALL DSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 ) 477 IF( ( KK.GT.1 ) .AND. ( KP.LT.(KK-1) ) ) 478 $ CALL DSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ), 479 $ LDA ) 480 T = A( KK, KK ) 481 A( KK, KK ) = A( KP, KP ) 482 A( KP, KP ) = T 483 IF( KSTEP.EQ.2 ) THEN 484 T = A( K-1, K ) 485 A( K-1, K ) = A( KP, K ) 486 A( KP, K ) = T 487 END IF 488* 489* Convert upper triangle of A into U form by applying 490* the interchanges in columns k+1:N. 491* 492 IF( K.LT.N ) 493 $ CALL DSWAP( N-K, A( KK, K+1 ), LDA, A( KP, K+1 ), 494 $ LDA ) 495* 496 END IF 497* 498* Update the leading submatrix 499* 500 IF( KSTEP.EQ.1 ) THEN 501* 502* 1-by-1 pivot block D(k): column k now holds 503* 504* W(k) = U(k)*D(k) 505* 506* where U(k) is the k-th column of U 507* 508 IF( K.GT.1 ) THEN 509* 510* Perform a rank-1 update of A(1:k-1,1:k-1) and 511* store U(k) in column k 512* 513 IF( ABS( A( K, K ) ).GE.SFMIN ) THEN 514* 515* Perform a rank-1 update of A(1:k-1,1:k-1) as 516* A := A - U(k)*D(k)*U(k)**T 517* = A - W(k)*1/D(k)*W(k)**T 518* 519 D11 = ONE / A( K, K ) 520 CALL DSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA ) 521* 522* Store U(k) in column k 523* 524 CALL DSCAL( K-1, D11, A( 1, K ), 1 ) 525 ELSE 526* 527* Store L(k) in column K 528* 529 D11 = A( K, K ) 530 DO 16 II = 1, K - 1 531 A( II, K ) = A( II, K ) / D11 532 16 CONTINUE 533* 534* Perform a rank-1 update of A(k+1:n,k+1:n) as 535* A := A - U(k)*D(k)*U(k)**T 536* = A - W(k)*(1/D(k))*W(k)**T 537* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T 538* 539 CALL DSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA ) 540 END IF 541* 542* Store the superdiagonal element of D in array E 543* 544 E( K ) = ZERO 545* 546 END IF 547* 548 ELSE 549* 550* 2-by-2 pivot block D(k): columns k and k-1 now hold 551* 552* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 553* 554* where U(k) and U(k-1) are the k-th and (k-1)-th columns 555* of U 556* 557* Perform a rank-2 update of A(1:k-2,1:k-2) as 558* 559* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T 560* = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T 561* 562* and store L(k) and L(k+1) in columns k and k+1 563* 564 IF( K.GT.2 ) THEN 565* 566 D12 = A( K-1, K ) 567 D22 = A( K-1, K-1 ) / D12 568 D11 = A( K, K ) / D12 569 T = ONE / ( D11*D22-ONE ) 570* 571 DO 30 J = K - 2, 1, -1 572* 573 WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) 574 WK = T*( D22*A( J, K )-A( J, K-1 ) ) 575* 576 DO 20 I = J, 1, -1 577 A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - 578 $ ( A( I, K-1 ) / D12 )*WKM1 579 20 CONTINUE 580* 581* Store U(k) and U(k-1) in cols k and k-1 for row J 582* 583 A( J, K ) = WK / D12 584 A( J, K-1 ) = WKM1 / D12 585* 586 30 CONTINUE 587* 588 END IF 589* 590* Copy superdiagonal elements of D(K) to E(K) and 591* ZERO out superdiagonal entry of A 592* 593 E( K ) = A( K-1, K ) 594 E( K-1 ) = ZERO 595 A( K-1, K ) = ZERO 596* 597 END IF 598* 599* End column K is nonsingular 600* 601 END IF 602* 603* Store details of the interchanges in IPIV 604* 605 IF( KSTEP.EQ.1 ) THEN 606 IPIV( K ) = KP 607 ELSE 608 IPIV( K ) = -P 609 IPIV( K-1 ) = -KP 610 END IF 611* 612* Decrease K and return to the start of the main loop 613* 614 K = K - KSTEP 615 GO TO 10 616* 617 34 CONTINUE 618* 619 ELSE 620* 621* Factorize A as L*D*L**T using the lower triangle of A 622* 623* Initialize the unused last entry of the subdiagonal array E. 624* 625 E( N ) = ZERO 626* 627* K is the main loop index, increasing from 1 to N in steps of 628* 1 or 2 629* 630 K = 1 631 40 CONTINUE 632* 633* If K > N, exit from loop 634* 635 IF( K.GT.N ) 636 $ GO TO 64 637 KSTEP = 1 638 P = K 639* 640* Determine rows and columns to be interchanged and whether 641* a 1-by-1 or 2-by-2 pivot block will be used 642* 643 ABSAKK = ABS( A( K, K ) ) 644* 645* IMAX is the row-index of the largest off-diagonal element in 646* column K, and COLMAX is its absolute value. 647* Determine both COLMAX and IMAX. 648* 649 IF( K.LT.N ) THEN 650 IMAX = K + IDAMAX( N-K, A( K+1, K ), 1 ) 651 COLMAX = ABS( A( IMAX, K ) ) 652 ELSE 653 COLMAX = ZERO 654 END IF 655* 656 IF( ( MAX( ABSAKK, COLMAX ).EQ.ZERO ) ) THEN 657* 658* Column K is zero or underflow: set INFO and continue 659* 660 IF( INFO.EQ.0 ) 661 $ INFO = K 662 KP = K 663* 664* Set E( K ) to zero 665* 666 IF( K.LT.N ) 667 $ E( K ) = ZERO 668* 669 ELSE 670* 671* Test for interchange 672* 673* Equivalent to testing for (used to handle NaN and Inf) 674* ABSAKK.GE.ALPHA*COLMAX 675* 676 IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN 677* 678* no interchange, use 1-by-1 pivot block 679* 680 KP = K 681* 682 ELSE 683* 684 DONE = .FALSE. 685* 686* Loop until pivot found 687* 688 42 CONTINUE 689* 690* Begin pivot search loop body 691* 692* JMAX is the column-index of the largest off-diagonal 693* element in row IMAX, and ROWMAX is its absolute value. 694* Determine both ROWMAX and JMAX. 695* 696 IF( IMAX.NE.K ) THEN 697 JMAX = K - 1 + IDAMAX( IMAX-K, A( IMAX, K ), LDA ) 698 ROWMAX = ABS( A( IMAX, JMAX ) ) 699 ELSE 700 ROWMAX = ZERO 701 END IF 702* 703 IF( IMAX.LT.N ) THEN 704 ITEMP = IMAX + IDAMAX( N-IMAX, A( IMAX+1, IMAX ), 705 $ 1 ) 706 DTEMP = ABS( A( ITEMP, IMAX ) ) 707 IF( DTEMP.GT.ROWMAX ) THEN 708 ROWMAX = DTEMP 709 JMAX = ITEMP 710 END IF 711 END IF 712* 713* Equivalent to testing for (used to handle NaN and Inf) 714* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX 715* 716 IF( .NOT.( ABS( A( IMAX, IMAX ) ).LT.ALPHA*ROWMAX ) ) 717 $ THEN 718* 719* interchange rows and columns K and IMAX, 720* use 1-by-1 pivot block 721* 722 KP = IMAX 723 DONE = .TRUE. 724* 725* Equivalent to testing for ROWMAX .EQ. COLMAX, 726* used to handle NaN and Inf 727* 728 ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN 729* 730* interchange rows and columns K+1 and IMAX, 731* use 2-by-2 pivot block 732* 733 KP = IMAX 734 KSTEP = 2 735 DONE = .TRUE. 736 ELSE 737* 738* Pivot NOT found, set variables and repeat 739* 740 P = IMAX 741 COLMAX = ROWMAX 742 IMAX = JMAX 743 END IF 744* 745* End pivot search loop body 746* 747 IF( .NOT. DONE ) GOTO 42 748* 749 END IF 750* 751* Swap TWO rows and TWO columns 752* 753* First swap 754* 755 IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN 756* 757* Interchange rows and column K and P in the trailing 758* submatrix A(k:n,k:n) if we have a 2-by-2 pivot 759* 760 IF( P.LT.N ) 761 $ CALL DSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 ) 762 IF( P.GT.(K+1) ) 763 $ CALL DSWAP( P-K-1, A( K+1, K ), 1, A( P, K+1 ), LDA ) 764 T = A( K, K ) 765 A( K, K ) = A( P, P ) 766 A( P, P ) = T 767* 768* Convert lower triangle of A into L form by applying 769* the interchanges in columns 1:k-1. 770* 771 IF ( K.GT.1 ) 772 $ CALL DSWAP( K-1, A( K, 1 ), LDA, A( P, 1 ), LDA ) 773* 774 END IF 775* 776* Second swap 777* 778 KK = K + KSTEP - 1 779 IF( KP.NE.KK ) THEN 780* 781* Interchange rows and columns KK and KP in the trailing 782* submatrix A(k:n,k:n) 783* 784 IF( KP.LT.N ) 785 $ CALL DSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 ) 786 IF( ( KK.LT.N ) .AND. ( KP.GT.(KK+1) ) ) 787 $ CALL DSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ), 788 $ LDA ) 789 T = A( KK, KK ) 790 A( KK, KK ) = A( KP, KP ) 791 A( KP, KP ) = T 792 IF( KSTEP.EQ.2 ) THEN 793 T = A( K+1, K ) 794 A( K+1, K ) = A( KP, K ) 795 A( KP, K ) = T 796 END IF 797* 798* Convert lower triangle of A into L form by applying 799* the interchanges in columns 1:k-1. 800* 801 IF ( K.GT.1 ) 802 $ CALL DSWAP( K-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA ) 803* 804 END IF 805* 806* Update the trailing submatrix 807* 808 IF( KSTEP.EQ.1 ) THEN 809* 810* 1-by-1 pivot block D(k): column k now holds 811* 812* W(k) = L(k)*D(k) 813* 814* where L(k) is the k-th column of L 815* 816 IF( K.LT.N ) THEN 817* 818* Perform a rank-1 update of A(k+1:n,k+1:n) and 819* store L(k) in column k 820* 821 IF( ABS( A( K, K ) ).GE.SFMIN ) THEN 822* 823* Perform a rank-1 update of A(k+1:n,k+1:n) as 824* A := A - L(k)*D(k)*L(k)**T 825* = A - W(k)*(1/D(k))*W(k)**T 826* 827 D11 = ONE / A( K, K ) 828 CALL DSYR( UPLO, N-K, -D11, A( K+1, K ), 1, 829 $ A( K+1, K+1 ), LDA ) 830* 831* Store L(k) in column k 832* 833 CALL DSCAL( N-K, D11, A( K+1, K ), 1 ) 834 ELSE 835* 836* Store L(k) in column k 837* 838 D11 = A( K, K ) 839 DO 46 II = K + 1, N 840 A( II, K ) = A( II, K ) / D11 841 46 CONTINUE 842* 843* Perform a rank-1 update of A(k+1:n,k+1:n) as 844* A := A - L(k)*D(k)*L(k)**T 845* = A - W(k)*(1/D(k))*W(k)**T 846* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T 847* 848 CALL DSYR( UPLO, N-K, -D11, A( K+1, K ), 1, 849 $ A( K+1, K+1 ), LDA ) 850 END IF 851* 852* Store the subdiagonal element of D in array E 853* 854 E( K ) = ZERO 855* 856 END IF 857* 858 ELSE 859* 860* 2-by-2 pivot block D(k): columns k and k+1 now hold 861* 862* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) 863* 864* where L(k) and L(k+1) are the k-th and (k+1)-th columns 865* of L 866* 867* 868* Perform a rank-2 update of A(k+2:n,k+2:n) as 869* 870* A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T 871* = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T 872* 873* and store L(k) and L(k+1) in columns k and k+1 874* 875 IF( K.LT.N-1 ) THEN 876* 877 D21 = A( K+1, K ) 878 D11 = A( K+1, K+1 ) / D21 879 D22 = A( K, K ) / D21 880 T = ONE / ( D11*D22-ONE ) 881* 882 DO 60 J = K + 2, N 883* 884* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J 885* 886 WK = T*( D11*A( J, K )-A( J, K+1 ) ) 887 WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) 888* 889* Perform a rank-2 update of A(k+2:n,k+2:n) 890* 891 DO 50 I = J, N 892 A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - 893 $ ( A( I, K+1 ) / D21 )*WKP1 894 50 CONTINUE 895* 896* Store L(k) and L(k+1) in cols k and k+1 for row J 897* 898 A( J, K ) = WK / D21 899 A( J, K+1 ) = WKP1 / D21 900* 901 60 CONTINUE 902* 903 END IF 904* 905* Copy subdiagonal elements of D(K) to E(K) and 906* ZERO out subdiagonal entry of A 907* 908 E( K ) = A( K+1, K ) 909 E( K+1 ) = ZERO 910 A( K+1, K ) = ZERO 911* 912 END IF 913* 914* End column K is nonsingular 915* 916 END IF 917* 918* Store details of the interchanges in IPIV 919* 920 IF( KSTEP.EQ.1 ) THEN 921 IPIV( K ) = KP 922 ELSE 923 IPIV( K ) = -P 924 IPIV( K+1 ) = -KP 925 END IF 926* 927* Increase K and return to the start of the main loop 928* 929 K = K + KSTEP 930 GO TO 40 931* 932 64 CONTINUE 933* 934 END IF 935* 936 RETURN 937* 938* End of DSYTF2_RK 939* 940 END 941