262 SUBROUTINE slaqr5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS,
263 $ SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U,
264 $ LDU, NV, WV, LDWV, NH, WH, LDWH )
272 INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
273 $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
277 REAL H( LDH, * ), SI( * ), SR( * ), U( LDU, * ),
278 $ V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ),
285 PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
288 REAL ALPHA, BETA, H11, H12, H21, H22, REFSUM,
289 $ SAFMAX, SAFMIN, SCL, SMLNUM, SWAP, T1, T2,
290 $ t3, tst1, tst2, ulp
291 INTEGER I, I2, I4, INCOL, J, JBOT, JCOL, JLEN,
292 $ JROW, JTOP, K, K1, KDU, KMS, KRCOL,
293 $ m, m22, mbot, mtop, nbmps, ndcol,
303 INTRINSIC abs, max, min, mod, real
330 DO 10 i = 1, nshfts - 2, 2
331 IF( si( i ).NE.-si( i+1 ) )
THEN
335 sr( i+1 ) = sr( i+2 )
340 si( i+1 ) = si( i+2 )
350 ns = nshfts - mod( nshfts, 2 )
354 safmin = slamch(
'SAFE MINIMUM' )
355 safmax = one / safmin
356 CALL slabad( safmin, safmax )
357 ulp = slamch(
'PRECISION' )
358 smlnum = safmin*( real( n ) / ulp )
363 accum = ( kacc22.EQ.1 ) .OR. ( kacc22.EQ.2 )
368 $ h( ktop+2, ktop ) = zero
380 DO 180 incol = ktop - 2*nbmps + 1, kbot - 2, 2*nbmps
385 jtop = max( ktop, incol )
386 ELSE IF( wantt )
THEN
394 $
CALL slaset(
'ALL', kdu, kdu, zero, one, u, ldu )
408 DO 145 krcol = incol, min( incol+2*nbmps-1, kbot-2 )
417 mtop = max( 1, ( ktop-krcol ) / 2+1 )
418 mbot = min( nbmps, ( kbot-krcol-1 ) / 2 )
420 bmp22 = ( mbot.LT.nbmps ) .AND. ( krcol+2*( m22-1 ) ).EQ.
431 k = krcol + 2*( m22-1 )
432 IF( k.EQ.ktop-1 )
THEN
433 CALL slaqr1( 2, h( k+1, k+1 ), ldh, sr( 2*m22-1 ),
434 $ si( 2*m22-1 ), sr( 2*m22 ), si( 2*m22 ),
437 CALL slarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
440 v( 2, m22 ) = h( k+2, k )
441 CALL slarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
452 DO 30 j = jtop, min( kbot, k+3 )
453 refsum = h( j, k+1 ) + v( 2, m22 )*h( j, k+2 )
454 h( j, k+1 ) = h( j, k+1 ) - refsum*t1
455 h( j, k+2 ) = h( j, k+2 ) - refsum*t2
462 jbot = min( ndcol, kbot )
463 ELSE IF( wantt )
THEN
471 refsum = h( k+1, j ) + v( 2, m22 )*h( k+2, j )
472 h( k+1, j ) = h( k+1, j ) - refsum*t1
473 h( k+2, j ) = h( k+2, j ) - refsum*t2
486 IF( h( k+1, k ).NE.zero )
THEN
487 tst1 = abs( h( k, k ) ) + abs( h( k+1, k+1 ) )
488 IF( tst1.EQ.zero )
THEN
490 $ tst1 = tst1 + abs( h( k, k-1 ) )
492 $ tst1 = tst1 + abs( h( k, k-2 ) )
494 $ tst1 = tst1 + abs( h( k, k-3 ) )
496 $ tst1 = tst1 + abs( h( k+2, k+1 ) )
498 $ tst1 = tst1 + abs( h( k+3, k+1 ) )
500 $ tst1 = tst1 + abs( h( k+4, k+1 ) )
502 IF( abs( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
504 h12 = max( abs( h( k+1, k ) ),
505 $ abs( h( k, k+1 ) ) )
506 h21 = min( abs( h( k+1, k ) ),
507 $ abs( h( k, k+1 ) ) )
508 h11 = max( abs( h( k+1, k+1 ) ),
509 $ abs( h( k, k )-h( k+1, k+1 ) ) )
510 h22 = min( abs( h( k+1, k+1 ) ),
511 $ abs( h( k, k )-h( k+1, k+1 ) ) )
513 tst2 = h22*( h11 / scl )
515 IF( tst2.EQ.zero .OR. h21*( h12 / scl ).LE.
516 $ max( smlnum, ulp*tst2 ) )
THEN
529 DO 50 j = max( 1, ktop-incol ), kdu
530 refsum = u( j, kms+1 ) + v( 2, m22 )*u( j, kms+2 )
531 u( j, kms+1 ) = u( j, kms+1 ) - refsum*t1
532 u( j, kms+2 ) = u( j, kms+2 ) - refsum*t2
534 ELSE IF( wantz )
THEN
538 refsum = z( j, k+1 )+v( 2, m22 )*z( j, k+2 )
539 z( j, k+1 ) = z( j, k+1 ) - refsum*t1
540 z( j, k+2 ) = z( j, k+2 ) - refsum*t2
547 DO 80 m = mbot, mtop, -1
548 k = krcol + 2*( m-1 )
549 IF( k.EQ.ktop-1 )
THEN
550 CALL slaqr1( 3, h( ktop, ktop ), ldh, sr( 2*m-1 ),
551 $ si( 2*m-1 ), sr( 2*m ), si( 2*m ),
554 CALL slarfg( 3, alpha, v( 2, m ), 1, v( 1, m ) )
561 refsum = v( 1, m )*v( 3, m )*h( k+3, k+2 )
562 h( k+3, k ) = -refsum
563 h( k+3, k+1 ) = -refsum*v( 2, m )
564 h( k+3, k+2 ) = h( k+3, k+2 ) - refsum*v( 3, m )
570 v( 2, m ) = h( k+2, k )
571 v( 3, m ) = h( k+3, k )
572 CALL slarfg( 3, beta, v( 2, m ), 1, v( 1, m ) )
579 IF( h( k+3, k ).NE.zero .OR. h( k+3, k+1 ).NE.
580 $ zero .OR. h( k+3, k+2 ).EQ.zero )
THEN
595 CALL slaqr1( 3, h( k+1, k+1 ), ldh, sr( 2*m-1 ),
596 $ si( 2*m-1 ), sr( 2*m ), si( 2*m ),
599 CALL slarfg( 3, alpha, vt( 2 ), 1, vt( 1 ) )
600 refsum = vt( 1 )*( h( k+1, k )+vt( 2 )*
603 IF( abs( h( k+2, k )-refsum*vt( 2 ) )+
604 $ abs( refsum*vt( 3 ) ).GT.ulp*
605 $ ( abs( h( k, k ) )+abs( h( k+1,
606 $ k+1 ) )+abs( h( k+2, k+2 ) ) ) )
THEN
622 h( k+1, k ) = h( k+1, k ) - refsum
641 DO 70 j = jtop, min( kbot, k+3 )
642 refsum = h( j, k+1 ) + v( 2, m )*h( j, k+2 )
643 $ + v( 3, m )*h( j, k+3 )
644 h( j, k+1 ) = h( j, k+1 ) - refsum*t1
645 h( j, k+2 ) = h( j, k+2 ) - refsum*t2
646 h( j, k+3 ) = h( j, k+3 ) - refsum*t3
652 refsum = h( k+1, k+1 ) + v( 2, m )*h( k+2, k+1 )
653 $ + v( 3, m )*h( k+3, k+1 )
654 h( k+1, k+1 ) = h( k+1, k+1 ) - refsum*t1
655 h( k+2, k+1 ) = h( k+2, k+1 ) - refsum*t2
656 h( k+3, k+1 ) = h( k+3, k+1 ) - refsum*t3
669 IF( h( k+1, k ).NE.zero )
THEN
670 tst1 = abs( h( k, k ) ) + abs( h( k+1, k+1 ) )
671 IF( tst1.EQ.zero )
THEN
673 $ tst1 = tst1 + abs( h( k, k-1 ) )
675 $ tst1 = tst1 + abs( h( k, k-2 ) )
677 $ tst1 = tst1 + abs( h( k, k-3 ) )
679 $ tst1 = tst1 + abs( h( k+2, k+1 ) )
681 $ tst1 = tst1 + abs( h( k+3, k+1 ) )
683 $ tst1 = tst1 + abs( h( k+4, k+1 ) )
685 IF( abs( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
687 h12 = max( abs( h( k+1, k ) ), abs( h( k, k+1 ) ) )
688 h21 = min( abs( h( k+1, k ) ), abs( h( k, k+1 ) ) )
689 h11 = max( abs( h( k+1, k+1 ) ),
690 $ abs( h( k, k )-h( k+1, k+1 ) ) )
691 h22 = min( abs( h( k+1, k+1 ) ),
692 $ abs( h( k, k )-h( k+1, k+1 ) ) )
694 tst2 = h22*( h11 / scl )
696 IF( tst2.EQ.zero .OR. h21*( h12 / scl ).LE.
697 $ max( smlnum, ulp*tst2 ) )
THEN
707 jbot = min( ndcol, kbot )
708 ELSE IF( wantt )
THEN
714 DO 100 m = mbot, mtop, -1
715 k = krcol + 2*( m-1 )
719 DO 90 j = max( ktop, krcol + 2*m ), jbot
720 refsum = h( k+1, j ) + v( 2, m )*h( k+2, j )
721 $ + v( 3, m )*h( k+3, j )
722 h( k+1, j ) = h( k+1, j ) - refsum*t1
723 h( k+2, j ) = h( k+2, j ) - refsum*t2
724 h( k+3, j ) = h( k+3, j ) - refsum*t3
736 DO 120 m = mbot, mtop, -1
737 k = krcol + 2*( m-1 )
739 i2 = max( 1, ktop-incol )
740 i2 = max( i2, kms-(krcol-incol)+1 )
741 i4 = min( kdu, krcol + 2*( mbot-1 ) - incol + 5 )
746 refsum = u( j, kms+1 ) + v( 2, m )*u( j, kms+2 )
747 $ + v( 3, m )*u( j, kms+3 )
748 u( j, kms+1 ) = u( j, kms+1 ) - refsum*t1
749 u( j, kms+2 ) = u( j, kms+2 ) - refsum*t2
750 u( j, kms+3 ) = u( j, kms+3 ) - refsum*t3
753 ELSE IF( wantz )
THEN
759 DO 140 m = mbot, mtop, -1
760 k = krcol + 2*( m-1 )
764 DO 130 j = iloz, ihiz
765 refsum = z( j, k+1 ) + v( 2, m )*z( j, k+2 )
766 $ + v( 3, m )*z( j, k+3 )
767 z( j, k+1 ) = z( j, k+1 ) - refsum*t1
768 z( j, k+2 ) = z( j, k+2 ) - refsum*t2
769 z( j, k+3 ) = z( j, k+3 ) - refsum*t3
790 k1 = max( 1, ktop-incol )
791 nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
795 DO 150 jcol = min( ndcol, kbot ) + 1, jbot, nh
796 jlen = min( nh, jbot-jcol+1 )
797 CALL sgemm(
'C',
'N', nu, jlen, nu, one, u( k1, k1 ),
798 $ ldu, h( incol+k1, jcol ), ldh, zero, wh,
800 CALL slacpy(
'ALL', nu, jlen, wh, ldwh,
801 $ h( incol+k1, jcol ), ldh )
806 DO 160 jrow = jtop, max( ktop, incol ) - 1, nv
807 jlen = min( nv, max( ktop, incol )-jrow )
808 CALL sgemm(
'N',
'N', jlen, nu, nu, one,
809 $ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
810 $ ldu, zero, wv, ldwv )
811 CALL slacpy(
'ALL', jlen, nu, wv, ldwv,
812 $ h( jrow, incol+k1 ), ldh )
818 DO 170 jrow = iloz, ihiz, nv
819 jlen = min( nv, ihiz-jrow+1 )
820 CALL sgemm(
'N',
'N', jlen, nu, nu, one,
821 $ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
822 $ ldu, zero, wv, ldwv )
823 CALL slacpy(
'ALL', jlen, nu, wv, ldwv,
824 $ z( jrow, incol+k1 ), ldz )
subroutine slabad(SMALL, LARGE)
SLABAD
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
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 strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM