258 SUBROUTINE slaqr5( 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 REAL H( ldh, * ), SI( * ), SR( * ), U( ldu, * ),
274 $ v( ldv, * ), wh( ldwh, * ), wv( ldwv, * ),
281 parameter ( zero = 0.0e0, one = 1.0e0 )
284 REAL 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
299 INTRINSIC abs, max, min, mod, real
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 = slamch(
'SAFE MINIMUM' )
351 safmax = one / safmin
352 CALL slabad( safmin, safmax )
353 ulp = slamch(
'PRECISION' )
354 smlnum = safmin*(
REAL( 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 slaset(
'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 slaqr1( 3, h( ktop, ktop ), ldh, sr( 2*m-1 ),
419 $ si( 2*m-1 ), sr( 2*m ), si( 2*m ),
422 CALL slarfg( 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 slarfg( 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 slaqr1( 3, h( k+1, k+1 ), ldh, sr( 2*m-1 ),
451 $ si( 2*m-1 ), sr( 2*m ), si( 2*m ),
454 CALL slarfg( 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 slaqr1( 2, h( k+1, k+1 ), ldh, sr( 2*m22-1 ),
494 $ si( 2*m22-1 ), sr( 2*m22 ), si( 2*m22 ),
497 CALL slarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
500 v( 2, m22 ) = h( k+2, k )
501 CALL slarfg( 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 ) - refsum*
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 sgemm(
'C',
'N', nu, jlen, nu, one, u( k1, k1 ),
724 $ ldu, h( incol+k1, jcol ), ldh, zero, wh,
726 CALL slacpy(
'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 sgemm(
'N',
'N', jlen, nu, nu, one,
735 $ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
736 $ ldu, zero, wv, ldwv )
737 CALL slacpy(
'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 sgemm(
'N',
'N', jlen, nu, nu, one,
747 $ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
748 $ ldu, zero, wv, ldwv )
749 CALL slacpy(
'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 slacpy(
'ALL', knz, jlen, h( incol+1+j2, jcol ),
780 $ ldh, wh( kzs+1, 1 ), ldwh )
784 CALL slaset(
'ALL', kzs, jlen, zero, zero, wh, ldwh )
785 CALL strmm(
'L',
'U',
'C',
'N', knz, jlen, one,
786 $ u( j2+1, 1+kzs ), ldu, wh( kzs+1, 1 ),
791 CALL sgemm(
'C',
'N', i2, jlen, j2, one, u, ldu,
792 $ h( incol+1, jcol ), ldh, one, wh, ldwh )
796 CALL slacpy(
'ALL', j2, jlen, h( incol+1, jcol ), ldh,
797 $ wh( i2+1, 1 ), ldwh )
801 CALL strmm(
'L',
'L',
'C',
'N', j2, jlen, one,
802 $ u( 1, i2+1 ), ldu, wh( i2+1, 1 ), ldwh )
806 CALL sgemm(
'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 slacpy(
'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 slacpy(
'ALL', jlen, knz, h( jrow, incol+1+j2 ),
826 $ ldh, wv( 1, 1+kzs ), ldwv )
830 CALL slaset(
'ALL', jlen, kzs, zero, zero, wv, ldwv )
831 CALL strmm(
'R',
'U',
'N',
'N', jlen, knz, one,
832 $ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
837 CALL sgemm(
'N',
'N', jlen, i2, j2, one,
838 $ h( jrow, incol+1 ), ldh, u, ldu, one, wv,
843 CALL slacpy(
'ALL', jlen, j2, h( jrow, incol+1 ), ldh,
844 $ wv( 1, 1+i2 ), ldwv )
848 CALL strmm(
'R',
'L',
'N',
'N', jlen, i4-i2, one,
849 $ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ), ldwv )
853 CALL sgemm(
'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 slacpy(
'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 slacpy(
'ALL', jlen, knz,
874 $ z( jrow, incol+1+j2 ), ldz,
875 $ wv( 1, 1+kzs ), ldwv )
879 CALL slaset(
'ALL', jlen, kzs, zero, zero, wv,
881 CALL strmm(
'R',
'U',
'N',
'N', jlen, knz, one,
882 $ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
887 CALL sgemm(
'N',
'N', jlen, i2, j2, one,
888 $ z( jrow, incol+1 ), ldz, u, ldu, one,
893 CALL slacpy(
'ALL', jlen, j2, z( jrow, incol+1 ),
894 $ ldz, wv( 1, 1+i2 ), ldwv )
898 CALL strmm(
'R',
'L',
'N',
'N', jlen, i4-i2, one,
899 $ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ),
904 CALL sgemm(
'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 slacpy(
'ALL', jlen, kdu, wv, ldwv,
912 $ z( jrow, incol+1 ), ldz )
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
subroutine slaqr1(N, H, LDH, SR1, SI1, SR2, SI2, V)
SLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and spe...
subroutine slaqr5(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)
SLAQR5 performs a single small-bulge multi-shift QR sweep.
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.