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