260 INTEGER ihiz, iloz, kacc22, kbot, ktop, ldh, ldu, ldv,
261 $ ldwh, ldwv, ldz, n, nh, nshfts, nv
265 COMPLEX*16 h( ldh, * ), s( * ), u( ldu, * ), v( ldv, * ),
266 $ wh( ldwh, * ), wv( ldwv, * ), z( ldz, * )
272 parameter ( zero = ( 0.0d0, 0.0d0 ),
273 $ one = ( 1.0d0, 0.0d0 ) )
274 DOUBLE PRECISION rzero, rone
275 parameter ( rzero = 0.0d0, rone = 1.0d0 )
278 COMPLEX*16 alpha, beta, cdum, refsum
279 DOUBLE PRECISION h11, h12, h21, h22, safmax, safmin, scl,
280 $ smlnum, tst1, tst2, ulp
281 INTEGER i2, i4, incol, j, j2, j4, jbot, jcol, jlen,
282 $ jrow, jtop, k, k1, kdu, kms, knz, krcol, kzs,
283 $ m, m22, mbot, mend, mstart, mtop, nbmps, ndcol,
285 LOGICAL accum, blk22, bmp22
293 INTRINSIC abs, dble, dconjg, dimag, max, min, mod
303 DOUBLE PRECISION cabs1
306 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
324 ns = nshfts - mod( nshfts, 2 )
328 safmin =
dlamch(
'SAFE MINIMUM' )
329 safmax = rone / safmin
330 CALL dlabad( safmin, safmax )
331 ulp =
dlamch(
'PRECISION' )
332 smlnum = safmin*( dble( n ) / ulp )
337 accum = ( kacc22.EQ.1 ) .OR. ( kacc22.EQ.2 )
341 blk22 = ( ns.GT.2 ) .AND. ( kacc22.EQ.2 )
346 $ h( ktop+2, ktop ) = zero
358 DO 210 incol = 3*( 1-nbmps ) + ktop - 1, kbot - 2, 3*nbmps - 2
361 $
CALL zlaset(
'ALL', kdu, kdu, zero, one, u, ldu )
375 DO 140 krcol = incol, min( incol+3*nbmps-3, kbot-2 )
384 mtop = max( 1, ( ( ktop-1 )-krcol+2 ) / 3+1 )
385 mbot = min( nbmps, ( kbot-krcol ) / 3 )
387 bmp22 = ( mbot.LT.nbmps ) .AND. ( krcol+3*( m22-1 ) ).EQ.
394 k = krcol + 3*( m-1 )
395 IF( k.EQ.ktop-1 )
THEN
396 CALL zlaqr1( 3, h( ktop, ktop ), ldh, s( 2*m-1 ),
397 $ s( 2*m ), v( 1, m ) )
399 CALL zlarfg( 3, alpha, v( 2, m ), 1, v( 1, m ) )
402 v( 2, m ) = h( k+2, k )
403 v( 3, m ) = h( k+3, k )
404 CALL zlarfg( 3, beta, v( 2, m ), 1, v( 1, m ) )
411 IF( h( k+3, k ).NE.zero .OR. h( k+3, k+1 ).NE.
412 $ zero .OR. h( k+3, k+2 ).EQ.zero )
THEN
427 CALL zlaqr1( 3, h( k+1, k+1 ), ldh, s( 2*m-1 ),
430 CALL zlarfg( 3, alpha, vt( 2 ), 1, vt( 1 ) )
431 refsum = dconjg( vt( 1 ) )*
432 $ ( h( k+1, k )+dconjg( vt( 2 ) )*
435 IF( cabs1( h( k+2, k )-refsum*vt( 2 ) )+
436 $ cabs1( refsum*vt( 3 ) ).GT.ulp*
437 $ ( cabs1( h( k, k ) )+cabs1( h( k+1,
438 $ k+1 ) )+cabs1( h( k+2, k+2 ) ) ) )
THEN
454 h( k+1, k ) = h( k+1, k ) - refsum
467 k = krcol + 3*( m22-1 )
469 IF( k.EQ.ktop-1 )
THEN
470 CALL zlaqr1( 2, h( k+1, k+1 ), ldh, s( 2*m22-1 ),
471 $ s( 2*m22 ), v( 1, m22 ) )
473 CALL zlarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
476 v( 2, m22 ) = h( k+2, k )
477 CALL zlarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
486 jbot = min( ndcol, kbot )
487 ELSE IF( wantt )
THEN
492 DO 30 j = max( ktop, krcol ), jbot
493 mend = min( mbot, ( j-krcol+2 ) / 3 )
495 k = krcol + 3*( m-1 )
496 refsum = dconjg( v( 1, m ) )*
497 $ ( h( k+1, j )+dconjg( v( 2, m ) )*
498 $ h( k+2, j )+dconjg( v( 3, m ) )*h( k+3, j ) )
499 h( k+1, j ) = h( k+1, j ) - refsum
500 h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m )
501 h( k+3, j ) = h( k+3, j ) - refsum*v( 3, m )
505 k = krcol + 3*( m22-1 )
506 DO 40 j = max( k+1, ktop ), jbot
507 refsum = dconjg( v( 1, m22 ) )*
508 $ ( h( k+1, j )+dconjg( v( 2, m22 ) )*
510 h( k+1, j ) = h( k+1, j ) - refsum
511 h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m22 )
520 jtop = max( ktop, incol )
521 ELSE IF( wantt )
THEN
527 IF( v( 1, m ).NE.zero )
THEN
528 k = krcol + 3*( m-1 )
529 DO 50 j = jtop, min( kbot, k+3 )
530 refsum = v( 1, m )*( h( j, k+1 )+v( 2, m )*
531 $ h( j, k+2 )+v( 3, m )*h( j, k+3 ) )
532 h( j, k+1 ) = h( j, k+1 ) - refsum
533 h( j, k+2 ) = h( j, k+2 ) -
534 $ refsum*dconjg( v( 2, m ) )
535 h( j, k+3 ) = h( j, k+3 ) -
536 $ refsum*dconjg( v( 3, m ) )
546 DO 60 j = max( 1, ktop-incol ), kdu
547 refsum = v( 1, m )*( u( j, kms+1 )+v( 2, m )*
548 $ u( j, kms+2 )+v( 3, m )*u( j, kms+3 ) )
549 u( j, kms+1 ) = u( j, kms+1 ) - refsum
550 u( j, kms+2 ) = u( j, kms+2 ) -
551 $ refsum*dconjg( v( 2, m ) )
552 u( j, kms+3 ) = u( j, kms+3 ) -
553 $ refsum*dconjg( v( 3, m ) )
555 ELSE IF( wantz )
THEN
562 refsum = v( 1, m )*( z( j, k+1 )+v( 2, m )*
563 $ z( j, k+2 )+v( 3, m )*z( j, k+3 ) )
564 z( j, k+1 ) = z( j, k+1 ) - refsum
565 z( j, k+2 ) = z( j, k+2 ) -
566 $ refsum*dconjg( v( 2, m ) )
567 z( j, k+3 ) = z( j, k+3 ) -
568 $ refsum*dconjg( v( 3, m ) )
576 k = krcol + 3*( m22-1 )
578 IF ( v( 1, m22 ).NE.zero )
THEN
579 DO 90 j = jtop, min( kbot, k+3 )
580 refsum = v( 1, m22 )*( h( j, k+1 )+v( 2, m22 )*
582 h( j, k+1 ) = h( j, k+1 ) - refsum
583 h( j, k+2 ) = h( j, k+2 ) -
584 $ refsum*dconjg( v( 2, m22 ) )
589 DO 100 j = max( 1, ktop-incol ), kdu
590 refsum = v( 1, m22 )*( u( j, kms+1 )+
591 $ v( 2, m22 )*u( j, kms+2 ) )
592 u( j, kms+1 ) = u( j, kms+1 ) - refsum
593 u( j, kms+2 ) = u( j, kms+2 ) -
594 $ refsum*dconjg( v( 2, m22 ) )
596 ELSE IF( wantz )
THEN
597 DO 110 j = iloz, ihiz
598 refsum = v( 1, m22 )*( z( j, k+1 )+v( 2, m22 )*
600 z( j, k+1 ) = z( j, k+1 ) - refsum
601 z( j, k+2 ) = z( j, k+2 ) -
602 $ refsum*dconjg( v( 2, m22 ) )
611 IF( krcol+3*( mstart-1 ).LT.ktop )
612 $ mstart = mstart + 1
616 IF( krcol.EQ.kbot-2 )
618 DO 120 m = mstart, mend
619 k = min( kbot-1, krcol+3*( m-1 ) )
630 IF( h( k+1, k ).NE.zero )
THEN
631 tst1 = cabs1( h( k, k ) ) + cabs1( h( k+1, k+1 ) )
632 IF( tst1.EQ.rzero )
THEN
634 $ tst1 = tst1 + cabs1( h( k, k-1 ) )
636 $ tst1 = tst1 + cabs1( h( k, k-2 ) )
638 $ tst1 = tst1 + cabs1( h( k, k-3 ) )
640 $ tst1 = tst1 + cabs1( h( k+2, k+1 ) )
642 $ tst1 = tst1 + cabs1( h( k+3, k+1 ) )
644 $ tst1 = tst1 + cabs1( h( k+4, k+1 ) )
646 IF( cabs1( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
648 h12 = max( cabs1( h( k+1, k ) ),
649 $ cabs1( h( k, k+1 ) ) )
650 h21 = min( cabs1( h( k+1, k ) ),
651 $ cabs1( h( k, k+1 ) ) )
652 h11 = max( cabs1( h( k+1, k+1 ) ),
653 $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
654 h22 = min( cabs1( h( k+1, k+1 ) ),
655 $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
657 tst2 = h22*( h11 / scl )
659 IF( tst2.EQ.rzero .OR. h21*( h12 / scl ).LE.
660 $ max( smlnum, ulp*tst2 ) )h( k+1, k ) = zero
667 mend = min( nbmps, ( kbot-krcol-1 ) / 3 )
668 DO 130 m = mtop, mend
669 k = krcol + 3*( m-1 )
670 refsum = v( 1, m )*v( 3, m )*h( k+4, k+3 )
671 h( k+4, k+1 ) = -refsum
672 h( k+4, k+2 ) = -refsum*dconjg( v( 2, m ) )
673 h( k+4, k+3 ) = h( k+4, k+3 ) -
674 $ refsum*dconjg( v( 3, m ) )
693 IF( ( .NOT.blk22 ) .OR. ( incol.LT.ktop ) .OR.
694 $ ( ndcol.GT.kbot ) .OR. ( ns.LE.2 ) )
THEN
705 k1 = max( 1, ktop-incol )
706 nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
710 DO 150 jcol = min( ndcol, kbot ) + 1, jbot, nh
711 jlen = min( nh, jbot-jcol+1 )
712 CALL zgemm(
'C',
'N', nu, jlen, nu, one, u( k1, k1 ),
713 $ ldu, h( incol+k1, jcol ), ldh, zero, wh,
715 CALL zlacpy(
'ALL', nu, jlen, wh, ldwh,
716 $ h( incol+k1, jcol ), ldh )
721 DO 160 jrow = jtop, max( ktop, incol ) - 1, nv
722 jlen = min( nv, max( ktop, incol )-jrow )
723 CALL zgemm(
'N',
'N', jlen, nu, nu, one,
724 $ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
725 $ ldu, zero, wv, ldwv )
726 CALL zlacpy(
'ALL', jlen, nu, wv, ldwv,
727 $ h( jrow, incol+k1 ), ldh )
733 DO 170 jrow = iloz, ihiz, nv
734 jlen = min( nv, ihiz-jrow+1 )
735 CALL zgemm(
'N',
'N', jlen, nu, nu, one,
736 $ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
737 $ ldu, zero, wv, ldwv )
738 CALL zlacpy(
'ALL', jlen, nu, wv, ldwv,
739 $ z( jrow, incol+k1 ), ldz )
757 kzs = ( j4-j2 ) - ( ns+1 )
762 DO 180 jcol = min( ndcol, kbot ) + 1, jbot, nh
763 jlen = min( nh, jbot-jcol+1 )
768 CALL zlacpy(
'ALL', knz, jlen, h( incol+1+j2, jcol ),
769 $ ldh, wh( kzs+1, 1 ), ldwh )
773 CALL zlaset(
'ALL', kzs, jlen, zero, zero, wh, ldwh )
774 CALL ztrmm(
'L',
'U',
'C',
'N', knz, jlen, one,
775 $ u( j2+1, 1+kzs ), ldu, wh( kzs+1, 1 ),
780 CALL zgemm(
'C',
'N', i2, jlen, j2, one, u, ldu,
781 $ h( incol+1, jcol ), ldh, one, wh, ldwh )
785 CALL zlacpy(
'ALL', j2, jlen, h( incol+1, jcol ), ldh,
786 $ wh( i2+1, 1 ), ldwh )
790 CALL ztrmm(
'L',
'L',
'C',
'N', j2, jlen, one,
791 $ u( 1, i2+1 ), ldu, wh( i2+1, 1 ), ldwh )
795 CALL zgemm(
'C',
'N', i4-i2, jlen, j4-j2, one,
796 $ u( j2+1, i2+1 ), ldu,
797 $ h( incol+1+j2, jcol ), ldh, one,
798 $ wh( i2+1, 1 ), ldwh )
802 CALL zlacpy(
'ALL', kdu, jlen, wh, ldwh,
803 $ h( incol+1, jcol ), ldh )
808 DO 190 jrow = jtop, max( incol, ktop ) - 1, nv
809 jlen = min( nv, max( incol, ktop )-jrow )
814 CALL zlacpy(
'ALL', jlen, knz, h( jrow, incol+1+j2 ),
815 $ ldh, wv( 1, 1+kzs ), ldwv )
819 CALL zlaset(
'ALL', jlen, kzs, zero, zero, wv, ldwv )
820 CALL ztrmm(
'R',
'U',
'N',
'N', jlen, knz, one,
821 $ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
826 CALL zgemm(
'N',
'N', jlen, i2, j2, one,
827 $ h( jrow, incol+1 ), ldh, u, ldu, one, wv,
832 CALL zlacpy(
'ALL', jlen, j2, h( jrow, incol+1 ), ldh,
833 $ wv( 1, 1+i2 ), ldwv )
837 CALL ztrmm(
'R',
'L',
'N',
'N', jlen, i4-i2, one,
838 $ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ), ldwv )
842 CALL zgemm(
'N',
'N', jlen, i4-i2, j4-j2, one,
843 $ h( jrow, incol+1+j2 ), ldh,
844 $ u( j2+1, i2+1 ), ldu, one, wv( 1, 1+i2 ),
849 CALL zlacpy(
'ALL', jlen, kdu, wv, ldwv,
850 $ h( jrow, incol+1 ), ldh )
856 DO 200 jrow = iloz, ihiz, nv
857 jlen = min( nv, ihiz-jrow+1 )
862 CALL zlacpy(
'ALL', jlen, knz,
863 $ z( jrow, incol+1+j2 ), ldz,
864 $ wv( 1, 1+kzs ), ldwv )
868 CALL zlaset(
'ALL', jlen, kzs, zero, zero, wv,
870 CALL ztrmm(
'R',
'U',
'N',
'N', jlen, knz, one,
871 $ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
876 CALL zgemm(
'N',
'N', jlen, i2, j2, one,
877 $ z( jrow, incol+1 ), ldz, u, ldu, one,
882 CALL zlacpy(
'ALL', jlen, j2, z( jrow, incol+1 ),
883 $ ldz, wv( 1, 1+i2 ), ldwv )
887 CALL ztrmm(
'R',
'L',
'N',
'N', jlen, i4-i2, one,
888 $ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ),
893 CALL zgemm(
'N',
'N', jlen, i4-i2, j4-j2, one,
894 $ z( jrow, incol+1+j2 ), ldz,
895 $ u( j2+1, i2+1 ), ldu, one,
896 $ wv( 1, 1+i2 ), ldwv )
900 CALL zlacpy(
'ALL', jlen, kdu, wv, ldwv,
901 $ z( jrow, incol+1 ), ldz )
double precision function dlamch(CMACH)
DLAMCH
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
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 dlabad(SMALL, LARGE)
DLABAD
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
subroutine zlaqr1(N, H, LDH, S1, S2, V)
ZLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and spe...