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*> \date November 2013 165* 166*> \ingroup realSYcomputational 167* 168*> \par Contributors: 169* ================== 170*> 171*> \verbatim 172*> 173*> November 2013, Igor Kozachenko, 174*> Computer Science Division, 175*> University of California, Berkeley 176*> 177*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, 178*> School of Mathematics, 179*> University of Manchester 180*> 181*> \endverbatim 182* 183* ===================================================================== 184 SUBROUTINE SLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, 185 $ INFO ) 186* 187* -- LAPACK computational routine (version 3.5.0) -- 188* -- LAPACK is a software package provided by Univ. of Tennessee, -- 189* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 190* November 2013 191* 192* .. Scalar Arguments .. 193 CHARACTER UPLO 194 INTEGER INFO, KB, LDA, LDW, N, NB 195* .. 196* .. Array Arguments .. 197 INTEGER IPIV( * ) 198 REAL A( LDA, * ), W( LDW, * ) 199* .. 200* 201* ===================================================================== 202* 203* .. Parameters .. 204 REAL ZERO, ONE 205 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 206 REAL EIGHT, SEVTEN 207 PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 ) 208* .. 209* .. Local Scalars .. 210 LOGICAL DONE 211 INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, JP1, JP2, K, KK, 212 $ KW, KKW, KP, KSTEP, P, II 213 214 REAL ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, 215 $ STEMP, R1, ROWMAX, T, SFMIN 216* .. 217* .. External Functions .. 218 LOGICAL LSAME 219 INTEGER ISAMAX 220 REAL SLAMCH 221 EXTERNAL LSAME, ISAMAX, SLAMCH 222* .. 223* .. External Subroutines .. 224 EXTERNAL SCOPY, SGEMM, SGEMV, SSCAL, SSWAP 225* .. 226* .. Intrinsic Functions .. 227 INTRINSIC ABS, MAX, MIN, SQRT 228* .. 229* .. Executable Statements .. 230* 231 INFO = 0 232* 233* Initialize ALPHA for use in choosing pivot block size. 234* 235 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 236* 237* Compute machine safe minimum 238* 239 SFMIN = SLAMCH( 'S' ) 240* 241 IF( LSAME( UPLO, 'U' ) ) THEN 242* 243* Factorize the trailing columns of A using the upper triangle 244* of A and working backwards, and compute the matrix W = U12*D 245* for use in updating A11 246* 247* K is the main loop index, decreasing from N in steps of 1 or 2 248* 249 K = N 250 10 CONTINUE 251* 252* KW is the column of W which corresponds to column K of A 253* 254 KW = NB + K - N 255* 256* Exit from loop 257* 258 IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 ) 259 $ GO TO 30 260* 261 KSTEP = 1 262 P = K 263* 264* Copy column K of A to column KW of W and update it 265* 266 CALL SCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 ) 267 IF( K.LT.N ) 268 $ CALL SGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ), 269 $ LDA, W( K, KW+1 ), LDW, ONE, W( 1, KW ), 1 ) 270* 271* Determine rows and columns to be interchanged and whether 272* a 1-by-1 or 2-by-2 pivot block will be used 273* 274 ABSAKK = ABS( W( K, KW ) ) 275* 276* IMAX is the row-index of the largest off-diagonal element in 277* column K, and COLMAX is its absolute value. 278* Determine both COLMAX and IMAX. 279* 280 IF( K.GT.1 ) THEN 281 IMAX = ISAMAX( K-1, W( 1, KW ), 1 ) 282 COLMAX = ABS( W( IMAX, KW ) ) 283 ELSE 284 COLMAX = ZERO 285 END IF 286* 287 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 288* 289* Column K is zero or underflow: set INFO and continue 290* 291 IF( INFO.EQ.0 ) 292 $ INFO = K 293 KP = K 294 CALL SCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 ) 295 ELSE 296* 297* ============================================================ 298* 299* Test for interchange 300* 301* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX 302* (used to handle NaN and Inf) 303* 304 IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN 305* 306* no interchange, use 1-by-1 pivot block 307* 308 KP = K 309* 310 ELSE 311* 312 DONE = .FALSE. 313* 314* Loop until pivot found 315* 316 12 CONTINUE 317* 318* Begin pivot search loop body 319* 320* 321* Copy column IMAX to column KW-1 of W and update it 322* 323 CALL SCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 ) 324 CALL SCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA, 325 $ W( IMAX+1, KW-1 ), 1 ) 326* 327 IF( K.LT.N ) 328 $ CALL SGEMV( 'No transpose', K, N-K, -ONE, 329 $ A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW, 330 $ ONE, W( 1, KW-1 ), 1 ) 331* 332* JMAX is the column-index of the largest off-diagonal 333* element in row IMAX, and ROWMAX is its absolute value. 334* Determine both ROWMAX and JMAX. 335* 336 IF( IMAX.NE.K ) THEN 337 JMAX = IMAX + ISAMAX( K-IMAX, W( IMAX+1, KW-1 ), 338 $ 1 ) 339 ROWMAX = ABS( W( JMAX, KW-1 ) ) 340 ELSE 341 ROWMAX = ZERO 342 END IF 343* 344 IF( IMAX.GT.1 ) THEN 345 ITEMP = ISAMAX( IMAX-1, W( 1, KW-1 ), 1 ) 346 STEMP = ABS( W( ITEMP, KW-1 ) ) 347 IF( STEMP.GT.ROWMAX ) THEN 348 ROWMAX = STEMP 349 JMAX = ITEMP 350 END IF 351 END IF 352* 353* Equivalent to testing for 354* ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX 355* (used to handle NaN and Inf) 356* 357 IF( .NOT.(ABS( W( IMAX, KW-1 ) ).LT.ALPHA*ROWMAX ) ) 358 $ THEN 359* 360* interchange rows and columns K and IMAX, 361* use 1-by-1 pivot block 362* 363 KP = IMAX 364* 365* copy column KW-1 of W to column KW of W 366* 367 CALL SCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 ) 368* 369 DONE = .TRUE. 370* 371* Equivalent to testing for ROWMAX.EQ.COLMAX, 372* (used to handle NaN and Inf) 373* 374 ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) ) 375 $ THEN 376* 377* interchange rows and columns K-1 and IMAX, 378* use 2-by-2 pivot block 379* 380 KP = IMAX 381 KSTEP = 2 382 DONE = .TRUE. 383 ELSE 384* 385* Pivot not found: set params and repeat 386* 387 P = IMAX 388 COLMAX = ROWMAX 389 IMAX = JMAX 390* 391* Copy updated JMAXth (next IMAXth) column to Kth of W 392* 393 CALL SCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 ) 394* 395 END IF 396* 397* End pivot search loop body 398* 399 IF( .NOT. DONE ) GOTO 12 400* 401 END IF 402* 403* ============================================================ 404* 405 KK = K - KSTEP + 1 406* 407* KKW is the column of W which corresponds to column KK of A 408* 409 KKW = NB + KK - N 410* 411 IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN 412* 413* Copy non-updated column K to column P 414* 415 CALL SCOPY( K-P, A( P+1, K ), 1, A( P, P+1 ), LDA ) 416 CALL SCOPY( P, A( 1, K ), 1, A( 1, P ), 1 ) 417* 418* Interchange rows K and P in last N-K+1 columns of A 419* and last N-K+2 columns of W 420* 421 CALL SSWAP( N-K+1, A( K, K ), LDA, A( P, K ), LDA ) 422 CALL SSWAP( N-KK+1, W( K, KKW ), LDW, W( P, KKW ), LDW ) 423 END IF 424* 425* Updated column KP is already stored in column KKW of W 426* 427 IF( KP.NE.KK ) THEN 428* 429* Copy non-updated column KK to column KP 430* 431 A( KP, K ) = A( KK, K ) 432 CALL SCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ), 433 $ LDA ) 434 CALL SCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 ) 435* 436* Interchange rows KK and KP in last N-KK+1 columns 437* of A and W 438* 439 CALL SSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA ) 440 CALL SSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ), 441 $ LDW ) 442 END IF 443* 444 IF( KSTEP.EQ.1 ) THEN 445* 446* 1-by-1 pivot block D(k): column KW of W now holds 447* 448* W(k) = U(k)*D(k) 449* 450* where U(k) is the k-th column of U 451* 452* Store U(k) in column k of A 453* 454 CALL SCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 ) 455 IF( K.GT.1 ) THEN 456 IF( ABS( A( K, K ) ).GE.SFMIN ) THEN 457 R1 = ONE / A( K, K ) 458 CALL SSCAL( K-1, R1, A( 1, K ), 1 ) 459 ELSE IF( A( K, K ).NE.ZERO ) THEN 460 DO 14 II = 1, K - 1 461 A( II, K ) = A( II, K ) / A( K, K ) 462 14 CONTINUE 463 END IF 464 END IF 465* 466 ELSE 467* 468* 2-by-2 pivot block D(k): columns KW and KW-1 of W now 469* hold 470* 471* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 472* 473* where U(k) and U(k-1) are the k-th and (k-1)-th columns 474* of U 475* 476 IF( K.GT.2 ) THEN 477* 478* Store U(k) and U(k-1) in columns k and k-1 of A 479* 480 D12 = W( K-1, KW ) 481 D11 = W( K, KW ) / D12 482 D22 = W( K-1, KW-1 ) / D12 483 T = ONE / ( D11*D22-ONE ) 484 DO 20 J = 1, K - 2 485 A( J, K-1 ) = T*( (D11*W( J, KW-1 )-W( J, KW ) ) / 486 $ D12 ) 487 A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) / 488 $ D12 ) 489 20 CONTINUE 490 END IF 491* 492* Copy D(k) to A 493* 494 A( K-1, K-1 ) = W( K-1, KW-1 ) 495 A( K-1, K ) = W( K-1, KW ) 496 A( K, K ) = W( K, KW ) 497 END IF 498 END IF 499* 500* Store details of the interchanges in IPIV 501* 502 IF( KSTEP.EQ.1 ) THEN 503 IPIV( K ) = KP 504 ELSE 505 IPIV( K ) = -P 506 IPIV( K-1 ) = -KP 507 END IF 508* 509* Decrease K and return to the start of the main loop 510* 511 K = K - KSTEP 512 GO TO 10 513* 514 30 CONTINUE 515* 516* Update the upper triangle of A11 (= A(1:k,1:k)) as 517* 518* A11 := A11 - U12*D*U12**T = A11 - U12*W**T 519* 520* computing blocks of NB columns at a time 521* 522 DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB 523 JB = MIN( NB, K-J+1 ) 524* 525* Update the upper triangle of the diagonal block 526* 527 DO 40 JJ = J, J + JB - 1 528 CALL SGEMV( 'No transpose', JJ-J+1, N-K, -ONE, 529 $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE, 530 $ A( J, JJ ), 1 ) 531 40 CONTINUE 532* 533* Update the rectangular superdiagonal block 534* 535 IF( J.GE.2 ) 536 $ CALL SGEMM( 'No transpose', 'Transpose', J-1, JB, 537 $ N-K, -ONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, 538 $ ONE, A( 1, J ), LDA ) 539 50 CONTINUE 540* 541* Put U12 in standard form by partially undoing the interchanges 542* in columns k+1:n 543* 544 J = K + 1 545 60 CONTINUE 546* 547 KSTEP = 1 548 JP1 = 1 549 JJ = J 550 JP2 = IPIV( J ) 551 IF( JP2.LT.0 ) THEN 552 JP2 = -JP2 553 J = J + 1 554 JP1 = -IPIV( J ) 555 KSTEP = 2 556 END IF 557* 558 J = J + 1 559 IF( JP2.NE.JJ .AND. J.LE.N ) 560 $ CALL SSWAP( N-J+1, A( JP2, J ), LDA, A( JJ, J ), LDA ) 561 JJ = J - 1 562 IF( JP1.NE.JJ .AND. KSTEP.EQ.2 ) 563 $ CALL SSWAP( N-J+1, A( JP1, J ), LDA, A( JJ, J ), LDA ) 564 IF( J.LE.N ) 565 $ GO TO 60 566* 567* Set KB to the number of columns factorized 568* 569 KB = N - K 570* 571 ELSE 572* 573* Factorize the leading columns of A using the lower triangle 574* of A and working forwards, and compute the matrix W = L21*D 575* for use in updating A22 576* 577* K is the main loop index, increasing from 1 in steps of 1 or 2 578* 579 K = 1 580 70 CONTINUE 581* 582* Exit from loop 583* 584 IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N ) 585 $ GO TO 90 586* 587 KSTEP = 1 588 P = K 589* 590* Copy column K of A to column K of W and update it 591* 592 CALL SCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 ) 593 IF( K.GT.1 ) 594 $ CALL SGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ), 595 $ LDA, W( K, 1 ), LDW, ONE, W( K, K ), 1 ) 596* 597* Determine rows and columns to be interchanged and whether 598* a 1-by-1 or 2-by-2 pivot block will be used 599* 600 ABSAKK = ABS( W( K, K ) ) 601* 602* IMAX is the row-index of the largest off-diagonal element in 603* column K, and COLMAX is its absolute value. 604* Determine both COLMAX and IMAX. 605* 606 IF( K.LT.N ) THEN 607 IMAX = K + ISAMAX( N-K, W( K+1, K ), 1 ) 608 COLMAX = ABS( W( IMAX, K ) ) 609 ELSE 610 COLMAX = ZERO 611 END IF 612* 613 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 614* 615* Column K is zero or underflow: set INFO and continue 616* 617 IF( INFO.EQ.0 ) 618 $ INFO = K 619 KP = K 620 CALL SCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 ) 621 ELSE 622* 623* ============================================================ 624* 625* Test for interchange 626* 627* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX 628* (used to handle NaN and Inf) 629* 630 IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN 631* 632* no interchange, use 1-by-1 pivot block 633* 634 KP = K 635* 636 ELSE 637* 638 DONE = .FALSE. 639* 640* Loop until pivot found 641* 642 72 CONTINUE 643* 644* Begin pivot search loop body 645* 646* 647* Copy column IMAX to column K+1 of W and update it 648* 649 CALL SCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1) 650 CALL SCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, 651 $ W( IMAX, K+1 ), 1 ) 652 IF( K.GT.1 ) 653 $ CALL SGEMV( 'No transpose', N-K+1, K-1, -ONE, 654 $ A( K, 1 ), LDA, W( IMAX, 1 ), LDW, 655 $ ONE, W( K, K+1 ), 1 ) 656* 657* JMAX is the column-index of the largest off-diagonal 658* element in row IMAX, and ROWMAX is its absolute value. 659* Determine both ROWMAX and JMAX. 660* 661 IF( IMAX.NE.K ) THEN 662 JMAX = K - 1 + ISAMAX( IMAX-K, W( K, K+1 ), 1 ) 663 ROWMAX = ABS( W( JMAX, K+1 ) ) 664 ELSE 665 ROWMAX = ZERO 666 END IF 667* 668 IF( IMAX.LT.N ) THEN 669 ITEMP = IMAX + ISAMAX( N-IMAX, W( IMAX+1, K+1 ), 1) 670 STEMP = ABS( W( ITEMP, K+1 ) ) 671 IF( STEMP.GT.ROWMAX ) THEN 672 ROWMAX = STEMP 673 JMAX = ITEMP 674 END IF 675 END IF 676* 677* Equivalent to testing for 678* ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX 679* (used to handle NaN and Inf) 680* 681 IF( .NOT.( ABS( W( IMAX, K+1 ) ).LT.ALPHA*ROWMAX ) ) 682 $ THEN 683* 684* interchange rows and columns K and IMAX, 685* use 1-by-1 pivot block 686* 687 KP = IMAX 688* 689* copy column K+1 of W to column K of W 690* 691 CALL SCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 ) 692* 693 DONE = .TRUE. 694* 695* Equivalent to testing for ROWMAX.EQ.COLMAX, 696* (used to handle NaN and Inf) 697* 698 ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) ) 699 $ THEN 700* 701* interchange rows and columns K+1 and IMAX, 702* use 2-by-2 pivot block 703* 704 KP = IMAX 705 KSTEP = 2 706 DONE = .TRUE. 707 ELSE 708* 709* Pivot not found: set params and repeat 710* 711 P = IMAX 712 COLMAX = ROWMAX 713 IMAX = JMAX 714* 715* Copy updated JMAXth (next IMAXth) column to Kth of W 716* 717 CALL SCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 ) 718* 719 END IF 720* 721* End pivot search loop body 722* 723 IF( .NOT. DONE ) GOTO 72 724* 725 END IF 726* 727* ============================================================ 728* 729 KK = K + KSTEP - 1 730* 731 IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN 732* 733* Copy non-updated column K to column P 734* 735 CALL SCOPY( P-K, A( K, K ), 1, A( P, K ), LDA ) 736 CALL SCOPY( N-P+1, A( P, K ), 1, A( P, P ), 1 ) 737* 738* Interchange rows K and P in first K columns of A 739* and first K+1 columns of W 740* 741 CALL SSWAP( K, A( K, 1 ), LDA, A( P, 1 ), LDA ) 742 CALL SSWAP( KK, W( K, 1 ), LDW, W( P, 1 ), LDW ) 743 END IF 744* 745* Updated column KP is already stored in column KK of W 746* 747 IF( KP.NE.KK ) THEN 748* 749* Copy non-updated column KK to column KP 750* 751 A( KP, K ) = A( KK, K ) 752 CALL SCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA ) 753 CALL SCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 ) 754* 755* Interchange rows KK and KP in first KK columns of A and W 756* 757 CALL SSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA ) 758 CALL SSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW ) 759 END IF 760* 761 IF( KSTEP.EQ.1 ) THEN 762* 763* 1-by-1 pivot block D(k): column k of W now holds 764* 765* W(k) = L(k)*D(k) 766* 767* where L(k) is the k-th column of L 768* 769* Store L(k) in column k of A 770* 771 CALL SCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 ) 772 IF( K.LT.N ) THEN 773 IF( ABS( A( K, K ) ).GE.SFMIN ) THEN 774 R1 = ONE / A( K, K ) 775 CALL SSCAL( N-K, R1, A( K+1, K ), 1 ) 776 ELSE IF( A( K, K ).NE.ZERO ) THEN 777 DO 74 II = K + 1, N 778 A( II, K ) = A( II, K ) / A( K, K ) 779 74 CONTINUE 780 END IF 781 END IF 782* 783 ELSE 784* 785* 2-by-2 pivot block D(k): columns k and k+1 of W now hold 786* 787* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) 788* 789* where L(k) and L(k+1) are the k-th and (k+1)-th columns 790* of L 791* 792 IF( K.LT.N-1 ) THEN 793* 794* Store L(k) and L(k+1) in columns k and k+1 of A 795* 796 D21 = W( K+1, K ) 797 D11 = W( K+1, K+1 ) / D21 798 D22 = W( K, K ) / D21 799 T = ONE / ( D11*D22-ONE ) 800 DO 80 J = K + 2, N 801 A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) / 802 $ D21 ) 803 A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) / 804 $ D21 ) 805 80 CONTINUE 806 END IF 807* 808* Copy D(k) to A 809* 810 A( K, K ) = W( K, K ) 811 A( K+1, K ) = W( K+1, K ) 812 A( K+1, K+1 ) = W( K+1, K+1 ) 813 END IF 814 END IF 815* 816* Store details of the interchanges in IPIV 817* 818 IF( KSTEP.EQ.1 ) THEN 819 IPIV( K ) = KP 820 ELSE 821 IPIV( K ) = -P 822 IPIV( K+1 ) = -KP 823 END IF 824* 825* Increase K and return to the start of the main loop 826* 827 K = K + KSTEP 828 GO TO 70 829* 830 90 CONTINUE 831* 832* Update the lower triangle of A22 (= A(k:n,k:n)) as 833* 834* A22 := A22 - L21*D*L21**T = A22 - L21*W**T 835* 836* computing blocks of NB columns at a time 837* 838 DO 110 J = K, N, NB 839 JB = MIN( NB, N-J+1 ) 840* 841* Update the lower triangle of the diagonal block 842* 843 DO 100 JJ = J, J + JB - 1 844 CALL SGEMV( 'No transpose', J+JB-JJ, K-1, -ONE, 845 $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE, 846 $ A( JJ, JJ ), 1 ) 847 100 CONTINUE 848* 849* Update the rectangular subdiagonal block 850* 851 IF( J+JB.LE.N ) 852 $ CALL SGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB, 853 $ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW, 854 $ ONE, A( J+JB, J ), LDA ) 855 110 CONTINUE 856* 857* Put L21 in standard form by partially undoing the interchanges 858* in columns 1:k-1 859* 860 J = K - 1 861 120 CONTINUE 862* 863 KSTEP = 1 864 JP1 = 1 865 JJ = J 866 JP2 = IPIV( J ) 867 IF( JP2.LT.0 ) THEN 868 JP2 = -JP2 869 J = J - 1 870 JP1 = -IPIV( J ) 871 KSTEP = 2 872 END IF 873* 874 J = J - 1 875 IF( JP2.NE.JJ .AND. J.GE.1 ) 876 $ CALL SSWAP( J, A( JP2, 1 ), LDA, A( JJ, 1 ), LDA ) 877 JJ = J + 1 878 IF( JP1.NE.JJ .AND. KSTEP.EQ.2 ) 879 $ CALL SSWAP( J, A( JP1, 1 ), LDA, A( JJ, 1 ), LDA ) 880 IF( J.GE.1 ) 881 $ GO TO 120 882* 883* Set KB to the number of columns factorized 884* 885 KB = K - 1 886* 887 END IF 888 RETURN 889* 890* End of SLASYF_ROOK 891* 892 END 893