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,
277 INTRINSIC REAL, cmplx, conjg, max
284 nb =
ilaenv( 1,
'CGGHD3',
' ', n, ilo, ihi, -1 )
285 lwkopt = max( 6*n*nb, 1 )
286 work( 1 ) = cmplx( lwkopt )
287 initq =
lsame( compq,
'I' )
288 wantq = initq .OR.
lsame( compq,
'V' )
289 initz =
lsame( compz,
'I' )
290 wantz = initz .OR.
lsame( compz,
'V' )
291 lquery = ( lwork.EQ.-1 )
293 IF( .NOT.
lsame( compq,
'N' ) .AND. .NOT.wantq )
THEN
295 ELSE IF( .NOT.
lsame( compz,
'N' ) .AND. .NOT.wantz )
THEN
297 ELSE IF( n.LT.0 )
THEN
299 ELSE IF( ilo.LT.1 )
THEN
301 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
303 ELSE IF( lda.LT.max( 1, n ) )
THEN
305 ELSE IF( ldb.LT.max( 1, n ) )
THEN
307 ELSE IF( ( wantq .AND. ldq.LT.n ) .OR. ldq.LT.1 )
THEN
309 ELSE IF( ( wantz .AND. ldz.LT.n ) .OR. ldz.LT.1 )
THEN
311 ELSE IF( lwork.LT.1 .AND. .NOT.lquery )
THEN
315 CALL xerbla(
'CGGHD3', -info )
317 ELSE IF( lquery )
THEN
324 $
CALL claset(
'All', n, n, czero, cone, q, ldq )
326 $
CALL claset(
'All', n, n, czero, cone, z, ldz )
331 $
CALL claset(
'Lower', n-1, n-1, czero, czero, b(2, 1), ldb )
343 nbmin =
ilaenv( 2,
'CGGHD3',
' ', n, ilo, ihi, -1 )
344 IF( nb.GT.1 .AND. nb.LT.nh )
THEN
348 nx = max( nb,
ilaenv( 3,
'CGGHD3',
' ', n, ilo, ihi, -1 ) )
353 IF( lwork.LT.lwkopt )
THEN
359 nbmin = max( 2,
ilaenv( 2,
'CGGHD3',
' ', n, ilo, ihi,
361 IF( lwork.GE.6*n*nbmin )
THEN
370 IF( nb.LT.nbmin .OR. nb.GE.nh )
THEN
380 kacc22 =
ilaenv( 16,
'CGGHD3',
' ', n, ilo, ihi, -1 )
382 DO jcol = ilo, ihi-2, nb
383 nnb = min( nb, ihi-jcol-1 )
391 n2nb = ( ihi-jcol-1 ) / nnb - 1
392 nblst = ihi - jcol - n2nb*nnb
393 CALL claset(
'All', nblst, nblst, czero, cone, work, nblst )
394 pw = nblst * nblst + 1
396 CALL claset(
'All', 2*nnb, 2*nnb, czero, cone,
397 $ work( pw ), 2*nnb )
403 DO j = jcol, jcol+nnb-1
410 CALL clartg( temp, a( i, j ), c, s, a( i-1, j ) )
411 a( i, j ) = cmplx( c )
417 ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
419 jrow = j + n2nb*nnb + 2
423 DO jj = ppw, ppw+len-1
424 temp = work( jj + nblst )
425 work( jj + nblst ) = ctemp*temp - s*work( jj )
426 work( jj ) = conjg( s )*temp + ctemp*work( jj )
429 ppw = ppw - nblst - 1
432 ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
434 DO jrow = j0, j+2, -nnb
437 DO i = jrow+nnb-1, jrow, -1
440 DO jj = ppw, ppw+len-1
441 temp = work( jj + 2*nnb )
442 work( jj + 2*nnb ) = ctemp*temp - s*work( jj )
443 work( jj ) = conjg( s )*temp + ctemp*work( jj )
446 ppw = ppw - 2*nnb - 1
448 ppwo = ppwo + 4*nnb*nnb
467 DO i = min( jj+1, ihi ), j+2, -1
471 b( i, jj ) = ctemp*temp - conjg( s )*b( i-1, jj )
472 b( i-1, jj ) = s*temp + ctemp*b( i-1, jj )
478 temp = b( jj+1, jj+1 )
479 CALL clartg( temp, b( jj+1, jj ), c, s,
481 b( jj+1, jj ) = czero
482 CALL crot( jj-top, b( top+1, jj+1 ), 1,
483 $ b( top+1, jj ), 1, c, s )
484 a( jj+1, j ) = cmplx( c )
485 b( jj+1, j ) = -conjg( s )
491 jj = mod( ihi-j-1, 3 )
492 DO i = ihi-j-3, jj+1, -3
493 ctemp = a( j+1+i, j )
502 temp1 = a( k, j+i+1 )
503 temp2 = a( k, j+i+2 )
504 temp3 = a( k, j+i+3 )
505 a( k, j+i+3 ) = c2*temp3 + conjg( s2 )*temp2
506 temp2 = -s2*temp3 + c2*temp2
507 a( k, j+i+2 ) = c1*temp2 + conjg( s1 )*temp1
508 temp1 = -s1*temp2 + c1*temp1
509 a( k, j+i+1 ) = ctemp*temp1 + conjg( s )*temp
510 a( k, j+i ) = -s*temp1 + ctemp*temp
516 c = dble( a( j+1+i, j ) )
517 CALL crot( ihi-top, a( top+1, j+i+1 ), 1,
518 $ a( top+1, j+i ), 1, c,
519 $ -conjg( b( j+1+i, j ) ) )
525 IF ( j .LT. jcol + nnb - 1 )
THEN
538 jrow = ihi - nblst + 1
539 CALL cgemv(
'Conjugate', nblst, len, cone, work,
540 $ nblst, a( jrow, j+1 ), 1, czero,
543 DO i = jrow, jrow+nblst-len-1
544 work( ppw ) = a( i, j+1 )
547 CALL ctrmv(
'Lower',
'Conjugate',
'Non-unit',
548 $ nblst-len, work( len*nblst + 1 ), nblst,
549 $ work( pw+len ), 1 )
550 CALL cgemv(
'Conjugate', len, nblst-len, cone,
551 $ work( (len+1)*nblst - len + 1 ), nblst,
552 $ a( jrow+nblst-len, j+1 ), 1, cone,
553 $ work( pw+len ), 1 )
555 DO i = jrow, jrow+nblst-1
556 a( i, j+1 ) = work( ppw )
573 ppwo = 1 + nblst*nblst
575 DO jrow = j0, jcol+1, -nnb
577 DO i = jrow, jrow+nnb-1
578 work( ppw ) = a( i, j+1 )
582 DO i = jrow+nnb, jrow+nnb+len-1
583 work( ppw ) = a( i, j+1 )
586 CALL ctrmv(
'Upper',
'Conjugate',
'Non-unit', len,
587 $ work( ppwo + nnb ), 2*nnb, work( pw ),
589 CALL ctrmv(
'Lower',
'Conjugate',
'Non-unit', nnb,
590 $ work( ppwo + 2*len*nnb ),
591 $ 2*nnb, work( pw + len ), 1 )
592 CALL cgemv(
'Conjugate', nnb, len, cone,
593 $ work( ppwo ), 2*nnb, a( jrow, j+1 ), 1,
594 $ cone, work( pw ), 1 )
595 CALL cgemv(
'Conjugate', len, nnb, cone,
596 $ work( ppwo + 2*len*nnb + nnb ), 2*nnb,
597 $ a( jrow+nnb, j+1 ), 1, cone,
598 $ work( pw+len ), 1 )
600 DO i = jrow, jrow+len+nnb-1
601 a( i, j+1 ) = work( ppw )
604 ppwo = ppwo + 4*nnb*nnb
611 cola = n - jcol - nnb + 1
613 CALL cgemm(
'Conjugate',
'No Transpose', nblst,
614 $ cola, nblst, cone, work, nblst,
615 $ a( j, jcol+nnb ), lda, czero, work( pw ),
617 CALL clacpy(
'All', nblst, cola, work( pw ), nblst,
618 $ a( j, jcol+nnb ), lda )
619 ppwo = nblst*nblst + 1
621 DO j = j0, jcol+1, -nnb
633 CALL cunm22(
'Left',
'Conjugate', 2*nnb, cola, nnb,
634 $ nnb, work( ppwo ), 2*nnb,
635 $ a( j, jcol+nnb ), lda, work( pw ),
641 CALL cgemm(
'Conjugate',
'No Transpose', 2*nnb,
642 $ cola, 2*nnb, cone, work( ppwo ), 2*nnb,
643 $ a( j, jcol+nnb ), lda, czero, work( pw ),
645 CALL clacpy(
'All', 2*nnb, cola, work( pw ), 2*nnb,
646 $ a( j, jcol+nnb ), lda )
648 ppwo = ppwo + 4*nnb*nnb
656 topq = max( 2, j - jcol + 1 )
662 CALL cgemm(
'No Transpose',
'No Transpose', nh,
663 $ nblst, nblst, cone, q( topq, j ), ldq,
664 $ work, nblst, czero, work( pw ), nh )
665 CALL clacpy(
'All', nh, nblst, work( pw ), nh,
666 $ q( topq, j ), ldq )
667 ppwo = nblst*nblst + 1
669 DO j = j0, jcol+1, -nnb
671 topq = max( 2, j - jcol + 1 )
678 CALL cunm22(
'Right',
'No Transpose', nh, 2*nnb,
679 $ nnb, nnb, work( ppwo ), 2*nnb,
680 $ q( topq, j ), ldq, work( pw ),
686 CALL cgemm(
'No Transpose',
'No Transpose', nh,
687 $ 2*nnb, 2*nnb, cone, q( topq, j ), ldq,
688 $ work( ppwo ), 2*nnb, czero, work( pw ),
690 CALL clacpy(
'All', nh, 2*nnb, work( pw ), nh,
691 $ q( topq, j ), ldq )
693 ppwo = ppwo + 4*nnb*nnb
699 IF ( wantz .OR. top.GT.0 )
THEN
704 CALL claset(
'All', nblst, nblst, czero, cone, work,
706 pw = nblst * nblst + 1
708 CALL claset(
'All', 2*nnb, 2*nnb, czero, cone,
709 $ work( pw ), 2*nnb )
715 DO j = jcol, jcol+nnb-1
716 ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
718 jrow = j + n2nb*nnb + 2
724 DO jj = ppw, ppw+len-1
725 temp = work( jj + nblst )
726 work( jj + nblst ) = ctemp*temp -
727 $ conjg( s )*work( jj )
728 work( jj ) = s*temp + ctemp*work( jj )
731 ppw = ppw - nblst - 1
734 ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
736 DO jrow = j0, j+2, -nnb
739 DO i = jrow+nnb-1, jrow, -1
744 DO jj = ppw, ppw+len-1
745 temp = work( jj + 2*nnb )
746 work( jj + 2*nnb ) = ctemp*temp -
747 $ conjg( s )*work( jj )
748 work( jj ) = s*temp + ctemp*work( jj )
751 ppw = ppw - 2*nnb - 1
753 ppwo = ppwo + 4*nnb*nnb
758 CALL claset(
'Lower', ihi - jcol - 1, nnb, czero, czero,
759 $ a( jcol + 2, jcol ), lda )
760 CALL claset(
'Lower', ihi - jcol - 1, nnb, czero, czero,
761 $ b( jcol + 2, jcol ), ldb )
768 CALL cgemm(
'No Transpose',
'No Transpose', top,
769 $ nblst, nblst, cone, a( 1, j ), lda,
770 $ work, nblst, czero, work( pw ), top )
771 CALL clacpy(
'All', top, nblst, work( pw ), top,
773 ppwo = nblst*nblst + 1
775 DO j = j0, jcol+1, -nnb
780 CALL cunm22(
'Right',
'No Transpose', top, 2*nnb,
781 $ nnb, nnb, work( ppwo ), 2*nnb,
782 $ a( 1, j ), lda, work( pw ),
788 CALL cgemm(
'No Transpose',
'No Transpose', top,
789 $ 2*nnb, 2*nnb, cone, a( 1, j ), lda,
790 $ work( ppwo ), 2*nnb, czero,
792 CALL clacpy(
'All', top, 2*nnb, work( pw ), top,
795 ppwo = ppwo + 4*nnb*nnb
799 CALL cgemm(
'No Transpose',
'No Transpose', top,
800 $ nblst, nblst, cone, b( 1, j ), ldb,
801 $ work, nblst, czero, work( pw ), top )
802 CALL clacpy(
'All', top, nblst, work( pw ), top,
804 ppwo = nblst*nblst + 1
806 DO j = j0, jcol+1, -nnb
811 CALL cunm22(
'Right',
'No Transpose', top, 2*nnb,
812 $ nnb, nnb, work( ppwo ), 2*nnb,
813 $ b( 1, j ), ldb, work( pw ),
819 CALL cgemm(
'No Transpose',
'No Transpose', top,
820 $ 2*nnb, 2*nnb, cone, b( 1, j ), ldb,
821 $ work( ppwo ), 2*nnb, czero,
823 CALL clacpy(
'All', top, 2*nnb, work( pw ), top,
826 ppwo = ppwo + 4*nnb*nnb
835 topq = max( 2, j - jcol + 1 )
841 CALL cgemm(
'No Transpose',
'No Transpose', nh,
842 $ nblst, nblst, cone, z( topq, j ), ldz,
843 $ work, nblst, czero, work( pw ), nh )
844 CALL clacpy(
'All', nh, nblst, work( pw ), nh,
845 $ z( topq, j ), ldz )
846 ppwo = nblst*nblst + 1
848 DO j = j0, jcol+1, -nnb
850 topq = max( 2, j - jcol + 1 )
857 CALL cunm22(
'Right',
'No Transpose', nh, 2*nnb,
858 $ nnb, nnb, work( ppwo ), 2*nnb,
859 $ z( topq, j ), ldz, work( pw ),
865 CALL cgemm(
'No Transpose',
'No Transpose', nh,
866 $ 2*nnb, 2*nnb, cone, z( topq, j ), ldz,
867 $ work( ppwo ), 2*nnb, czero, work( pw ),
869 CALL clacpy(
'All', nh, 2*nnb, work( pw ), nh,
870 $ z( topq, j ), ldz )
872 ppwo = ppwo + 4*nnb*nnb
883 IF ( jcol.NE.ilo )
THEN
891 $
CALL cgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
892 $ ldq, z, ldz, ierr )
893 work( 1 ) = cmplx( lwkopt )
subroutine clartg(F, G, CS, SN, R)
CLARTG generates a plane rotation with real cosine and complex sine.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
CGGHRD
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine ctrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRMV
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
logical function lsame(CA, CB)
LSAME
subroutine cunm22(SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC, WORK, LWORK, INFO)
CUNM22 multiplies a general matrix by a banded unitary matrix.
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...