243 CHARACTER COMPQ, COMPZ
244 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
247 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
248 $ Z( LDZ, * ), WORK( * )
255 parameter( cone = ( 1.0e+0, 0.0e+0 ),
256 $ czero = ( 0.0e+0, 0.0e+0 ) )
259 LOGICAL BLK22, INITQ, INITZ, LQUERY, WANTQ, WANTZ
260 CHARACTER*1 COMPQ2, COMPZ2
261 INTEGER COLA, I, IERR, J, J0, JCOL, JJ, JROW, K,
262 $ KACC22, LEN, LWKOPT, N2NB, NB, NBLST, NBMIN,
263 $ NH, NNB, NX, PPW, PPWO, PW, TOP, TOPQ
265 COMPLEX C1, C2, CTEMP, S, S1, S2, TEMP, TEMP1, TEMP2,
278 INTRINSIC real, cmplx, conjg, max
285 nb =
ilaenv( 1,
'CGGHD3',
' ', n, ilo, ihi, -1 )
286 lwkopt = max( 6*n*nb, 1 )
287 work( 1 ) = cmplx( lwkopt )
288 initq =
lsame( compq,
'I' )
289 wantq = initq .OR.
lsame( compq,
'V' )
290 initz =
lsame( compz,
'I' )
291 wantz = initz .OR.
lsame( compz,
'V' )
292 lquery = ( lwork.EQ.-1 )
294 IF( .NOT.
lsame( compq,
'N' ) .AND. .NOT.wantq )
THEN
296 ELSE IF( .NOT.
lsame( compz,
'N' ) .AND. .NOT.wantz )
THEN
298 ELSE IF( n.LT.0 )
THEN
300 ELSE IF( ilo.LT.1 )
THEN
302 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
304 ELSE IF( lda.LT.max( 1, n ) )
THEN
306 ELSE IF( ldb.LT.max( 1, n ) )
THEN
308 ELSE IF( ( wantq .AND. ldq.LT.n ) .OR. ldq.LT.1 )
THEN
310 ELSE IF( ( wantz .AND. ldz.LT.n ) .OR. ldz.LT.1 )
THEN
312 ELSE IF( lwork.LT.1 .AND. .NOT.lquery )
THEN
316 CALL xerbla(
'CGGHD3', -info )
318 ELSE IF( lquery )
THEN
325 $
CALL claset(
'All', n, n, czero, cone, q, ldq )
327 $
CALL claset(
'All', n, n, czero, cone, z, ldz )
332 $
CALL claset(
'Lower', n-1, n-1, czero, czero, b(2, 1), ldb )
344 nbmin =
ilaenv( 2,
'CGGHD3',
' ', n, ilo, ihi, -1 )
345 IF( nb.GT.1 .AND. nb.LT.nh )
THEN
349 nx = max( nb,
ilaenv( 3,
'CGGHD3',
' ', n, ilo, ihi, -1 ) )
354 IF( lwork.LT.lwkopt )
THEN
360 nbmin = max( 2,
ilaenv( 2,
'CGGHD3',
' ', n, ilo, ihi,
362 IF( lwork.GE.6*n*nbmin )
THEN
371 IF( nb.LT.nbmin .OR. nb.GE.nh )
THEN
381 kacc22 =
ilaenv( 16,
'CGGHD3',
' ', n, ilo, ihi, -1 )
383 DO jcol = ilo, ihi-2, nb
384 nnb = min( nb, ihi-jcol-1 )
392 n2nb = ( ihi-jcol-1 ) / nnb - 1
393 nblst = ihi - jcol - n2nb*nnb
394 CALL claset(
'All', nblst, nblst, czero, cone, work, nblst )
395 pw = nblst * nblst + 1
397 CALL claset(
'All', 2*nnb, 2*nnb, czero, cone,
398 $ work( pw ), 2*nnb )
404 DO j = jcol, jcol+nnb-1
411 CALL clartg( temp, a( i, j ), c, s, a( i-1, j ) )
412 a( i, j ) = cmplx( c )
418 ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
420 jrow = j + n2nb*nnb + 2
424 DO jj = ppw, ppw+len-1
425 temp = work( jj + nblst )
426 work( jj + nblst ) = ctemp*temp - s*work( jj )
427 work( jj ) = conjg( s )*temp + ctemp*work( jj )
430 ppw = ppw - nblst - 1
433 ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
435 DO jrow = j0, j+2, -nnb
438 DO i = jrow+nnb-1, jrow, -1
441 DO jj = ppw, ppw+len-1
442 temp = work( jj + 2*nnb )
443 work( jj + 2*nnb ) = ctemp*temp - s*work( jj )
444 work( jj ) = conjg( s )*temp + ctemp*work( jj )
447 ppw = ppw - 2*nnb - 1
449 ppwo = ppwo + 4*nnb*nnb
468 DO i = min( jj+1, ihi ), j+2, -1
472 b( i, jj ) = ctemp*temp - conjg( s )*b( i-1, jj )
473 b( i-1, jj ) = s*temp + ctemp*b( i-1, jj )
479 temp = b( jj+1, jj+1 )
480 CALL clartg( temp, b( jj+1, jj ), c, s,
482 b( jj+1, jj ) = czero
483 CALL crot( jj-top, b( top+1, jj+1 ), 1,
484 $ b( top+1, jj ), 1, c, s )
485 a( jj+1, j ) = cmplx( c )
486 b( jj+1, j ) = -conjg( s )
492 jj = mod( ihi-j-1, 3 )
493 DO i = ihi-j-3, jj+1, -3
494 ctemp = a( j+1+i, j )
503 temp1 = a( k, j+i+1 )
504 temp2 = a( k, j+i+2 )
505 temp3 = a( k, j+i+3 )
506 a( k, j+i+3 ) = c2*temp3 + conjg( s2 )*temp2
507 temp2 = -s2*temp3 + c2*temp2
508 a( k, j+i+2 ) = c1*temp2 + conjg( s1 )*temp1
509 temp1 = -s1*temp2 + c1*temp1
510 a( k, j+i+1 ) = ctemp*temp1 + conjg( s )*temp
511 a( k, j+i ) = -s*temp1 + ctemp*temp
517 c = dble( a( j+1+i, j ) )
518 CALL crot( ihi-top, a( top+1, j+i+1 ), 1,
519 $ a( top+1, j+i ), 1, c,
520 $ -conjg( b( j+1+i, j ) ) )
526 IF ( j .LT. jcol + nnb - 1 )
THEN
539 jrow = ihi - nblst + 1
540 CALL cgemv(
'Conjugate', nblst, len, cone, work,
541 $ nblst, a( jrow, j+1 ), 1, czero,
544 DO i = jrow, jrow+nblst-len-1
545 work( ppw ) = a( i, j+1 )
548 CALL ctrmv(
'Lower',
'Conjugate',
'Non-unit',
549 $ nblst-len, work( len*nblst + 1 ), nblst,
550 $ work( pw+len ), 1 )
551 CALL cgemv(
'Conjugate', len, nblst-len, cone,
552 $ work( (len+1)*nblst - len + 1 ), nblst,
553 $ a( jrow+nblst-len, j+1 ), 1, cone,
554 $ work( pw+len ), 1 )
556 DO i = jrow, jrow+nblst-1
557 a( i, j+1 ) = work( ppw )
574 ppwo = 1 + nblst*nblst
576 DO jrow = j0, jcol+1, -nnb
578 DO i = jrow, jrow+nnb-1
579 work( ppw ) = a( i, j+1 )
583 DO i = jrow+nnb, jrow+nnb+len-1
584 work( ppw ) = a( i, j+1 )
587 CALL ctrmv(
'Upper',
'Conjugate',
'Non-unit', len,
588 $ work( ppwo + nnb ), 2*nnb, work( pw ),
590 CALL ctrmv(
'Lower',
'Conjugate',
'Non-unit', nnb,
591 $ work( ppwo + 2*len*nnb ),
592 $ 2*nnb, work( pw + len ), 1 )
593 CALL cgemv(
'Conjugate', nnb, len, cone,
594 $ work( ppwo ), 2*nnb, a( jrow, j+1 ), 1,
595 $ cone, work( pw ), 1 )
596 CALL cgemv(
'Conjugate', len, nnb, cone,
597 $ work( ppwo + 2*len*nnb + nnb ), 2*nnb,
598 $ a( jrow+nnb, j+1 ), 1, cone,
599 $ work( pw+len ), 1 )
601 DO i = jrow, jrow+len+nnb-1
602 a( i, j+1 ) = work( ppw )
605 ppwo = ppwo + 4*nnb*nnb
612 cola = n - jcol - nnb + 1
614 CALL cgemm(
'Conjugate',
'No Transpose', nblst,
615 $ cola, nblst, cone, work, nblst,
616 $ a( j, jcol+nnb ), lda, czero, work( pw ),
618 CALL clacpy(
'All', nblst, cola, work( pw ), nblst,
619 $ a( j, jcol+nnb ), lda )
620 ppwo = nblst*nblst + 1
622 DO j = j0, jcol+1, -nnb
634 CALL cunm22(
'Left',
'Conjugate', 2*nnb, cola, nnb,
635 $ nnb, work( ppwo ), 2*nnb,
636 $ a( j, jcol+nnb ), lda, work( pw ),
642 CALL cgemm(
'Conjugate',
'No Transpose', 2*nnb,
643 $ cola, 2*nnb, cone, work( ppwo ), 2*nnb,
644 $ a( j, jcol+nnb ), lda, czero, work( pw ),
646 CALL clacpy(
'All', 2*nnb, cola, work( pw ), 2*nnb,
647 $ a( j, jcol+nnb ), lda )
649 ppwo = ppwo + 4*nnb*nnb
657 topq = max( 2, j - jcol + 1 )
663 CALL cgemm(
'No Transpose',
'No Transpose', nh,
664 $ nblst, nblst, cone, q( topq, j ), ldq,
665 $ work, nblst, czero, work( pw ), nh )
666 CALL clacpy(
'All', nh, nblst, work( pw ), nh,
667 $ q( topq, j ), ldq )
668 ppwo = nblst*nblst + 1
670 DO j = j0, jcol+1, -nnb
672 topq = max( 2, j - jcol + 1 )
679 CALL cunm22(
'Right',
'No Transpose', nh, 2*nnb,
680 $ nnb, nnb, work( ppwo ), 2*nnb,
681 $ q( topq, j ), ldq, work( pw ),
687 CALL cgemm(
'No Transpose',
'No Transpose', nh,
688 $ 2*nnb, 2*nnb, cone, q( topq, j ), ldq,
689 $ work( ppwo ), 2*nnb, czero, work( pw ),
691 CALL clacpy(
'All', nh, 2*nnb, work( pw ), nh,
692 $ q( topq, j ), ldq )
694 ppwo = ppwo + 4*nnb*nnb
700 IF ( wantz .OR. top.GT.0 )
THEN
705 CALL claset(
'All', nblst, nblst, czero, cone, work,
707 pw = nblst * nblst + 1
709 CALL claset(
'All', 2*nnb, 2*nnb, czero, cone,
710 $ work( pw ), 2*nnb )
716 DO j = jcol, jcol+nnb-1
717 ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
719 jrow = j + n2nb*nnb + 2
725 DO jj = ppw, ppw+len-1
726 temp = work( jj + nblst )
727 work( jj + nblst ) = ctemp*temp -
728 $ conjg( s )*work( jj )
729 work( jj ) = s*temp + ctemp*work( jj )
732 ppw = ppw - nblst - 1
735 ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
737 DO jrow = j0, j+2, -nnb
740 DO i = jrow+nnb-1, jrow, -1
745 DO jj = ppw, ppw+len-1
746 temp = work( jj + 2*nnb )
747 work( jj + 2*nnb ) = ctemp*temp -
748 $ conjg( s )*work( jj )
749 work( jj ) = s*temp + ctemp*work( jj )
752 ppw = ppw - 2*nnb - 1
754 ppwo = ppwo + 4*nnb*nnb
759 CALL claset(
'Lower', ihi - jcol - 1, nnb, czero, czero,
760 $ a( jcol + 2, jcol ), lda )
761 CALL claset(
'Lower', ihi - jcol - 1, nnb, czero, czero,
762 $ b( jcol + 2, jcol ), ldb )
769 CALL cgemm(
'No Transpose',
'No Transpose', top,
770 $ nblst, nblst, cone, a( 1, j ), lda,
771 $ work, nblst, czero, work( pw ), top )
772 CALL clacpy(
'All', top, nblst, work( pw ), top,
774 ppwo = nblst*nblst + 1
776 DO j = j0, jcol+1, -nnb
781 CALL cunm22(
'Right',
'No Transpose', top, 2*nnb,
782 $ nnb, nnb, work( ppwo ), 2*nnb,
783 $ a( 1, j ), lda, work( pw ),
789 CALL cgemm(
'No Transpose',
'No Transpose', top,
790 $ 2*nnb, 2*nnb, cone, a( 1, j ), lda,
791 $ work( ppwo ), 2*nnb, czero,
793 CALL clacpy(
'All', top, 2*nnb, work( pw ), top,
796 ppwo = ppwo + 4*nnb*nnb
800 CALL cgemm(
'No Transpose',
'No Transpose', top,
801 $ nblst, nblst, cone, b( 1, j ), ldb,
802 $ work, nblst, czero, work( pw ), top )
803 CALL clacpy(
'All', top, nblst, work( pw ), top,
805 ppwo = nblst*nblst + 1
807 DO j = j0, jcol+1, -nnb
812 CALL cunm22(
'Right',
'No Transpose', top, 2*nnb,
813 $ nnb, nnb, work( ppwo ), 2*nnb,
814 $ b( 1, j ), ldb, work( pw ),
820 CALL cgemm(
'No Transpose',
'No Transpose', top,
821 $ 2*nnb, 2*nnb, cone, b( 1, j ), ldb,
822 $ work( ppwo ), 2*nnb, czero,
824 CALL clacpy(
'All', top, 2*nnb, work( pw ), top,
827 ppwo = ppwo + 4*nnb*nnb
836 topq = max( 2, j - jcol + 1 )
842 CALL cgemm(
'No Transpose',
'No Transpose', nh,
843 $ nblst, nblst, cone, z( topq, j ), ldz,
844 $ work, nblst, czero, work( pw ), nh )
845 CALL clacpy(
'All', nh, nblst, work( pw ), nh,
846 $ z( topq, j ), ldz )
847 ppwo = nblst*nblst + 1
849 DO j = j0, jcol+1, -nnb
851 topq = max( 2, j - jcol + 1 )
858 CALL cunm22(
'Right',
'No Transpose', nh, 2*nnb,
859 $ nnb, nnb, work( ppwo ), 2*nnb,
860 $ z( topq, j ), ldz, work( pw ),
866 CALL cgemm(
'No Transpose',
'No Transpose', nh,
867 $ 2*nnb, 2*nnb, cone, z( topq, j ), ldz,
868 $ work( ppwo ), 2*nnb, czero, work( pw ),
870 CALL clacpy(
'All', nh, 2*nnb, work( pw ), nh,
871 $ z( topq, j ), ldz )
873 ppwo = ppwo + 4*nnb*nnb
884 IF ( jcol.NE.ilo )
THEN
892 $
CALL cgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
893 $ ldq, z, ldz, ierr )
894 work( 1 ) = cmplx( lwkopt )