1*> \brief \b ZHETRF_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 ZHETRF_AA_2STAGE + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZHETRF_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* COMPLEX*16 A( LDA, * ), TB( * ), WORK( * ) 31* .. 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> ZHETRF_AA_2STAGE computes the factorization of a double hermitian matrix A 39*> using the Aasen's algorithm. The form of the factorization is 40*> 41*> A = U**H*T*U or A = L*T*L**H 42*> 43*> where U (or L) is a product of permutation and unit upper (lower) 44*> triangular matrices, and T is a hermitian 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 COMPLEX*16 array, dimension (LDA,N) 70*> On entry, the hermitian 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 COMPLEX*16 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 COMPLEX*16 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*> \ingroup complex16SYcomputational 156* 157* ===================================================================== 158 SUBROUTINE ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, 159 $ IPIV2, WORK, LWORK, INFO ) 160* 161* -- LAPACK computational routine -- 162* -- LAPACK is a software package provided by Univ. of Tennessee, -- 163* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 164* 165 IMPLICIT NONE 166* 167* .. Scalar Arguments .. 168 CHARACTER UPLO 169 INTEGER N, LDA, LTB, LWORK, INFO 170* .. 171* .. Array Arguments .. 172 INTEGER IPIV( * ), IPIV2( * ) 173 COMPLEX*16 A( LDA, * ), TB( * ), WORK( * ) 174* .. 175* 176* ===================================================================== 177* .. Parameters .. 178 COMPLEX*16 ZERO, ONE 179 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ), 180 $ ONE = ( 1.0E+0, 0.0E+0 ) ) 181* 182* .. Local Scalars .. 183 LOGICAL UPPER, TQUERY, WQUERY 184 INTEGER I, J, K, I1, I2, TD 185 INTEGER LDTB, NB, KB, JB, NT, IINFO 186 COMPLEX*16 PIV 187* .. 188* .. External Functions .. 189 LOGICAL LSAME 190 INTEGER ILAENV 191 EXTERNAL LSAME, ILAENV 192* .. 193* .. External Subroutines .. 194 EXTERNAL XERBLA, ZCOPY, ZLACGV, ZLACPY, 195 $ ZLASET, ZGBTRF, ZGEMM, ZGETRF, 196 $ ZHEGST, ZSWAP, ZTRSM 197* .. 198* .. Intrinsic Functions .. 199 INTRINSIC DCONJG, MIN, MAX 200* .. 201* .. Executable Statements .. 202* 203* Test the input parameters. 204* 205 INFO = 0 206 UPPER = LSAME( UPLO, 'U' ) 207 WQUERY = ( LWORK.EQ.-1 ) 208 TQUERY = ( LTB.EQ.-1 ) 209 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 210 INFO = -1 211 ELSE IF( N.LT.0 ) THEN 212 INFO = -2 213 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 214 INFO = -4 215 ELSE IF ( LTB .LT. 4*N .AND. .NOT.TQUERY ) THEN 216 INFO = -6 217 ELSE IF ( LWORK .LT. N .AND. .NOT.WQUERY ) THEN 218 INFO = -10 219 END IF 220* 221 IF( INFO.NE.0 ) THEN 222 CALL XERBLA( 'ZHETRF_AA_2STAGE', -INFO ) 223 RETURN 224 END IF 225* 226* Answer the query 227* 228 NB = ILAENV( 1, 'ZHETRF_AA_2STAGE', UPLO, N, -1, -1, -1 ) 229 IF( INFO.EQ.0 ) THEN 230 IF( TQUERY ) THEN 231 TB( 1 ) = (3*NB+1)*N 232 END IF 233 IF( WQUERY ) THEN 234 WORK( 1 ) = N*NB 235 END IF 236 END IF 237 IF( TQUERY .OR. WQUERY ) THEN 238 RETURN 239 END IF 240* 241* Quick return 242* 243 IF ( N.EQ.0 ) THEN 244 RETURN 245 ENDIF 246* 247* Determine the number of the block size 248* 249 LDTB = LTB/N 250 IF( LDTB .LT. 3*NB+1 ) THEN 251 NB = (LDTB-1)/3 252 END IF 253 IF( LWORK .LT. NB*N ) THEN 254 NB = LWORK/N 255 END IF 256* 257* Determine the number of the block columns 258* 259 NT = (N+NB-1)/NB 260 TD = 2*NB 261 KB = MIN(NB, N) 262* 263* Initialize vectors/matrices 264* 265 DO J = 1, KB 266 IPIV( J ) = J 267 END DO 268* 269* Save NB 270* 271 TB( 1 ) = NB 272* 273 IF( UPPER ) THEN 274* 275* ..................................................... 276* Factorize A as U**H*D*U using the upper triangle of A 277* ..................................................... 278* 279 DO J = 0, NT-1 280* 281* Generate Jth column of W and H 282* 283 KB = MIN(NB, N-J*NB) 284 DO I = 1, J-1 285 IF( I.EQ.1 ) THEN 286* H(I,J) = T(I,I)*U(I,J) + T(I+1,I)*U(I+1,J) 287 IF( I .EQ. (J-1) ) THEN 288 JB = NB+KB 289 ELSE 290 JB = 2*NB 291 END IF 292 CALL ZGEMM( 'NoTranspose', 'NoTranspose', 293 $ NB, KB, JB, 294 $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1, 295 $ A( (I-1)*NB+1, J*NB+1 ), LDA, 296 $ ZERO, WORK( I*NB+1 ), N ) 297 ELSE 298* 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) 299 IF( I .EQ. (J-1) ) THEN 300 JB = 2*NB+KB 301 ELSE 302 JB = 3*NB 303 END IF 304 CALL ZGEMM( 'NoTranspose', 'NoTranspose', 305 $ NB, KB, JB, 306 $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ), 307 $ LDTB-1, 308 $ A( (I-2)*NB+1, J*NB+1 ), LDA, 309 $ ZERO, WORK( I*NB+1 ), N ) 310 END IF 311 END DO 312* 313* Compute T(J,J) 314* 315 CALL ZLACPY( 'Upper', KB, KB, A( J*NB+1, J*NB+1 ), LDA, 316 $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 317 IF( J.GT.1 ) THEN 318* T(J,J) = U(1:J,J)'*H(1:J) 319 CALL ZGEMM( 'Conjugate transpose', 'NoTranspose', 320 $ KB, KB, (J-1)*NB, 321 $ -ONE, A( 1, J*NB+1 ), LDA, 322 $ WORK( NB+1 ), N, 323 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 324* T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J) 325 CALL ZGEMM( 'Conjugate transpose', 'NoTranspose', 326 $ KB, NB, KB, 327 $ ONE, A( (J-1)*NB+1, J*NB+1 ), LDA, 328 $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1, 329 $ ZERO, WORK( 1 ), N ) 330 CALL ZGEMM( 'NoTranspose', 'NoTranspose', 331 $ KB, KB, NB, 332 $ -ONE, WORK( 1 ), N, 333 $ A( (J-2)*NB+1, J*NB+1 ), LDA, 334 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 335 END IF 336 IF( J.GT.0 ) THEN 337 CALL ZHEGST( 1, 'Upper', KB, 338 $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1, 339 $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO ) 340 END IF 341* 342* Expand T(J,J) into full format 343* 344 DO I = 1, KB 345 TB( TD+1 + (J*NB+I-1)*LDTB ) 346 $ = REAL( TB( TD+1 + (J*NB+I-1)*LDTB ) ) 347 DO K = I+1, KB 348 TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB ) 349 $ = DCONJG( 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 ZGEMM( '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 ZGEMM( '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 ZGEMM( 'Conjugate 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 ZGETRF 383* 384 DO K = 1, NB 385 CALL ZCOPY( 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 ZGETRF( 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* 403* Copy only L-factor 404* 405 CALL ZCOPY( N-K-(J+1)*NB, 406 $ WORK( K+1+(K-1)*N ), 1, 407 $ A( J*NB+K, (J+1)*NB+K+1 ), LDA ) 408* 409* Transpose U-factor to be copied back into T(J+1, J) 410* 411 CALL ZLACGV( K, WORK( 1+(K-1)*N ), 1 ) 412 END DO 413* 414* Compute T(J+1, J), zero out for GEMM update 415* 416 KB = MIN(NB, N-(J+1)*NB) 417 CALL ZLASET( 'Full', KB, NB, ZERO, ZERO, 418 $ TB( TD+NB+1 + (J*NB)*LDTB) , LDTB-1 ) 419 CALL ZLACPY( 'Upper', KB, NB, 420 $ WORK, N, 421 $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 ) 422 IF( J.GT.0 ) THEN 423 CALL ZTRSM( 'R', 'U', 'N', 'U', KB, NB, ONE, 424 $ A( (J-1)*NB+1, J*NB+1 ), LDA, 425 $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 ) 426 END IF 427* 428* Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM 429* updates 430* 431 DO K = 1, NB 432 DO I = 1, KB 433 TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB ) 434 $ = DCONJG( TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB ) ) 435 END DO 436 END DO 437 CALL ZLASET( 'Lower', KB, NB, ZERO, ONE, 438 $ A( J*NB+1, (J+1)*NB+1), LDA ) 439* 440* Apply pivots to trailing submatrix of A 441* 442 DO K = 1, KB 443* > Adjust ipiv 444 IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB 445* 446 I1 = (J+1)*NB+K 447 I2 = IPIV( (J+1)*NB+K ) 448 IF( I1.NE.I2 ) THEN 449* > Apply pivots to previous columns of L 450 CALL ZSWAP( K-1, A( (J+1)*NB+1, I1 ), 1, 451 $ A( (J+1)*NB+1, I2 ), 1 ) 452* > Swap A(I1+1:M, I1) with A(I2, I1+1:M) 453 IF( I2.GT.(I1+1) ) THEN 454 CALL ZSWAP( I2-I1-1, A( I1, I1+1 ), LDA, 455 $ A( I1+1, I2 ), 1 ) 456 CALL ZLACGV( I2-I1-1, A( I1+1, I2 ), 1 ) 457 END IF 458 CALL ZLACGV( I2-I1, A( I1, I1+1 ), LDA ) 459* > Swap A(I2+1:M, I1) with A(I2+1:M, I2) 460 IF( I2.LT.N ) 461 $ CALL ZSWAP( N-I2, A( I1, I2+1 ), LDA, 462 $ A( I2, I2+1 ), LDA ) 463* > Swap A(I1, I1) with A(I2, I2) 464 PIV = A( I1, I1 ) 465 A( I1, I1 ) = A( I2, I2 ) 466 A( I2, I2 ) = PIV 467* > Apply pivots to previous columns of L 468 IF( J.GT.0 ) THEN 469 CALL ZSWAP( J*NB, A( 1, I1 ), 1, 470 $ A( 1, I2 ), 1 ) 471 END IF 472 ENDIF 473 END DO 474 END IF 475 END DO 476 ELSE 477* 478* ..................................................... 479* Factorize A as L*D*L**H using the lower triangle of A 480* ..................................................... 481* 482 DO J = 0, NT-1 483* 484* Generate Jth column of W and H 485* 486 KB = MIN(NB, N-J*NB) 487 DO I = 1, J-1 488 IF( I.EQ.1 ) THEN 489* H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)' 490 IF( I .EQ. (J-1) ) THEN 491 JB = NB+KB 492 ELSE 493 JB = 2*NB 494 END IF 495 CALL ZGEMM( 'NoTranspose', 'Conjugate transpose', 496 $ NB, KB, JB, 497 $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1, 498 $ A( J*NB+1, (I-1)*NB+1 ), LDA, 499 $ ZERO, WORK( I*NB+1 ), N ) 500 ELSE 501* 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)' 502 IF( I .EQ. (J-1) ) THEN 503 JB = 2*NB+KB 504 ELSE 505 JB = 3*NB 506 END IF 507 CALL ZGEMM( 'NoTranspose', 'Conjugate transpose', 508 $ NB, KB, JB, 509 $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ), 510 $ LDTB-1, 511 $ A( J*NB+1, (I-2)*NB+1 ), LDA, 512 $ ZERO, WORK( I*NB+1 ), N ) 513 END IF 514 END DO 515* 516* Compute T(J,J) 517* 518 CALL ZLACPY( 'Lower', KB, KB, A( J*NB+1, J*NB+1 ), LDA, 519 $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 520 IF( J.GT.1 ) THEN 521* T(J,J) = L(J,1:J)*H(1:J) 522 CALL ZGEMM( 'NoTranspose', 'NoTranspose', 523 $ KB, KB, (J-1)*NB, 524 $ -ONE, A( J*NB+1, 1 ), LDA, 525 $ WORK( NB+1 ), N, 526 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 527* T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)' 528 CALL ZGEMM( 'NoTranspose', 'NoTranspose', 529 $ KB, NB, KB, 530 $ ONE, A( J*NB+1, (J-1)*NB+1 ), LDA, 531 $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1, 532 $ ZERO, WORK( 1 ), N ) 533 CALL ZGEMM( 'NoTranspose', 'Conjugate transpose', 534 $ KB, KB, NB, 535 $ -ONE, WORK( 1 ), N, 536 $ A( J*NB+1, (J-2)*NB+1 ), LDA, 537 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 ) 538 END IF 539 IF( J.GT.0 ) THEN 540 CALL ZHEGST( 1, 'Lower', KB, 541 $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1, 542 $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO ) 543 END IF 544* 545* Expand T(J,J) into full format 546* 547 DO I = 1, KB 548 TB( TD+1 + (J*NB+I-1)*LDTB ) 549 $ = REAL( TB( TD+1 + (J*NB+I-1)*LDTB ) ) 550 DO K = I+1, KB 551 TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB ) 552 $ = DCONJG( TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB ) ) 553 END DO 554 END DO 555* 556 IF( J.LT.NT-1 ) THEN 557 IF( J.GT.0 ) THEN 558* 559* Compute H(J,J) 560* 561 IF( J.EQ.1 ) THEN 562 CALL ZGEMM( 'NoTranspose', 'Conjugate transpose', 563 $ KB, KB, KB, 564 $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1, 565 $ A( J*NB+1, (J-1)*NB+1 ), LDA, 566 $ ZERO, WORK( J*NB+1 ), N ) 567 ELSE 568 CALL ZGEMM( 'NoTranspose', 'Conjugate transpose', 569 $ KB, KB, NB+KB, 570 $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ), 571 $ LDTB-1, 572 $ A( J*NB+1, (J-2)*NB+1 ), LDA, 573 $ ZERO, WORK( J*NB+1 ), N ) 574 END IF 575* 576* Update with the previous column 577* 578 CALL ZGEMM( 'NoTranspose', 'NoTranspose', 579 $ N-(J+1)*NB, NB, J*NB, 580 $ -ONE, A( (J+1)*NB+1, 1 ), LDA, 581 $ WORK( NB+1 ), N, 582 $ ONE, A( (J+1)*NB+1, J*NB+1 ), LDA ) 583 END IF 584* 585* Factorize panel 586* 587 CALL ZGETRF( N-(J+1)*NB, NB, 588 $ A( (J+1)*NB+1, J*NB+1 ), LDA, 589 $ IPIV( (J+1)*NB+1 ), IINFO ) 590c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN 591c INFO = IINFO+(J+1)*NB 592c END IF 593* 594* Compute T(J+1, J), zero out for GEMM update 595* 596 KB = MIN(NB, N-(J+1)*NB) 597 CALL ZLASET( 'Full', KB, NB, ZERO, ZERO, 598 $ TB( TD+NB+1 + (J*NB)*LDTB) , LDTB-1 ) 599 CALL ZLACPY( 'Upper', KB, NB, 600 $ A( (J+1)*NB+1, J*NB+1 ), LDA, 601 $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 ) 602 IF( J.GT.0 ) THEN 603 CALL ZTRSM( 'R', 'L', 'C', 'U', KB, NB, ONE, 604 $ A( J*NB+1, (J-1)*NB+1 ), LDA, 605 $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 ) 606 END IF 607* 608* Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM 609* updates 610* 611 DO K = 1, NB 612 DO I = 1, KB 613 TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB ) 614 $ = DCONJG( TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB ) ) 615 END DO 616 END DO 617 CALL ZLASET( 'Upper', KB, NB, ZERO, ONE, 618 $ A( (J+1)*NB+1, J*NB+1), LDA ) 619* 620* Apply pivots to trailing submatrix of A 621* 622 DO K = 1, KB 623* > Adjust ipiv 624 IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB 625* 626 I1 = (J+1)*NB+K 627 I2 = IPIV( (J+1)*NB+K ) 628 IF( I1.NE.I2 ) THEN 629* > Apply pivots to previous columns of L 630 CALL ZSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA, 631 $ A( I2, (J+1)*NB+1 ), LDA ) 632* > Swap A(I1+1:M, I1) with A(I2, I1+1:M) 633 IF( I2.GT.(I1+1) ) THEN 634 CALL ZSWAP( I2-I1-1, A( I1+1, I1 ), 1, 635 $ A( I2, I1+1 ), LDA ) 636 CALL ZLACGV( I2-I1-1, A( I2, I1+1 ), LDA ) 637 END IF 638 CALL ZLACGV( I2-I1, A( I1+1, I1 ), 1 ) 639* > Swap A(I2+1:M, I1) with A(I2+1:M, I2) 640 IF( I2.LT.N ) 641 $ CALL ZSWAP( N-I2, A( I2+1, I1 ), 1, 642 $ A( I2+1, I2 ), 1 ) 643* > Swap A(I1, I1) with A(I2, I2) 644 PIV = A( I1, I1 ) 645 A( I1, I1 ) = A( I2, I2 ) 646 A( I2, I2 ) = PIV 647* > Apply pivots to previous columns of L 648 IF( J.GT.0 ) THEN 649 CALL ZSWAP( J*NB, A( I1, 1 ), LDA, 650 $ A( I2, 1 ), LDA ) 651 END IF 652 ENDIF 653 END DO 654* 655* Apply pivots to previous columns of L 656* 657c CALL ZLASWP( J*NB, A( 1, 1 ), LDA, 658c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 ) 659 END IF 660 END DO 661 END IF 662* 663* Factor the band matrix 664 CALL ZGBTRF( N, N, NB, NB, TB, LDTB, IPIV2, INFO ) 665* 666 RETURN 667* 668* End of ZHETRF_AA_2STAGE 669* 670 END 671