1*> \brief \b SSYTRF_AA_2STAGE 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SSYTRF_AA_2STAGE + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytrf_aa_2stage.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytrf_aa_2stage.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytrf_aa_2stage.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, 22* IPIV2, WORK, LWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER UPLO 26* INTEGER N, LDA, LTB, LWORK, INFO 27* .. 28* .. Array Arguments .. 29* INTEGER IPIV( * ), IPIV2( * ) 30* REAL A( LDA, * ), TB( * ), WORK( * ) 31* .. 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> SSYTRF_AA_2STAGE computes the factorization of a real symmetric matrix A 39*> using the Aasen's algorithm. The form of the factorization is 40*> 41*> A = U**T*T*U or A = L*T*L**T 42*> 43*> where U (or L) is a product of permutation and unit upper (lower) 44*> triangular matrices, and T is a symmetric band matrix with the 45*> bandwidth of NB (NB is internally selected and stored in TB( 1 ), and T is 46*> LU factorized with partial pivoting). 47*> 48*> This is the blocked version of the algorithm, calling Level 3 BLAS. 49*> \endverbatim 50* 51* Arguments: 52* ========== 53* 54*> \param[in] UPLO 55*> \verbatim 56*> UPLO is CHARACTER*1 57*> = 'U': Upper triangle of A is stored; 58*> = 'L': Lower triangle of A is stored. 59*> \endverbatim 60*> 61*> \param[in] N 62*> \verbatim 63*> N is INTEGER 64*> The order of the matrix A. N >= 0. 65*> \endverbatim 66*> 67*> \param[in,out] A 68*> \verbatim 69*> A is REAL array, dimension (LDA,N) 70*> On entry, the symmetric matrix A. If UPLO = 'U', the leading 71*> N-by-N upper triangular part of A contains the upper 72*> triangular part of the matrix A, and the strictly lower 73*> triangular part of A is not referenced. If UPLO = 'L', the 74*> leading N-by-N lower triangular part of A contains the lower 75*> triangular part of the matrix A, and the strictly upper 76*> triangular part of A is not referenced. 77*> 78*> On exit, L is stored below (or above) the subdiaonal blocks, 79*> when UPLO is 'L' (or 'U'). 80*> \endverbatim 81*> 82*> \param[in] LDA 83*> \verbatim 84*> LDA is INTEGER 85*> The leading dimension of the array A. LDA >= max(1,N). 86*> \endverbatim 87*> 88*> \param[out] TB 89*> \verbatim 90*> TB is REAL array, dimension (LTB) 91*> On exit, details of the LU factorization of the band matrix. 92*> \endverbatim 93*> 94*> \param[in] LTB 95*> \verbatim 96*> LTB is INTEGER 97*> The size of the array TB. LTB >= 4*N, internally 98*> used to select NB such that LTB >= (3*NB+1)*N. 99*> 100*> If LTB = -1, then a workspace query is assumed; the 101*> routine only calculates the optimal size of LTB, 102*> returns this value as the first entry of TB, and 103*> no error message related to LTB is issued by XERBLA. 104*> \endverbatim 105*> 106*> \param[out] IPIV 107*> \verbatim 108*> IPIV is INTEGER array, dimension (N) 109*> On exit, it contains the details of the interchanges, i.e., 110*> the row and column k of A were interchanged with the 111*> row and column IPIV(k). 112*> \endverbatim 113*> 114*> \param[out] IPIV2 115*> \verbatim 116*> IPIV2 is INTEGER array, dimension (N) 117*> On exit, it contains the details of the interchanges, i.e., 118*> the row and column k of T were interchanged with the 119*> row and column IPIV(k). 120*> \endverbatim 121*> 122*> \param[out] WORK 123*> \verbatim 124*> WORK is REAL workspace of size LWORK 125*> \endverbatim 126*> 127*> \param[in] LWORK 128*> \verbatim 129*> LWORK is INTEGER 130*> The size of WORK. LWORK >= N, internally used to select NB 131*> such that LWORK >= N*NB. 132*> 133*> If LWORK = -1, then a workspace query is assumed; the 134*> routine only calculates the optimal size of the WORK array, 135*> returns this value as the first entry of the WORK array, and 136*> no error message related to LWORK is issued by XERBLA. 137*> \endverbatim 138*> 139*> \param[out] INFO 140*> \verbatim 141*> INFO is INTEGER 142*> = 0: successful exit 143*> < 0: if INFO = -i, the i-th argument had an illegal value. 144*> > 0: if INFO = i, band LU factorization failed on i-th column 145*> \endverbatim 146* 147* Authors: 148* ======== 149* 150*> \author Univ. of Tennessee 151*> \author Univ. of California Berkeley 152*> \author Univ. of Colorado Denver 153*> \author NAG Ltd. 154* 155*> \date November 2017 156* 157*> \ingroup realSYcomputational 158* 159* ===================================================================== 160 SUBROUTINE SSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, 161 $ IPIV2, WORK, LWORK, INFO ) 162* 163* -- LAPACK computational routine (version 3.8.0) -- 164* -- LAPACK is a software package provided by Univ. of Tennessee, -- 165* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 166* November 2017 167* 168 IMPLICIT NONE 169* 170* .. Scalar Arguments .. 171 CHARACTER UPLO 172 INTEGER N, LDA, LTB, LWORK, INFO 173* .. 174* .. Array Arguments .. 175 INTEGER IPIV( * ), IPIV2( * ) 176 REAL A( LDA, * ), TB( * ), WORK( * ) 177* .. 178* 179* ===================================================================== 180* .. Parameters .. 181 REAL ZERO, ONE 182 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 183* 184* .. Local Scalars .. 185 LOGICAL UPPER, TQUERY, WQUERY 186 INTEGER I, J, K, I1, I2, TD 187 INTEGER LDTB, NB, KB, JB, NT, IINFO 188 REAL PIV 189* .. 190* .. External Functions .. 191 LOGICAL LSAME 192 INTEGER ILAENV 193 EXTERNAL LSAME, ILAENV 194* .. 195* .. External Subroutines .. 196 EXTERNAL XERBLA, SCOPY, SLACPY, 197 $ SLASET, SGBTRF, SGEMM, SGETRF, 198 $ SSYGST, SSWAP, STRSM 199* .. 200* .. Intrinsic Functions .. 201 INTRINSIC MIN, MAX 202* .. 203* .. Executable Statements .. 204* 205* Test the input parameters. 206* 207 INFO = 0 208 UPPER = LSAME( UPLO, 'U' ) 209 WQUERY = ( LWORK.EQ.-1 ) 210 TQUERY = ( LTB.EQ.-1 ) 211 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 212 INFO = -1 213 ELSE IF( N.LT.0 ) THEN 214 INFO = -2 215 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 216 INFO = -4 217 ELSE IF ( LTB .LT. 4*N .AND. .NOT.TQUERY ) THEN 218 INFO = -6 219 ELSE IF ( LWORK .LT. N .AND. .NOT.WQUERY ) THEN 220 INFO = -10 221 END IF 222* 223 IF( INFO.NE.0 ) THEN 224 CALL XERBLA( 'SSYTRF_AA_2STAGE', -INFO ) 225 RETURN 226 END IF 227* 228* Answer the query 229* 230 NB = ILAENV( 1, 'SSYTRF_AA_2STAGE', UPLO, N, -1, -1, -1 ) 231 IF( INFO.EQ.0 ) THEN 232 IF( TQUERY ) THEN 233 TB( 1 ) = (3*NB+1)*N 234 END IF 235 IF( WQUERY ) THEN 236 WORK( 1 ) = N*NB 237 END IF 238 END IF 239 IF( TQUERY .OR. WQUERY ) THEN 240 RETURN 241 END IF 242* 243* Quick return 244* 245 IF ( N.EQ.0 ) THEN 246 RETURN 247 ENDIF 248* 249* Determine the number of the block size 250* 251 LDTB = LTB/N 252 IF( LDTB .LT. 3*NB+1 ) THEN 253 NB = (LDTB-1)/3 254 END IF 255 IF( LWORK .LT. NB*N ) THEN 256 NB = LWORK/N 257 END IF 258* 259* Determine the number of the block columns 260* 261 NT = (N+NB-1)/NB 262 TD = 2*NB 263 KB = MIN(NB, N) 264* 265* Initialize vectors/matrices 266* 267 DO J = 1, KB 268 IPIV( J ) = J 269 END DO 270* 271* Save NB 272* 273 TB( 1 ) = NB 274* 275 IF( UPPER ) THEN 276* 277* ..................................................... 278* Factorize A as U**T*D*U using the upper triangle of A 279* ..................................................... 280* 281 DO J = 0, NT-1 282* 283* Generate Jth column of W and H 284* 285 KB = MIN(NB, N-J*NB) 286 DO I = 1, J-1 287 IF( I.EQ.1 ) THEN 288* H(I,J) = T(I,I)*U(I,J) + T(I+1,I)*U(I+1,J) 289 IF( I .EQ. (J-1) ) THEN 290 JB = NB+KB 291 ELSE 292 JB = 2*NB 293 END IF 294 CALL SGEMM( 'NoTranspose', 'NoTranspose', 295 $ NB, KB, JB, 296 $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1, 297 $ A( (I-1)*NB+1, J*NB+1 ), LDA, 298 $ ZERO, WORK( I*NB+1 ), N ) 299 ELSE 300* H(I,J) = T(I,I-1)*U(I-1,J) + T(I,I)*U(I,J) + T(I,I+1)*U(I+1,J) 301 IF( I .EQ. J-1) THEN 302 JB = 2*NB+KB 303 ELSE 304 JB = 3*NB 305 END IF 306 CALL SGEMM( 'NoTranspose', 'NoTranspose', 307 $ NB, KB, JB, 308 $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ), 309 $ LDTB-1, 310 $ A( (I-2)*NB+1, J*NB+1 ), LDA, 311 $ ZERO, WORK( I*NB+1 ), N ) 312 END IF 313 END DO 314* 315* Compute T(J,J) 316* 317 CALL SLACPY( 'Upper', KB, KB, A( J*NB+1, J*NB+1 ), LDA, 318 $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 319 IF( J.GT.1 ) THEN 320* T(J,J) = U(1:J,J)'*H(1:J) 321 CALL SGEMM( 'Transpose', 'NoTranspose', 322 $ KB, KB, (J-1)*NB, 323 $ -ONE, A( 1, J*NB+1 ), LDA, 324 $ WORK( NB+1 ), N, 325 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 326* T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J) 327 CALL SGEMM( 'Transpose', 'NoTranspose', 328 $ KB, NB, KB, 329 $ ONE, A( (J-1)*NB+1, J*NB+1 ), LDA, 330 $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1, 331 $ ZERO, WORK( 1 ), N ) 332 CALL SGEMM( 'NoTranspose', 'NoTranspose', 333 $ KB, KB, NB, 334 $ -ONE, WORK( 1 ), N, 335 $ A( (J-2)*NB+1, J*NB+1 ), LDA, 336 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 337 END IF 338 IF( J.GT.0 ) THEN 339 CALL SSYGST( 1, 'Upper', KB, 340 $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1, 341 $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO ) 342 END IF 343* 344* Expand T(J,J) into full format 345* 346 DO I = 1, KB 347 DO K = I+1, KB 348 TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB ) 349 $ = TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB ) 350 END DO 351 END DO 352* 353 IF( J.LT.NT-1 ) THEN 354 IF( J.GT.0 ) THEN 355* 356* Compute H(J,J) 357* 358 IF( J.EQ.1 ) THEN 359 CALL SGEMM( 'NoTranspose', 'NoTranspose', 360 $ KB, KB, KB, 361 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1, 362 $ A( (J-1)*NB+1, J*NB+1 ), LDA, 363 $ ZERO, WORK( J*NB+1 ), N ) 364 ELSE 365 CALL SGEMM( 'NoTranspose', 'NoTranspose', 366 $ KB, KB, NB+KB, 367 $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ), 368 $ LDTB-1, 369 $ A( (J-2)*NB+1, J*NB+1 ), LDA, 370 $ ZERO, WORK( J*NB+1 ), N ) 371 END IF 372* 373* Update with the previous column 374* 375 CALL SGEMM( 'Transpose', 'NoTranspose', 376 $ NB, N-(J+1)*NB, J*NB, 377 $ -ONE, WORK( NB+1 ), N, 378 $ A( 1, (J+1)*NB+1 ), LDA, 379 $ ONE, A( J*NB+1, (J+1)*NB+1 ), LDA ) 380 END IF 381* 382* Copy panel to workspace to call SGETRF 383* 384 DO K = 1, NB 385 CALL SCOPY( N-(J+1)*NB, 386 $ A( J*NB+K, (J+1)*NB+1 ), LDA, 387 $ WORK( 1+(K-1)*N ), 1 ) 388 END DO 389* 390* Factorize panel 391* 392 CALL SGETRF( N-(J+1)*NB, NB, 393 $ WORK, N, 394 $ IPIV( (J+1)*NB+1 ), IINFO ) 395c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN 396c INFO = IINFO+(J+1)*NB 397c END IF 398* 399* Copy panel back 400* 401 DO K = 1, NB 402 CALL SCOPY( N-(J+1)*NB, 403 $ WORK( 1+(K-1)*N ), 1, 404 $ A( J*NB+K, (J+1)*NB+1 ), LDA ) 405 END DO 406* 407* Compute T(J+1, J), zero out for GEMM update 408* 409 KB = MIN(NB, N-(J+1)*NB) 410 CALL SLASET( 'Full', KB, NB, ZERO, ZERO, 411 $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 ) 412 CALL SLACPY( 'Upper', KB, NB, 413 $ WORK, N, 414 $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 ) 415 IF( J.GT.0 ) THEN 416 CALL STRSM( 'R', 'U', 'N', 'U', KB, NB, ONE, 417 $ A( (J-1)*NB+1, J*NB+1 ), LDA, 418 $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 ) 419 END IF 420* 421* Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM 422* updates 423* 424 DO K = 1, NB 425 DO I = 1, KB 426 TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB ) 427 $ = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB ) 428 END DO 429 END DO 430 CALL SLASET( 'Lower', KB, NB, ZERO, ONE, 431 $ A( J*NB+1, (J+1)*NB+1), LDA ) 432* 433* Apply pivots to trailing submatrix of A 434* 435 DO K = 1, KB 436* > Adjust ipiv 437 IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB 438* 439 I1 = (J+1)*NB+K 440 I2 = IPIV( (J+1)*NB+K ) 441 IF( I1.NE.I2 ) THEN 442* > Apply pivots to previous columns of L 443 CALL SSWAP( K-1, A( (J+1)*NB+1, I1 ), 1, 444 $ A( (J+1)*NB+1, I2 ), 1 ) 445* > Swap A(I1+1:M, I1) with A(I2, I1+1:M) 446 IF( I2.GT.(I1+1) ) 447 $ CALL SSWAP( I2-I1-1, A( I1, I1+1 ), LDA, 448 $ A( I1+1, I2 ), 1 ) 449* > Swap A(I2+1:M, I1) with A(I2+1:M, I2) 450 IF( I2.LT.N ) 451 $ CALL SSWAP( N-I2, A( I1, I2+1 ), LDA, 452 $ A( I2, I2+1 ), LDA ) 453* > Swap A(I1, I1) with A(I2, I2) 454 PIV = A( I1, I1 ) 455 A( I1, I1 ) = A( I2, I2 ) 456 A( I2, I2 ) = PIV 457* > Apply pivots to previous columns of L 458 IF( J.GT.0 ) THEN 459 CALL SSWAP( J*NB, A( 1, I1 ), 1, 460 $ A( 1, I2 ), 1 ) 461 END IF 462 ENDIF 463 END DO 464 END IF 465 END DO 466 ELSE 467* 468* ..................................................... 469* Factorize A as L*D*L**T using the lower triangle of A 470* ..................................................... 471* 472 DO J = 0, NT-1 473* 474* Generate Jth column of W and H 475* 476 KB = MIN(NB, N-J*NB) 477 DO I = 1, J-1 478 IF( I.EQ.1 ) THEN 479* H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)' 480 IF( I .EQ. (J-1) ) THEN 481 JB = NB+KB 482 ELSE 483 JB = 2*NB 484 END IF 485 CALL SGEMM( 'NoTranspose', 'Transpose', 486 $ NB, KB, JB, 487 $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1, 488 $ A( J*NB+1, (I-1)*NB+1 ), LDA, 489 $ ZERO, WORK( I*NB+1 ), N ) 490 ELSE 491* H(I,J) = T(I,I-1)*L(J,I-1)' + T(I,I)*L(J,I)' + T(I,I+1)*L(J,I+1)' 492 IF( I .EQ. J-1) THEN 493 JB = 2*NB+KB 494 ELSE 495 JB = 3*NB 496 END IF 497 CALL SGEMM( 'NoTranspose', 'Transpose', 498 $ NB, KB, JB, 499 $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ), 500 $ LDTB-1, 501 $ A( J*NB+1, (I-2)*NB+1 ), LDA, 502 $ ZERO, WORK( I*NB+1 ), N ) 503 END IF 504 END DO 505* 506* Compute T(J,J) 507* 508 CALL SLACPY( 'Lower', KB, KB, A( J*NB+1, J*NB+1 ), LDA, 509 $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 510 IF( J.GT.1 ) THEN 511* T(J,J) = L(J,1:J)*H(1:J) 512 CALL SGEMM( 'NoTranspose', 'NoTranspose', 513 $ KB, KB, (J-1)*NB, 514 $ -ONE, A( J*NB+1, 1 ), LDA, 515 $ WORK( NB+1 ), N, 516 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 517* T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)' 518 CALL SGEMM( 'NoTranspose', 'NoTranspose', 519 $ KB, NB, KB, 520 $ ONE, A( J*NB+1, (J-1)*NB+1 ), LDA, 521 $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1, 522 $ ZERO, WORK( 1 ), N ) 523 CALL SGEMM( 'NoTranspose', 'Transpose', 524 $ KB, KB, NB, 525 $ -ONE, WORK( 1 ), N, 526 $ A( J*NB+1, (J-2)*NB+1 ), LDA, 527 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 528 END IF 529 IF( J.GT.0 ) THEN 530 CALL SSYGST( 1, 'Lower', KB, 531 $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1, 532 $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO ) 533 END IF 534* 535* Expand T(J,J) into full format 536* 537 DO I = 1, KB 538 DO K = I+1, KB 539 TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB ) 540 $ = TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB ) 541 END DO 542 END DO 543* 544 IF( J.LT.NT-1 ) THEN 545 IF( J.GT.0 ) THEN 546* 547* Compute H(J,J) 548* 549 IF( J.EQ.1 ) THEN 550 CALL SGEMM( 'NoTranspose', 'Transpose', 551 $ KB, KB, KB, 552 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1, 553 $ A( J*NB+1, (J-1)*NB+1 ), LDA, 554 $ ZERO, WORK( J*NB+1 ), N ) 555 ELSE 556 CALL SGEMM( 'NoTranspose', 'Transpose', 557 $ KB, KB, NB+KB, 558 $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ), 559 $ LDTB-1, 560 $ A( J*NB+1, (J-2)*NB+1 ), LDA, 561 $ ZERO, WORK( J*NB+1 ), N ) 562 END IF 563* 564* Update with the previous column 565* 566 CALL SGEMM( 'NoTranspose', 'NoTranspose', 567 $ N-(J+1)*NB, NB, J*NB, 568 $ -ONE, A( (J+1)*NB+1, 1 ), LDA, 569 $ WORK( NB+1 ), N, 570 $ ONE, A( (J+1)*NB+1, J*NB+1 ), LDA ) 571 END IF 572* 573* Factorize panel 574* 575 CALL SGETRF( N-(J+1)*NB, NB, 576 $ A( (J+1)*NB+1, J*NB+1 ), LDA, 577 $ IPIV( (J+1)*NB+1 ), IINFO ) 578c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN 579c INFO = IINFO+(J+1)*NB 580c END IF 581* 582* Compute T(J+1, J), zero out for GEMM update 583* 584 KB = MIN(NB, N-(J+1)*NB) 585 CALL SLASET( 'Full', KB, NB, ZERO, ZERO, 586 $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 ) 587 CALL SLACPY( 'Upper', KB, NB, 588 $ A( (J+1)*NB+1, J*NB+1 ), LDA, 589 $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 ) 590 IF( J.GT.0 ) THEN 591 CALL STRSM( 'R', 'L', 'T', 'U', KB, NB, ONE, 592 $ A( J*NB+1, (J-1)*NB+1 ), LDA, 593 $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 ) 594 END IF 595* 596* Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM 597* updates 598* 599 DO K = 1, NB 600 DO I = 1, KB 601 TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB ) = 602 $ TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB ) 603 END DO 604 END DO 605 CALL SLASET( 'Upper', KB, NB, ZERO, ONE, 606 $ A( (J+1)*NB+1, J*NB+1), LDA ) 607* 608* Apply pivots to trailing submatrix of A 609* 610 DO K = 1, KB 611* > Adjust ipiv 612 IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB 613* 614 I1 = (J+1)*NB+K 615 I2 = IPIV( (J+1)*NB+K ) 616 IF( I1.NE.I2 ) THEN 617* > Apply pivots to previous columns of L 618 CALL SSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA, 619 $ A( I2, (J+1)*NB+1 ), LDA ) 620* > Swap A(I1+1:M, I1) with A(I2, I1+1:M) 621 IF( I2.GT.(I1+1) ) 622 $ CALL SSWAP( I2-I1-1, A( I1+1, I1 ), 1, 623 $ A( I2, I1+1 ), LDA ) 624* > Swap A(I2+1:M, I1) with A(I2+1:M, I2) 625 IF( I2.LT.N ) 626 $ CALL SSWAP( N-I2, A( I2+1, I1 ), 1, 627 $ A( I2+1, I2 ), 1 ) 628* > Swap A(I1, I1) with A(I2, I2) 629 PIV = A( I1, I1 ) 630 A( I1, I1 ) = A( I2, I2 ) 631 A( I2, I2 ) = PIV 632* > Apply pivots to previous columns of L 633 IF( J.GT.0 ) THEN 634 CALL SSWAP( J*NB, A( I1, 1 ), LDA, 635 $ A( I2, 1 ), LDA ) 636 END IF 637 ENDIF 638 END DO 639* 640* Apply pivots to previous columns of L 641* 642c CALL SLASWP( J*NB, A( 1, 1 ), LDA, 643c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 ) 644 END IF 645 END DO 646 END IF 647* 648* Factor the band matrix 649 CALL SGBTRF( N, N, NB, NB, TB, LDTB, IPIV2, INFO ) 650* 651 RETURN 652* 653* End of SSYTRF_AA_2STAGE 654* 655 END 656