258 SUBROUTINE dlaqr5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS,
259 $ sr, si, h, ldh, iloz, ihiz, z, ldz, v, ldv, u,
260 $ ldu, nv, wv, ldwv, nh, wh, ldwh )
268 INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
269 $ ldwh, ldwv, ldz, n, nh, nshfts, nv
273 DOUBLE PRECISION H( ldh, * ), SI( * ), SR( * ), U( ldu, * ),
274 $ v( ldv, * ), wh( ldwh, * ), wv( ldwv, * ),
280 DOUBLE PRECISION ZERO, ONE
281 parameter( zero = 0.0d0, one = 1.0d0 )
284 DOUBLE PRECISION ALPHA, BETA, H11, H12, H21, H22, REFSUM,
285 $ safmax, safmin, scl, smlnum, swap, tst1, tst2,
287 INTEGER I, I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN,
288 $ jrow, jtop, k, k1, kdu, kms, knz, krcol, kzs,
289 $ m, m22, mbot, mend, mstart, mtop, nbmps, ndcol,
291 LOGICAL ACCUM, BLK22, BMP22
294 DOUBLE PRECISION DLAMCH
299 INTRINSIC abs, dble, max, min, mod
302 DOUBLE PRECISION VT( 3 )
326 DO 10 i = 1, nshfts - 2, 2
327 IF( si( i ).NE.-si( i+1 ) )
THEN
331 sr( i+1 ) = sr( i+2 )
336 si( i+1 ) = si( i+2 )
346 ns = nshfts - mod( nshfts, 2 )
350 safmin = dlamch(
'SAFE MINIMUM' )
351 safmax = one / safmin
352 CALL
dlabad( safmin, safmax )
353 ulp = dlamch(
'PRECISION' )
354 smlnum = safmin*( dble( n ) / ulp )
359 accum = ( kacc22.EQ.1 ) .OR. ( kacc22.EQ.2 )
363 blk22 = ( ns.GT.2 ) .AND. ( kacc22.EQ.2 )
368 $ h( ktop+2, ktop ) = zero
380 DO 220 incol = 3*( 1-nbmps ) + ktop - 1, kbot - 2, 3*nbmps - 2
383 $ CALL
dlaset(
'ALL', kdu, kdu, zero, one, u, ldu )
397 DO 150 krcol = incol, min( incol+3*nbmps-3, kbot-2 )
406 mtop = max( 1, ( ( ktop-1 )-krcol+2 ) / 3+1 )
407 mbot = min( nbmps, ( kbot-krcol ) / 3 )
409 bmp22 = ( mbot.LT.nbmps ) .AND. ( krcol+3*( m22-1 ) ).EQ.
416 k = krcol + 3*( m-1 )
417 IF( k.EQ.ktop-1 )
THEN
418 CALL
dlaqr1( 3, h( ktop, ktop ), ldh, sr( 2*m-1 ),
419 $ si( 2*m-1 ), sr( 2*m ), si( 2*m ),
422 CALL
dlarfg( 3, alpha, v( 2, m ), 1, v( 1, m ) )
425 v( 2, m ) = h( k+2, k )
426 v( 3, m ) = h( k+3, k )
427 CALL
dlarfg( 3, beta, v( 2, m ), 1, v( 1, m ) )
434 IF( h( k+3, k ).NE.zero .OR. h( k+3, k+1 ).NE.
435 $ zero .OR. h( k+3, k+2 ).EQ.zero )
THEN
450 CALL
dlaqr1( 3, h( k+1, k+1 ), ldh, sr( 2*m-1 ),
451 $ si( 2*m-1 ), sr( 2*m ), si( 2*m ),
454 CALL
dlarfg( 3, alpha, vt( 2 ), 1, vt( 1 ) )
455 refsum = vt( 1 )*( h( k+1, k )+vt( 2 )*
458 IF( abs( h( k+2, k )-refsum*vt( 2 ) )+
459 $ abs( refsum*vt( 3 ) ).GT.ulp*
460 $ ( abs( h( k, k ) )+abs( h( k+1,
461 $ k+1 ) )+abs( h( k+2, k+2 ) ) ) )
THEN
477 h( k+1, k ) = h( k+1, k ) - refsum
490 k = krcol + 3*( m22-1 )
492 IF( k.EQ.ktop-1 )
THEN
493 CALL
dlaqr1( 2, h( k+1, k+1 ), ldh, sr( 2*m22-1 ),
494 $ si( 2*m22-1 ), sr( 2*m22 ), si( 2*m22 ),
497 CALL
dlarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
500 v( 2, m22 ) = h( k+2, k )
501 CALL
dlarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
510 jbot = min( ndcol, kbot )
511 ELSE IF( wantt )
THEN
516 DO 40 j = max( ktop, krcol ), jbot
517 mend = min( mbot, ( j-krcol+2 ) / 3 )
519 k = krcol + 3*( m-1 )
520 refsum = v( 1, m )*( h( k+1, j )+v( 2, m )*
521 $ h( k+2, j )+v( 3, m )*h( k+3, j ) )
522 h( k+1, j ) = h( k+1, j ) - refsum
523 h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m )
524 h( k+3, j ) = h( k+3, j ) - refsum*v( 3, m )
528 k = krcol + 3*( m22-1 )
529 DO 50 j = max( k+1, ktop ), jbot
530 refsum = v( 1, m22 )*( h( k+1, j )+v( 2, m22 )*
532 h( k+1, j ) = h( k+1, j ) - refsum
533 h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m22 )
542 jtop = max( ktop, incol )
543 ELSE IF( wantt )
THEN
549 IF( v( 1, m ).NE.zero )
THEN
550 k = krcol + 3*( m-1 )
551 DO 60 j = jtop, min( kbot, k+3 )
552 refsum = v( 1, m )*( h( j, k+1 )+v( 2, m )*
553 $ h( j, k+2 )+v( 3, m )*h( j, k+3 ) )
554 h( j, k+1 ) = h( j, k+1 ) - refsum
555 h( j, k+2 ) = h( j, k+2 ) - refsum*v( 2, m )
556 h( j, k+3 ) = h( j, k+3 ) - refsum*v( 3, m )
566 DO 70 j = max( 1, ktop-incol ), kdu
567 refsum = v( 1, m )*( u( j, kms+1 )+v( 2, m )*
568 $ u( j, kms+2 )+v( 3, m )*u( j, kms+3 ) )
569 u( j, kms+1 ) = u( j, kms+1 ) - refsum
570 u( j, kms+2 ) = u( j, kms+2 ) - refsum*v( 2, m )
571 u( j, kms+3 ) = u( j, kms+3 ) - refsum*v( 3, m )
573 ELSE IF( wantz )
THEN
580 refsum = v( 1, m )*( z( j, k+1 )+v( 2, m )*
581 $ z( j, k+2 )+v( 3, m )*z( j, k+3 ) )
582 z( j, k+1 ) = z( j, k+1 ) - refsum
583 z( j, k+2 ) = z( j, k+2 ) - refsum*v( 2, m )
584 z( j, k+3 ) = z( j, k+3 ) - refsum*v( 3, m )
592 k = krcol + 3*( m22-1 )
594 IF ( v( 1, m22 ).NE.zero )
THEN
595 DO 100 j = jtop, min( kbot, k+3 )
596 refsum = v( 1, m22 )*( h( j, k+1 )+v( 2, m22 )*
598 h( j, k+1 ) = h( j, k+1 ) - refsum
599 h( j, k+2 ) = h( j, k+2 ) - refsum*v( 2, m22 )
604 DO 110 j = max( 1, ktop-incol ), kdu
605 refsum = v( 1, m22 )*( u( j, kms+1 )+
606 $ v( 2, m22 )*u( j, kms+2 ) )
607 u( j, kms+1 ) = u( j, kms+1 ) - refsum
608 u( j, kms+2 ) = u( j, kms+2 ) -
611 ELSE IF( wantz )
THEN
612 DO 120 j = iloz, ihiz
613 refsum = v( 1, m22 )*( z( j, k+1 )+v( 2, m22 )*
615 z( j, k+1 ) = z( j, k+1 ) - refsum
616 z( j, k+2 ) = z( j, k+2 ) - refsum*v( 2, m22 )
625 IF( krcol+3*( mstart-1 ).LT.ktop )
626 $ mstart = mstart + 1
630 IF( krcol.EQ.kbot-2 )
632 DO 130 m = mstart, mend
633 k = min( kbot-1, krcol+3*( m-1 ) )
644 IF( h( k+1, k ).NE.zero )
THEN
645 tst1 = abs( h( k, k ) ) + abs( h( k+1, k+1 ) )
646 IF( tst1.EQ.zero )
THEN
648 $ tst1 = tst1 + abs( h( k, k-1 ) )
650 $ tst1 = tst1 + abs( h( k, k-2 ) )
652 $ tst1 = tst1 + abs( h( k, k-3 ) )
654 $ tst1 = tst1 + abs( h( k+2, k+1 ) )
656 $ tst1 = tst1 + abs( h( k+3, k+1 ) )
658 $ tst1 = tst1 + abs( h( k+4, k+1 ) )
660 IF( abs( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
662 h12 = max( abs( h( k+1, k ) ), abs( h( k, k+1 ) ) )
663 h21 = min( abs( h( k+1, k ) ), abs( h( k, k+1 ) ) )
664 h11 = max( abs( h( k+1, k+1 ) ),
665 $ abs( h( k, k )-h( k+1, k+1 ) ) )
666 h22 = min( abs( h( k+1, k+1 ) ),
667 $ abs( h( k, k )-h( k+1, k+1 ) ) )
669 tst2 = h22*( h11 / scl )
671 IF( tst2.EQ.zero .OR. h21*( h12 / scl ).LE.
672 $ max( smlnum, ulp*tst2 ) )h( k+1, k ) = zero
679 mend = min( nbmps, ( kbot-krcol-1 ) / 3 )
680 DO 140 m = mtop, mend
681 k = krcol + 3*( m-1 )
682 refsum = v( 1, m )*v( 3, m )*h( k+4, k+3 )
683 h( k+4, k+1 ) = -refsum
684 h( k+4, k+2 ) = -refsum*v( 2, m )
685 h( k+4, k+3 ) = h( k+4, k+3 ) - refsum*v( 3, m )
704 IF( ( .NOT.blk22 ) .OR. ( incol.LT.ktop ) .OR.
705 $ ( ndcol.GT.kbot ) .OR. ( ns.LE.2 ) )
THEN
716 k1 = max( 1, ktop-incol )
717 nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
721 DO 160 jcol = min( ndcol, kbot ) + 1, jbot, nh
722 jlen = min( nh, jbot-jcol+1 )
723 CALL
dgemm(
'C',
'N', nu, jlen, nu, one, u( k1, k1 ),
724 $ ldu, h( incol+k1, jcol ), ldh, zero, wh,
726 CALL
dlacpy(
'ALL', nu, jlen, wh, ldwh,
727 $ h( incol+k1, jcol ), ldh )
732 DO 170 jrow = jtop, max( ktop, incol ) - 1, nv
733 jlen = min( nv, max( ktop, incol )-jrow )
734 CALL
dgemm(
'N',
'N', jlen, nu, nu, one,
735 $ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
736 $ ldu, zero, wv, ldwv )
737 CALL
dlacpy(
'ALL', jlen, nu, wv, ldwv,
738 $ h( jrow, incol+k1 ), ldh )
744 DO 180 jrow = iloz, ihiz, nv
745 jlen = min( nv, ihiz-jrow+1 )
746 CALL
dgemm(
'N',
'N', jlen, nu, nu, one,
747 $ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
748 $ ldu, zero, wv, ldwv )
749 CALL
dlacpy(
'ALL', jlen, nu, wv, ldwv,
750 $ z( jrow, incol+k1 ), ldz )
768 kzs = ( j4-j2 ) - ( ns+1 )
773 DO 190 jcol = min( ndcol, kbot ) + 1, jbot, nh
774 jlen = min( nh, jbot-jcol+1 )
779 CALL
dlacpy(
'ALL', knz, jlen, h( incol+1+j2, jcol ),
780 $ ldh, wh( kzs+1, 1 ), ldwh )
784 CALL
dlaset(
'ALL', kzs, jlen, zero, zero, wh, ldwh )
785 CALL
dtrmm(
'L',
'U',
'C',
'N', knz, jlen, one,
786 $ u( j2+1, 1+kzs ), ldu, wh( kzs+1, 1 ),
791 CALL
dgemm(
'C',
'N', i2, jlen, j2, one, u, ldu,
792 $ h( incol+1, jcol ), ldh, one, wh, ldwh )
796 CALL
dlacpy(
'ALL', j2, jlen, h( incol+1, jcol ), ldh,
797 $ wh( i2+1, 1 ), ldwh )
801 CALL
dtrmm(
'L',
'L',
'C',
'N', j2, jlen, one,
802 $ u( 1, i2+1 ), ldu, wh( i2+1, 1 ), ldwh )
806 CALL
dgemm(
'C',
'N', i4-i2, jlen, j4-j2, one,
807 $ u( j2+1, i2+1 ), ldu,
808 $ h( incol+1+j2, jcol ), ldh, one,
809 $ wh( i2+1, 1 ), ldwh )
813 CALL
dlacpy(
'ALL', kdu, jlen, wh, ldwh,
814 $ h( incol+1, jcol ), ldh )
819 DO 200 jrow = jtop, max( incol, ktop ) - 1, nv
820 jlen = min( nv, max( incol, ktop )-jrow )
825 CALL
dlacpy(
'ALL', jlen, knz, h( jrow, incol+1+j2 ),
826 $ ldh, wv( 1, 1+kzs ), ldwv )
830 CALL
dlaset(
'ALL', jlen, kzs, zero, zero, wv, ldwv )
831 CALL
dtrmm(
'R',
'U',
'N',
'N', jlen, knz, one,
832 $ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
837 CALL
dgemm(
'N',
'N', jlen, i2, j2, one,
838 $ h( jrow, incol+1 ), ldh, u, ldu, one, wv,
843 CALL
dlacpy(
'ALL', jlen, j2, h( jrow, incol+1 ), ldh,
844 $ wv( 1, 1+i2 ), ldwv )
848 CALL
dtrmm(
'R',
'L',
'N',
'N', jlen, i4-i2, one,
849 $ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ), ldwv )
853 CALL
dgemm(
'N',
'N', jlen, i4-i2, j4-j2, one,
854 $ h( jrow, incol+1+j2 ), ldh,
855 $ u( j2+1, i2+1 ), ldu, one, wv( 1, 1+i2 ),
860 CALL
dlacpy(
'ALL', jlen, kdu, wv, ldwv,
861 $ h( jrow, incol+1 ), ldh )
867 DO 210 jrow = iloz, ihiz, nv
868 jlen = min( nv, ihiz-jrow+1 )
873 CALL
dlacpy(
'ALL', jlen, knz,
874 $ z( jrow, incol+1+j2 ), ldz,
875 $ wv( 1, 1+kzs ), ldwv )
879 CALL
dlaset(
'ALL', jlen, kzs, zero, zero, wv,
881 CALL
dtrmm(
'R',
'U',
'N',
'N', jlen, knz, one,
882 $ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
887 CALL
dgemm(
'N',
'N', jlen, i2, j2, one,
888 $ z( jrow, incol+1 ), ldz, u, ldu, one,
893 CALL
dlacpy(
'ALL', jlen, j2, z( jrow, incol+1 ),
894 $ ldz, wv( 1, 1+i2 ), ldwv )
898 CALL
dtrmm(
'R',
'L',
'N',
'N', jlen, i4-i2, one,
899 $ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ),
904 CALL
dgemm(
'N',
'N', jlen, i4-i2, j4-j2, one,
905 $ z( jrow, incol+1+j2 ), ldz,
906 $ u( j2+1, i2+1 ), ldu, one,
907 $ wv( 1, 1+i2 ), ldwv )
911 CALL
dlacpy(
'ALL', jlen, kdu, wv, ldwv,
912 $ z( jrow, incol+1 ), ldz )
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
subroutine dlaqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
DLAQR5 performs a single small-bulge multi-shift QR sweep.
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dlaqr1(N, H, LDH, SR1, SI1, SR2, SI2, V)
DLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and spe...
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM