172 INTEGER N, LDA, LTB, LWORK, INFO
175 INTEGER IPIV( * ), IPIV2( * )
176 COMPLEX 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
203 INTRINSIC conjg, min, max
210 upper =
lsame( uplo,
'U' )
211 wquery = ( lwork.EQ.-1 )
212 tquery = ( ltb.EQ.-1 )
213 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
215 ELSE IF( n.LT.0 )
THEN
217 ELSE IF( lda.LT.max( 1, n ) )
THEN
219 ELSE IF ( ltb .LT. 4*n .AND. .NOT.tquery )
THEN
221 ELSE IF ( lwork .LT. n .AND. .NOT.wquery )
THEN
226 CALL xerbla(
'CHETRF_AA_2STAGE', -info )
232 nb =
ilaenv( 1,
'CHETRF_AA_2STAGE', uplo, n, -1, -1, -1 )
241 IF( tquery .OR. wquery )
THEN
254 IF( ldtb .LT. 3*nb+1 )
THEN
257 IF( lwork .LT. nb*n )
THEN
291 IF( i .EQ. (j-1) )
THEN
296 CALL cgemm(
'NoTranspose',
'NoTranspose',
298 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
299 $ a( (i-1)*nb+1, j*nb+1 ), lda,
300 $ zero, work( i*nb+1 ), n )
303 IF( i .EQ. (j-1) )
THEN
308 CALL cgemm(
'NoTranspose',
'NoTranspose',
310 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
312 $ a( (i-2)*nb+1, j*nb+1 ), lda,
313 $ zero, work( i*nb+1 ), n )
319 CALL clacpy(
'Upper', kb, kb, a( j*nb+1, j*nb+1 ), lda,
320 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
323 CALL cgemm(
'Conjugate transpose',
'NoTranspose',
325 $ -one, a( 1, j*nb+1 ), lda,
327 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
329 CALL cgemm(
'Conjugate transpose',
'NoTranspose',
331 $ one, a( (j-1)*nb+1, j*nb+1 ), lda,
332 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
333 $ zero, work( 1 ), n )
334 CALL cgemm(
'NoTranspose',
'NoTranspose',
336 $ -one, work( 1 ), n,
337 $ a( (j-2)*nb+1, j*nb+1 ), lda,
338 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
341 CALL chegst( 1,
'Upper', kb,
342 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
343 $ a( (j-1)*nb+1, j*nb+1 ), lda, iinfo )
349 tb( td+1 + (j*nb+i-1)*ldtb )
350 $ = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
352 tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
353 $ = conjg( tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb ) )
363 CALL cgemm(
'NoTranspose',
'NoTranspose',
365 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
366 $ a( (j-1)*nb+1, j*nb+1 ), lda,
367 $ zero, work( j*nb+1 ), n )
369 CALL cgemm(
'NoTranspose',
'NoTranspose',
371 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
373 $ a( (j-2)*nb+1, j*nb+1 ), lda,
374 $ zero, work( j*nb+1 ), n )
379 CALL cgemm(
'Conjugate transpose',
'NoTranspose',
380 $ nb, n-(j+1)*nb, j*nb,
381 $ -one, work( nb+1 ), n,
382 $ a( 1, (j+1)*nb+1 ), lda,
383 $ one, a( j*nb+1, (j+1)*nb+1 ), lda )
389 CALL ccopy( n-(j+1)*nb,
390 $ a( j*nb+k, (j+1)*nb+1 ), lda,
391 $ work( 1+(k-1)*n ), 1 )
396 CALL cgetrf( n-(j+1)*nb, nb,
398 $ ipiv( (j+1)*nb+1 ), iinfo )
409 CALL ccopy( n-k-(j+1)*nb,
410 $ work( k+1+(k-1)*n ), 1,
411 $ a( j*nb+k, (j+1)*nb+k+1 ), lda )
415 CALL clacgv( k, work( 1+(k-1)*n ), 1 )
420 kb = min(nb, n-(j+1)*nb)
421 CALL claset(
'Full', kb, nb, zero, zero,
422 $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
423 CALL clacpy(
'Upper', kb, nb,
425 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
427 CALL ctrsm(
'R',
'U',
'N',
'U', kb, nb, one,
428 $ a( (j-1)*nb+1, j*nb+1 ), lda,
429 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
437 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
438 $ = conjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
441 CALL claset(
'Lower', kb, nb, zero, one,
442 $ a( j*nb+1, (j+1)*nb+1), lda )
448 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
451 i2 = ipiv( (j+1)*nb+k )
454 CALL cswap( k-1, a( (j+1)*nb+1, i1 ), 1,
455 $ a( (j+1)*nb+1, i2 ), 1 )
457 IF( i2.GT.(i1+1) )
THEN
458 CALL cswap( i2-i1-1, a( i1, i1+1 ), lda,
460 CALL clacgv( i2-i1-1, a( i1+1, i2 ), 1 )
462 CALL clacgv( i2-i1, a( i1, i1+1 ), lda )
465 $
CALL cswap( n-i2, a( i1, i2+1 ), lda,
466 $ a( i2, i2+1 ), lda )
469 a( i1, i1 ) = a( i2, i2 )
473 CALL cswap( j*nb, a( 1, i1 ), 1,
494 IF( i .EQ. (j-1) )
THEN
499 CALL cgemm(
'NoTranspose',
'Conjugate transpose',
501 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
502 $ a( j*nb+1, (i-1)*nb+1 ), lda,
503 $ zero, work( i*nb+1 ), n )
506 IF( i .EQ. (j-1) )
THEN
511 CALL cgemm(
'NoTranspose',
'Conjugate transpose',
513 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
515 $ a( j*nb+1, (i-2)*nb+1 ), lda,
516 $ zero, work( i*nb+1 ), n )
522 CALL clacpy(
'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
523 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
526 CALL cgemm(
'NoTranspose',
'NoTranspose',
528 $ -one, a( j*nb+1, 1 ), lda,
530 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
532 CALL cgemm(
'NoTranspose',
'NoTranspose',
534 $ one, a( j*nb+1, (j-1)*nb+1 ), lda,
535 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
536 $ zero, work( 1 ), n )
537 CALL cgemm(
'NoTranspose',
'Conjugate transpose',
539 $ -one, work( 1 ), n,
540 $ a( j*nb+1, (j-2)*nb+1 ), lda,
541 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
544 CALL chegst( 1,
'Lower', kb,
545 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
546 $ a( j*nb+1, (j-1)*nb+1 ), lda, iinfo )
552 tb( td+1 + (j*nb+i-1)*ldtb )
553 $ = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
555 tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
556 $ = conjg( tb( td+(k-i)+1 + (j*nb+i-1)*ldtb ) )
566 CALL cgemm(
'NoTranspose',
'Conjugate transpose',
568 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
569 $ a( j*nb+1, (j-1)*nb+1 ), lda,
570 $ zero, work( j*nb+1 ), n )
572 CALL cgemm(
'NoTranspose',
'Conjugate transpose',
574 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
576 $ a( j*nb+1, (j-2)*nb+1 ), lda,
577 $ zero, work( j*nb+1 ), n )
582 CALL cgemm(
'NoTranspose',
'NoTranspose',
583 $ n-(j+1)*nb, nb, j*nb,
584 $ -one, a( (j+1)*nb+1, 1 ), lda,
586 $ one, a( (j+1)*nb+1, j*nb+1 ), lda )
591 CALL cgetrf( n-(j+1)*nb, nb,
592 $ a( (j+1)*nb+1, j*nb+1 ), lda,
593 $ ipiv( (j+1)*nb+1 ), iinfo )
600 kb = min(nb, n-(j+1)*nb)
601 CALL claset(
'Full', kb, nb, zero, zero,
602 $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
603 CALL clacpy(
'Upper', kb, nb,
604 $ a( (j+1)*nb+1, j*nb+1 ), lda,
605 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
607 CALL ctrsm(
'R',
'L',
'C',
'U', kb, nb, one,
608 $ a( j*nb+1, (j-1)*nb+1 ), lda,
609 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
617 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
618 $ = conjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
621 CALL claset(
'Upper', kb, nb, zero, one,
622 $ a( (j+1)*nb+1, j*nb+1), lda )
628 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
631 i2 = ipiv( (j+1)*nb+k )
634 CALL cswap( k-1, a( i1, (j+1)*nb+1 ), lda,
635 $ a( i2, (j+1)*nb+1 ), lda )
637 IF( i2.GT.(i1+1) )
THEN
638 CALL cswap( i2-i1-1, a( i1+1, i1 ), 1,
639 $ a( i2, i1+1 ), lda )
640 CALL clacgv( i2-i1-1, a( i2, i1+1 ), lda )
642 CALL clacgv( i2-i1, a( i1+1, i1 ), 1 )
645 $
CALL cswap( n-i2, a( i2+1, i1 ), 1,
649 a( i1, i1 ) = a( i2, i2 )
653 CALL cswap( j*nb, a( i1, 1 ), lda,
668 CALL cgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )