238 CHARACTER compq, compz
239 INTEGER ihi, ilo, info, lda, ldb, ldq, ldz, n, lwork
242 COMPLEX*16 a( lda, * ), b( ldb, * ), q( ldq, * ),
243 $ z( ldz, * ), work( * )
249 COMPLEX*16 cone, czero
250 parameter ( cone = ( 1.0d+0, 0.0d+0 ),
251 $ czero = ( 0.0d+0, 0.0d+0 ) )
254 LOGICAL blk22, initq, initz, lquery, wantq, wantz
255 CHARACTER*1 compq2, compz2
256 INTEGER cola, i, ierr, j, j0, jcol, jj, jrow, k,
257 $ kacc22, len, lwkopt, n2nb, nb, nblst, nbmin,
258 $ nh, nnb, nx, ppw, ppwo, pw, top, topq
260 COMPLEX*16 c1, c2, ctemp, s, s1, s2, temp, temp1, temp2,
272 INTRINSIC dble, dcmplx, dconjg, max
279 nb =
ilaenv( 1,
'ZGGHD3',
' ', n, ilo, ihi, -1 )
280 lwkopt = max( 6*n*nb, 1 )
281 work( 1 ) = dcmplx( lwkopt )
282 initq =
lsame( compq,
'I' )
283 wantq = initq .OR.
lsame( compq,
'V' )
284 initz =
lsame( compz,
'I' )
285 wantz = initz .OR.
lsame( compz,
'V' )
286 lquery = ( lwork.EQ.-1 )
288 IF( .NOT.
lsame( compq,
'N' ) .AND. .NOT.wantq )
THEN
290 ELSE IF( .NOT.
lsame( compz,
'N' ) .AND. .NOT.wantz )
THEN
292 ELSE IF( n.LT.0 )
THEN
294 ELSE IF( ilo.LT.1 )
THEN
296 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
298 ELSE IF( lda.LT.max( 1, n ) )
THEN
300 ELSE IF( ldb.LT.max( 1, n ) )
THEN
302 ELSE IF( ( wantq .AND. ldq.LT.n ) .OR. ldq.LT.1 )
THEN
304 ELSE IF( ( wantz .AND. ldz.LT.n ) .OR. ldz.LT.1 )
THEN
306 ELSE IF( lwork.LT.1 .AND. .NOT.lquery )
THEN
310 CALL xerbla(
'ZGGHD3', -info )
312 ELSE IF( lquery )
THEN
319 $
CALL zlaset(
'All', n, n, czero, cone, q, ldq )
321 $
CALL zlaset(
'All', n, n, czero, cone, z, ldz )
326 $
CALL zlaset(
'Lower', n-1, n-1, czero, czero, b(2, 1), ldb )
338 nbmin =
ilaenv( 2,
'ZGGHD3',
' ', n, ilo, ihi, -1 )
339 IF( nb.GT.1 .AND. nb.LT.nh )
THEN
343 nx = max( nb,
ilaenv( 3,
'ZGGHD3',
' ', n, ilo, ihi, -1 ) )
348 IF( lwork.LT.lwkopt )
THEN
354 nbmin = max( 2,
ilaenv( 2,
'ZGGHD3',
' ', n, ilo, ihi,
356 IF( lwork.GE.6*n*nbmin )
THEN
365 IF( nb.LT.nbmin .OR. nb.GE.nh )
THEN
375 kacc22 =
ilaenv( 16,
'ZGGHD3',
' ', n, ilo, ihi, -1 )
377 DO jcol = ilo, ihi-2, nb
378 nnb = min( nb, ihi-jcol-1 )
386 n2nb = ( ihi-jcol-1 ) / nnb - 1
387 nblst = ihi - jcol - n2nb*nnb
388 CALL zlaset(
'All', nblst, nblst, czero, cone, work, nblst )
389 pw = nblst * nblst + 1
391 CALL zlaset(
'All', 2*nnb, 2*nnb, czero, cone,
392 $ work( pw ), 2*nnb )
398 DO j = jcol, jcol+nnb-1
405 CALL zlartg( temp, a( i, j ), c, s, a( i-1, j ) )
406 a( i, j ) = dcmplx( c )
412 ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
414 jrow = j + n2nb*nnb + 2
418 DO jj = ppw, ppw+len-1
419 temp = work( jj + nblst )
420 work( jj + nblst ) = ctemp*temp - s*work( jj )
421 work( jj ) = dconjg( s )*temp + ctemp*work( jj )
424 ppw = ppw - nblst - 1
427 ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
429 DO jrow = j0, j+2, -nnb
432 DO i = jrow+nnb-1, jrow, -1
435 DO jj = ppw, ppw+len-1
436 temp = work( jj + 2*nnb )
437 work( jj + 2*nnb ) = ctemp*temp - s*work( jj )
438 work( jj ) = dconjg( s )*temp + ctemp*work( jj )
441 ppw = ppw - 2*nnb - 1
443 ppwo = ppwo + 4*nnb*nnb
462 DO i = min( jj+1, ihi ), j+2, -1
466 b( i, jj ) = ctemp*temp - dconjg( s )*b( i-1, jj )
467 b( i-1, jj ) = s*temp + ctemp*b( i-1, jj )
473 temp = b( jj+1, jj+1 )
474 CALL zlartg( temp, b( jj+1, jj ), c, s,
476 b( jj+1, jj ) = czero
477 CALL zrot( jj-top, b( top+1, jj+1 ), 1,
478 $ b( top+1, jj ), 1, c, s )
479 a( jj+1, j ) = dcmplx( c )
480 b( jj+1, j ) = -dconjg( s )
486 jj = mod( ihi-j-1, 3 )
487 DO i = ihi-j-3, jj+1, -3
488 ctemp = a( j+1+i, j )
497 temp1 = a( k, j+i+1 )
498 temp2 = a( k, j+i+2 )
499 temp3 = a( k, j+i+3 )
500 a( k, j+i+3 ) = c2*temp3 + dconjg( s2 )*temp2
501 temp2 = -s2*temp3 + c2*temp2
502 a( k, j+i+2 ) = c1*temp2 + dconjg( s1 )*temp1
503 temp1 = -s1*temp2 + c1*temp1
504 a( k, j+i+1 ) = ctemp*temp1 + dconjg( s )*temp
505 a( k, j+i ) = -s*temp1 + ctemp*temp
511 c = dble( a( j+1+i, j ) )
512 CALL zrot( ihi-top, a( top+1, j+i+1 ), 1,
513 $ a( top+1, j+i ), 1, c,
514 $ -dconjg( b( j+1+i, j ) ) )
520 IF ( j .LT. jcol + nnb - 1 )
THEN
533 jrow = ihi - nblst + 1
534 CALL zgemv(
'Conjugate', nblst, len, cone, work,
535 $ nblst, a( jrow, j+1 ), 1, czero,
538 DO i = jrow, jrow+nblst-len-1
539 work( ppw ) = a( i, j+1 )
542 CALL ztrmv(
'Lower',
'Conjugate',
'Non-unit',
543 $ nblst-len, work( len*nblst + 1 ), nblst,
544 $ work( pw+len ), 1 )
545 CALL zgemv(
'Conjugate', len, nblst-len, cone,
546 $ work( (len+1)*nblst - len + 1 ), nblst,
547 $ a( jrow+nblst-len, j+1 ), 1, cone,
548 $ work( pw+len ), 1 )
550 DO i = jrow, jrow+nblst-1
551 a( i, j+1 ) = work( ppw )
568 ppwo = 1 + nblst*nblst
570 DO jrow = j0, jcol+1, -nnb
572 DO i = jrow, jrow+nnb-1
573 work( ppw ) = a( i, j+1 )
577 DO i = jrow+nnb, jrow+nnb+len-1
578 work( ppw ) = a( i, j+1 )
581 CALL ztrmv(
'Upper',
'Conjugate',
'Non-unit', len,
582 $ work( ppwo + nnb ), 2*nnb, work( pw ),
584 CALL ztrmv(
'Lower',
'Conjugate',
'Non-unit', nnb,
585 $ work( ppwo + 2*len*nnb ),
586 $ 2*nnb, work( pw + len ), 1 )
587 CALL zgemv(
'Conjugate', nnb, len, cone,
588 $ work( ppwo ), 2*nnb, a( jrow, j+1 ), 1,
589 $ cone, work( pw ), 1 )
590 CALL zgemv(
'Conjugate', len, nnb, cone,
591 $ work( ppwo + 2*len*nnb + nnb ), 2*nnb,
592 $ a( jrow+nnb, j+1 ), 1, cone,
593 $ work( pw+len ), 1 )
595 DO i = jrow, jrow+len+nnb-1
596 a( i, j+1 ) = work( ppw )
599 ppwo = ppwo + 4*nnb*nnb
606 cola = n - jcol - nnb + 1
608 CALL zgemm(
'Conjugate',
'No Transpose', nblst,
609 $ cola, nblst, cone, work, nblst,
610 $ a( j, jcol+nnb ), lda, czero, work( pw ),
612 CALL zlacpy(
'All', nblst, cola, work( pw ), nblst,
613 $ a( j, jcol+nnb ), lda )
614 ppwo = nblst*nblst + 1
616 DO j = j0, jcol+1, -nnb
628 CALL zunm22(
'Left',
'Conjugate', 2*nnb, cola, nnb,
629 $ nnb, work( ppwo ), 2*nnb,
630 $ a( j, jcol+nnb ), lda, work( pw ),
636 CALL zgemm(
'Conjugate',
'No Transpose', 2*nnb,
637 $ cola, 2*nnb, cone, work( ppwo ), 2*nnb,
638 $ a( j, jcol+nnb ), lda, czero, work( pw ),
640 CALL zlacpy(
'All', 2*nnb, cola, work( pw ), 2*nnb,
641 $ a( j, jcol+nnb ), lda )
643 ppwo = ppwo + 4*nnb*nnb
651 topq = max( 2, j - jcol + 1 )
657 CALL zgemm(
'No Transpose',
'No Transpose', nh,
658 $ nblst, nblst, cone, q( topq, j ), ldq,
659 $ work, nblst, czero, work( pw ), nh )
660 CALL zlacpy(
'All', nh, nblst, work( pw ), nh,
661 $ q( topq, j ), ldq )
662 ppwo = nblst*nblst + 1
664 DO j = j0, jcol+1, -nnb
666 topq = max( 2, j - jcol + 1 )
673 CALL zunm22(
'Right',
'No Transpose', nh, 2*nnb,
674 $ nnb, nnb, work( ppwo ), 2*nnb,
675 $ q( topq, j ), ldq, work( pw ),
681 CALL zgemm(
'No Transpose',
'No Transpose', nh,
682 $ 2*nnb, 2*nnb, cone, q( topq, j ), ldq,
683 $ work( ppwo ), 2*nnb, czero, work( pw ),
685 CALL zlacpy(
'All', nh, 2*nnb, work( pw ), nh,
686 $ q( topq, j ), ldq )
688 ppwo = ppwo + 4*nnb*nnb
694 IF ( wantz .OR. top.GT.0 )
THEN
699 CALL zlaset(
'All', nblst, nblst, czero, cone, work,
701 pw = nblst * nblst + 1
703 CALL zlaset(
'All', 2*nnb, 2*nnb, czero, cone,
704 $ work( pw ), 2*nnb )
710 DO j = jcol, jcol+nnb-1
711 ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
713 jrow = j + n2nb*nnb + 2
719 DO jj = ppw, ppw+len-1
720 temp = work( jj + nblst )
721 work( jj + nblst ) = ctemp*temp -
722 $ dconjg( s )*work( jj )
723 work( jj ) = s*temp + ctemp*work( jj )
726 ppw = ppw - nblst - 1
729 ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
731 DO jrow = j0, j+2, -nnb
734 DO i = jrow+nnb-1, jrow, -1
739 DO jj = ppw, ppw+len-1
740 temp = work( jj + 2*nnb )
741 work( jj + 2*nnb ) = ctemp*temp -
742 $ dconjg( s )*work( jj )
743 work( jj ) = s*temp + ctemp*work( jj )
746 ppw = ppw - 2*nnb - 1
748 ppwo = ppwo + 4*nnb*nnb
753 CALL zlaset(
'Lower', ihi - jcol - 1, nnb, czero, czero,
754 $ a( jcol + 2, jcol ), lda )
755 CALL zlaset(
'Lower', ihi - jcol - 1, nnb, czero, czero,
756 $ b( jcol + 2, jcol ), ldb )
763 CALL zgemm(
'No Transpose',
'No Transpose', top,
764 $ nblst, nblst, cone, a( 1, j ), lda,
765 $ work, nblst, czero, work( pw ), top )
766 CALL zlacpy(
'All', top, nblst, work( pw ), top,
768 ppwo = nblst*nblst + 1
770 DO j = j0, jcol+1, -nnb
775 CALL zunm22(
'Right',
'No Transpose', top, 2*nnb,
776 $ nnb, nnb, work( ppwo ), 2*nnb,
777 $ a( 1, j ), lda, work( pw ),
783 CALL zgemm(
'No Transpose',
'No Transpose', top,
784 $ 2*nnb, 2*nnb, cone, a( 1, j ), lda,
785 $ work( ppwo ), 2*nnb, czero,
787 CALL zlacpy(
'All', top, 2*nnb, work( pw ), top,
790 ppwo = ppwo + 4*nnb*nnb
794 CALL zgemm(
'No Transpose',
'No Transpose', top,
795 $ nblst, nblst, cone, b( 1, j ), ldb,
796 $ work, nblst, czero, work( pw ), top )
797 CALL zlacpy(
'All', top, nblst, work( pw ), top,
799 ppwo = nblst*nblst + 1
801 DO j = j0, jcol+1, -nnb
806 CALL zunm22(
'Right',
'No Transpose', top, 2*nnb,
807 $ nnb, nnb, work( ppwo ), 2*nnb,
808 $ b( 1, j ), ldb, work( pw ),
814 CALL zgemm(
'No Transpose',
'No Transpose', top,
815 $ 2*nnb, 2*nnb, cone, b( 1, j ), ldb,
816 $ work( ppwo ), 2*nnb, czero,
818 CALL zlacpy(
'All', top, 2*nnb, work( pw ), top,
821 ppwo = ppwo + 4*nnb*nnb
830 topq = max( 2, j - jcol + 1 )
836 CALL zgemm(
'No Transpose',
'No Transpose', nh,
837 $ nblst, nblst, cone, z( topq, j ), ldz,
838 $ work, nblst, czero, work( pw ), nh )
839 CALL zlacpy(
'All', nh, nblst, work( pw ), nh,
840 $ z( topq, j ), ldz )
841 ppwo = nblst*nblst + 1
843 DO j = j0, jcol+1, -nnb
845 topq = max( 2, j - jcol + 1 )
852 CALL zunm22(
'Right',
'No Transpose', nh, 2*nnb,
853 $ nnb, nnb, work( ppwo ), 2*nnb,
854 $ z( topq, j ), ldz, work( pw ),
860 CALL zgemm(
'No Transpose',
'No Transpose', nh,
861 $ 2*nnb, 2*nnb, cone, z( topq, j ), ldz,
862 $ work( ppwo ), 2*nnb, czero, work( pw ),
864 CALL zlacpy(
'All', nh, 2*nnb, work( pw ), nh,
865 $ z( topq, j ), ldz )
867 ppwo = ppwo + 4*nnb*nnb
878 IF ( jcol.NE.ilo )
THEN
886 $
CALL zgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
887 $ ldq, z, ldz, ierr )
888 work( 1 ) = dcmplx( lwkopt )
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
ZGGHRD
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zunm22(SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC, WORK, LWORK, INFO)
ZUNM22 multiplies a general matrix by a banded unitary matrix.
subroutine ztrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRMV
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine zrot(N, CX, INCX, CY, INCY, C, S)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
logical function lsame(CA, CB)
LSAME
subroutine zlartg(F, G, CS, SN, R)
ZLARTG generates a plane rotation with real cosine and complex sine.