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 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 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 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.
REAL function slamch(CMACH)
SLAMCH
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slabad(SMALL, LARGE)
SLABAD