250 SUBROUTINE zlaqr5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S,
251 $ h, ldh, iloz, ihiz, z, ldz, v, ldv, u, ldu, nv,
252 $ wv, ldwv, nh, wh, ldwh )
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 )
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 dlabad(SMALL, LARGE)
DLABAD
subroutine zlaqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
ZLAQR5 performs a single small-bulge multi-shift QR sweep.
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...
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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...
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM