172 INTEGER N, LDA, LTB, LWORK, INFO
175 INTEGER IPIV( * ), IPIV2( * )
176 COMPLEX*16 A( LDA, * ), TB( * ), WORK( * )
182 parameter( zero = ( 0.0e+0, 0.0e+0 ),
183 $ one = ( 1.0e+0, 0.0e+0 ) )
186 LOGICAL UPPER, TQUERY, WQUERY
187 INTEGER I, J, K, I1, I2, TD
188 INTEGER LDTB, NB, KB, JB, NT, IINFO
202 INTRINSIC dconjg, min, max
209 upper =
lsame( uplo,
'U' )
210 wquery = ( lwork.EQ.-1 )
211 tquery = ( ltb.EQ.-1 )
212 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
214 ELSE IF( n.LT.0 )
THEN
216 ELSE IF( lda.LT.max( 1, n ) )
THEN
218 ELSE IF ( ltb .LT. 4*n .AND. .NOT.tquery )
THEN
220 ELSE IF ( lwork .LT. n .AND. .NOT.wquery )
THEN
225 CALL xerbla(
'ZHETRF_AA_2STAGE', -info )
231 nb =
ilaenv( 1,
'ZHETRF_AA_2STAGE', uplo, n, -1, -1, -1 )
240 IF( tquery .OR. wquery )
THEN
253 IF( ldtb .LT. 3*nb+1 )
THEN
256 IF( lwork .LT. nb*n )
THEN
290 IF( i .EQ. (j-1) )
THEN
295 CALL zgemm(
'NoTranspose',
'NoTranspose',
297 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
298 $ a( (i-1)*nb+1, j*nb+1 ), lda,
299 $ zero, work( i*nb+1 ), n )
302 IF( i .EQ. (j-1) )
THEN
307 CALL zgemm(
'NoTranspose',
'NoTranspose',
309 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
311 $ a( (i-2)*nb+1, j*nb+1 ), lda,
312 $ zero, work( i*nb+1 ), n )
318 CALL zlacpy(
'Upper', kb, kb, a( j*nb+1, j*nb+1 ), lda,
319 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
322 CALL zgemm(
'Conjugate transpose',
'NoTranspose',
324 $ -one, a( 1, j*nb+1 ), lda,
326 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
328 CALL zgemm(
'Conjugate transpose',
'NoTranspose',
330 $ one, a( (j-1)*nb+1, j*nb+1 ), lda,
331 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
332 $ zero, work( 1 ), n )
333 CALL zgemm(
'NoTranspose',
'NoTranspose',
335 $ -one, work( 1 ), n,
336 $ a( (j-2)*nb+1, j*nb+1 ), lda,
337 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
340 CALL zhegst( 1,
'Upper', kb,
341 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
342 $ a( (j-1)*nb+1, j*nb+1 ), lda, iinfo )
348 tb( td+1 + (j*nb+i-1)*ldtb )
349 $ = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
351 tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
352 $ = dconjg( tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb ) )
362 CALL zgemm(
'NoTranspose',
'NoTranspose',
364 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
365 $ a( (j-1)*nb+1, j*nb+1 ), lda,
366 $ zero, work( j*nb+1 ), n )
368 CALL zgemm(
'NoTranspose',
'NoTranspose',
370 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
372 $ a( (j-2)*nb+1, j*nb+1 ), lda,
373 $ zero, work( j*nb+1 ), n )
378 CALL zgemm(
'Conjugate transpose',
'NoTranspose',
379 $ nb, n-(j+1)*nb, j*nb,
380 $ -one, work( nb+1 ), n,
381 $ a( 1, (j+1)*nb+1 ), lda,
382 $ one, a( j*nb+1, (j+1)*nb+1 ), lda )
388 CALL zcopy( n-(j+1)*nb,
389 $ a( j*nb+k, (j+1)*nb+1 ), lda,
390 $ work( 1+(k-1)*n ), 1 )
395 CALL zgetrf( n-(j+1)*nb, nb,
397 $ ipiv( (j+1)*nb+1 ), iinfo )
408 CALL zcopy( n-k-(j+1)*nb,
409 $ work( k+1+(k-1)*n ), 1,
410 $ a( j*nb+k, (j+1)*nb+k+1 ), lda )
414 CALL zlacgv( k, work( 1+(k-1)*n ), 1 )
419 kb = min(nb, n-(j+1)*nb)
420 CALL zlaset(
'Full', kb, nb, zero, zero,
421 $ tb( td+nb+1 + (j*nb)*ldtb) , ldtb-1 )
422 CALL zlacpy(
'Upper', kb, nb,
424 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
426 CALL ztrsm(
'R',
'U',
'N',
'U', kb, nb, one,
427 $ a( (j-1)*nb+1, j*nb+1 ), lda,
428 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
436 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
437 $ = dconjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
440 CALL zlaset(
'Lower', kb, nb, zero, one,
441 $ a( j*nb+1, (j+1)*nb+1), lda )
447 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
450 i2 = ipiv( (j+1)*nb+k )
453 CALL zswap( k-1, a( (j+1)*nb+1, i1 ), 1,
454 $ a( (j+1)*nb+1, i2 ), 1 )
456 IF( i2.GT.(i1+1) )
THEN
457 CALL zswap( i2-i1-1, a( i1, i1+1 ), lda,
459 CALL zlacgv( i2-i1-1, a( i1+1, i2 ), 1 )
461 CALL zlacgv( i2-i1, a( i1, i1+1 ), lda )
464 $
CALL zswap( n-i2, a( i1, i2+1 ), lda,
465 $ a( i2, i2+1 ), lda )
468 a( i1, i1 ) = a( i2, i2 )
472 CALL zswap( j*nb, a( 1, i1 ), 1,
493 IF( i .EQ. (j-1) )
THEN
498 CALL zgemm(
'NoTranspose',
'Conjugate transpose',
500 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
501 $ a( j*nb+1, (i-1)*nb+1 ), lda,
502 $ zero, work( i*nb+1 ), n )
505 IF( i .EQ. (j-1) )
THEN
510 CALL zgemm(
'NoTranspose',
'Conjugate transpose',
512 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
514 $ a( j*nb+1, (i-2)*nb+1 ), lda,
515 $ zero, work( i*nb+1 ), n )
521 CALL zlacpy(
'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
522 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
525 CALL zgemm(
'NoTranspose',
'NoTranspose',
527 $ -one, a( j*nb+1, 1 ), lda,
529 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
531 CALL zgemm(
'NoTranspose',
'NoTranspose',
533 $ one, a( j*nb+1, (j-1)*nb+1 ), lda,
534 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
535 $ zero, work( 1 ), n )
536 CALL zgemm(
'NoTranspose',
'Conjugate transpose',
538 $ -one, work( 1 ), n,
539 $ a( j*nb+1, (j-2)*nb+1 ), lda,
540 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
543 CALL zhegst( 1,
'Lower', kb,
544 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
545 $ a( j*nb+1, (j-1)*nb+1 ), lda, iinfo )
551 tb( td+1 + (j*nb+i-1)*ldtb )
552 $ = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
554 tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
555 $ = dconjg( tb( td+(k-i)+1 + (j*nb+i-1)*ldtb ) )
565 CALL zgemm(
'NoTranspose',
'Conjugate transpose',
567 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
568 $ a( j*nb+1, (j-1)*nb+1 ), lda,
569 $ zero, work( j*nb+1 ), n )
571 CALL zgemm(
'NoTranspose',
'Conjugate transpose',
573 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
575 $ a( j*nb+1, (j-2)*nb+1 ), lda,
576 $ zero, work( j*nb+1 ), n )
581 CALL zgemm(
'NoTranspose',
'NoTranspose',
582 $ n-(j+1)*nb, nb, j*nb,
583 $ -one, a( (j+1)*nb+1, 1 ), lda,
585 $ one, a( (j+1)*nb+1, j*nb+1 ), lda )
590 CALL zgetrf( n-(j+1)*nb, nb,
591 $ a( (j+1)*nb+1, j*nb+1 ), lda,
592 $ ipiv( (j+1)*nb+1 ), iinfo )
599 kb = min(nb, n-(j+1)*nb)
600 CALL zlaset(
'Full', kb, nb, zero, zero,
601 $ tb( td+nb+1 + (j*nb)*ldtb) , ldtb-1 )
602 CALL zlacpy(
'Upper', kb, nb,
603 $ a( (j+1)*nb+1, j*nb+1 ), lda,
604 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
606 CALL ztrsm(
'R',
'L',
'C',
'U', kb, nb, one,
607 $ a( j*nb+1, (j-1)*nb+1 ), lda,
608 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
616 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
617 $ = dconjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
620 CALL zlaset(
'Upper', kb, nb, zero, one,
621 $ a( (j+1)*nb+1, j*nb+1), lda )
627 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
630 i2 = ipiv( (j+1)*nb+k )
633 CALL zswap( k-1, a( i1, (j+1)*nb+1 ), lda,
634 $ a( i2, (j+1)*nb+1 ), lda )
636 IF( i2.GT.(i1+1) )
THEN
637 CALL zswap( i2-i1-1, a( i1+1, i1 ), 1,
638 $ a( i2, i1+1 ), lda )
639 CALL zlacgv( i2-i1-1, a( i2, i1+1 ), lda )
641 CALL zlacgv( i2-i1, a( i1+1, i1 ), 1 )
644 $
CALL zswap( n-i2, a( i2+1, i1 ), 1,
648 a( i1, i1 ) = a( i2, i2 )
652 CALL zswap( j*nb, a( i1, 1 ), lda,
667 CALL zgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )