SYNOPSIS¶
- 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 )
INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV, LDWH, LDWV,
LDZ, N, NH, NSHFTS, NV LOGICAL WANTT, WANTZ DOUBLE PRECISION H( LDH, * ),
SI( * ), SR( * ), U( LDU, * ), V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ), Z(
LDZ, * ) DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
DOUBLE PRECISION ALPHA, BETA, H11, H12, H21, H22, REFSUM, SAFMAX, SAFMIN,
SCL, SMLNUM, SWAP, TST1, TST2, ULP INTEGER I, I2, I4, INCOL, J, J2, J4,
JBOT, JCOL, JLEN, JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS, M, M22,
MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL, NS, NU LOGICAL ACCUM, BLK22, BMP22
DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH INTRINSIC ABS, DBLE, MAX, MIN, MOD
DOUBLE PRECISION VT( 3 ) EXTERNAL DGEMM, DLABAD, DLACPY, DLAQR1, DLARFG,
DLASET, DTRMM IF( NSHFTS.LT.2 ) RETURN IF( KTOP.GE.KBOT ) RETURN DO 10 I =
1, NSHFTS - 2, 2 IF( SI( I ).NE.-SI( I+1 ) ) THEN SWAP = SR( I ) SR( I ) =
SR( I+1 ) SR( I+1 ) = SR( I+2 ) SR( I+2 ) = SWAP SWAP = SI( I ) SI( I ) =
SI( I+1 ) SI( I+1 ) = SI( I+2 ) SI( I+2 ) = SWAP END IF 10 CONTINUE NS =
NSHFTS - MOD( NSHFTS, 2 ) SAFMIN = DLAMCH( 'SAFE MINIMUM' ) SAFMAX = ONE /
SAFMIN CALL DLABAD( SAFMIN, SAFMAX ) ULP = DLAMCH( 'PRECISION' ) SMLNUM =
SAFMIN*( DBLE( N ) / ULP ) ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 )
BLK22 = ( NS.GT.2 ) .AND. ( KACC22.EQ.2 ) IF( KTOP+2.LE.KBOT ) H( KTOP+2,
KTOP ) = ZERO NBMPS = NS / 2 KDU = 6*NBMPS - 3 DO 220 INCOL = 3*( 1-NBMPS )
+ KTOP - 1, KBOT - 2, 3*NBMPS - 2 NDCOL = INCOL + KDU IF( ACCUM ) CALL
DLASET( 'ALL', KDU, KDU, ZERO, ONE, U, LDU ) DO 150 KRCOL = INCOL, MIN(
INCOL+3*NBMPS-3, KBOT-2 ) MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 ) MBOT
= MIN( NBMPS, ( KBOT-KRCOL ) / 3 ) M22 = MBOT + 1 BMP22 = ( MBOT.LT.NBMPS )
.AND. ( KRCOL+3*( M22-1 ) ).EQ. ( KBOT-2 ) DO 20 M = MTOP, MBOT K = KRCOL +
3*( M-1 ) IF( K.EQ.KTOP-1 ) THEN CALL DLAQR1( 3, H( KTOP, KTOP ), LDH, SR(
2*M-1 ), SI( 2*M-1 ), SR( 2*M ), SI( 2*M ), V( 1, M ) ) ALPHA = V( 1, M )
CALL DLARFG( 3, ALPHA, V( 2, M ), 1, V( 1, M ) ) ELSE BETA = H( K+1, K ) V(
2, M ) = H( K+2, K ) V( 3, M ) = H( K+3, K ) CALL DLARFG( 3, BETA, V( 2, M
), 1, V( 1, M ) ) IF( H( K+3, K ).NE.ZERO .OR. H( K+3, K+1 ).NE. ZERO .OR.
H( K+3, K+2 ).EQ.ZERO ) THEN H( K+1, K ) = BETA H( K+2, K ) = ZERO H( K+3, K
) = ZERO ELSE CALL DLAQR1( 3, H( K+1, K+1 ), LDH, SR( 2*M-1 ), SI( 2*M-1 ),
SR( 2*M ), SI( 2*M ), VT ) ALPHA = VT( 1 ) CALL DLARFG( 3, ALPHA, VT( 2 ),
1, VT( 1 ) ) REFSUM = VT( 1 )*( H( K+1, K )+VT( 2 )* H( K+2, K ) ) IF( ABS(
H( K+2, K )-REFSUM*VT( 2 ) )+ ABS( REFSUM*VT( 3 ) ).GT.ULP* ( ABS( H( K, K )
)+ABS( H( K+1, K+1 ) )+ABS( H( K+2, K+2 ) ) ) ) THEN H( K+1, K ) = BETA H(
K+2, K ) = ZERO H( K+3, K ) = ZERO ELSE H( K+1, K ) = H( K+1, K ) - REFSUM
H( K+2, K ) = ZERO H( K+3, K ) = ZERO V( 1, M ) = VT( 1 ) V( 2, M ) = VT( 2
) V( 3, M ) = VT( 3 ) END IF END IF END IF 20 CONTINUE K = KRCOL + 3*( M22-1
) IF( BMP22 ) THEN IF( K.EQ.KTOP-1 ) THEN CALL DLAQR1( 2, H( K+1, K+1 ),
LDH, SR( 2*M22-1 ), SI( 2*M22-1 ), SR( 2*M22 ), SI( 2*M22 ), V( 1, M22 ) )
BETA = V( 1, M22 ) CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) ) ELSE
BETA = H( K+1, K ) V( 2, M22 ) = H( K+2, K ) CALL DLARFG( 2, BETA, V( 2, M22
), 1, V( 1, M22 ) ) H( K+1, K ) = BETA H( K+2, K ) = ZERO END IF END IF IF(
ACCUM ) THEN JBOT = MIN( NDCOL, KBOT ) ELSE IF( WANTT ) THEN JBOT = N ELSE
JBOT = KBOT END IF DO 40 J = MAX( KTOP, KRCOL ), JBOT MEND = MIN( MBOT, (
J-KRCOL+2 ) / 3 ) DO 30 M = MTOP, MEND K = KRCOL + 3*( M-1 ) REFSUM = V( 1,
M )*( H( K+1, J )+V( 2, M )* H( K+2, J )+V( 3, M )*H( K+3, J ) ) H( K+1, J )
= H( K+1, J ) - REFSUM H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M ) H( K+3,
J ) = H( K+3, J ) - REFSUM*V( 3, M ) 30 CONTINUE 40 CONTINUE IF( BMP22 )
THEN K = KRCOL + 3*( M22-1 ) DO 50 J = MAX( K+1, KTOP ), JBOT REFSUM = V( 1,
M22 )*( H( K+1, J )+V( 2, M22 )* H( K+2, J ) ) H( K+1, J ) = H( K+1, J ) -
REFSUM H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M22 ) 50 CONTINUE END IF IF(
ACCUM ) THEN JTOP = MAX( KTOP, INCOL ) ELSE IF( WANTT ) THEN JTOP = 1 ELSE
JTOP = KTOP END IF DO 90 M = MTOP, MBOT IF( V( 1, M ).NE.ZERO ) THEN K =
KRCOL + 3*( M-1 ) DO 60 J = JTOP, MIN( KBOT, K+3 ) REFSUM = V( 1, M )*( H(
J, K+1 )+V( 2, M )* H( J, K+2 )+V( 3, M )*H( J, K+3 ) ) H( J, K+1 ) = H( J,
K+1 ) - REFSUM H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M ) H( J, K+3 ) = H(
J, K+3 ) - REFSUM*V( 3, M ) 60 CONTINUE IF( ACCUM ) THEN KMS = K - INCOL DO
70 J = MAX( 1, KTOP-INCOL ), KDU REFSUM = V( 1, M )*( U( J, KMS+1 )+V( 2, M
)* U( J, KMS+2 )+V( 3, M )*U( J, KMS+3 ) ) U( J, KMS+1 ) = U( J, KMS+1 ) -
REFSUM U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*V( 2, M ) U( J, KMS+3 ) = U(
J, KMS+3 ) - REFSUM*V( 3, M ) 70 CONTINUE ELSE IF( WANTZ ) THEN DO 80 J =
ILOZ, IHIZ REFSUM = V( 1, M )*( Z( J, K+1 )+V( 2, M )* Z( J, K+2 )+V( 3, M
)*Z( J, K+3 ) ) Z( J, K+1 ) = Z( J, K+1 ) - REFSUM Z( J, K+2 ) = Z( J, K+2 )
- REFSUM*V( 2, M ) Z( J, K+3 ) = Z( J, K+3 ) - REFSUM*V( 3, M ) 80 CONTINUE
END IF END IF 90 CONTINUE K = KRCOL + 3*( M22-1 ) IF( BMP22 .AND. ( V( 1,
M22 ).NE.ZERO ) ) THEN DO 100 J = JTOP, MIN( KBOT, K+3 ) REFSUM = V( 1, M22
)*( H( J, K+1 )+V( 2, M22 )* H( J, K+2 ) ) H( J, K+1 ) = H( J, K+1 ) -
REFSUM H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M22 ) 100 CONTINUE IF( ACCUM
) THEN KMS = K - INCOL DO 110 J = MAX( 1, KTOP-INCOL ), KDU REFSUM = V( 1,
M22 )*( U( J, KMS+1 )+V( 2, M22 )* U( J, KMS+2 ) ) U( J, KMS+1 ) = U( J,
KMS+1 ) - REFSUM U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*V( 2, M22 ) 110
CONTINUE ELSE IF( WANTZ ) THEN DO 120 J = ILOZ, IHIZ REFSUM = V( 1, M22 )*(
Z( J, K+1 )+V( 2, M22 )* Z( J, K+2 ) ) Z( J, K+1 ) = Z( J, K+1 ) - REFSUM Z(
J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M22 ) 120 CONTINUE END IF END IF
MSTART = MTOP IF( KRCOL+3*( MSTART-1 ).LT.KTOP ) MSTART = MSTART + 1 MEND =
MBOT IF( BMP22 ) MEND = MEND + 1 IF( KRCOL.EQ.KBOT-2 ) MEND = MEND + 1 DO
130 M = MSTART, MEND K = MIN( KBOT-1, KRCOL+3*( M-1 ) ) IF( H( K+1, K
).NE.ZERO ) THEN TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) ) IF(
TST1.EQ.ZERO ) THEN IF( K.GE.KTOP+1 ) TST1 = TST1 + ABS( H( K, K-1 ) ) IF(
K.GE.KTOP+2 ) TST1 = TST1 + ABS( H( K, K-2 ) ) IF( K.GE.KTOP+3 ) TST1 = TST1
+ ABS( H( K, K-3 ) ) IF( K.LE.KBOT-2 ) TST1 = TST1 + ABS( H( K+2, K+1 ) )
IF( K.LE.KBOT-3 ) TST1 = TST1 + ABS( H( K+3, K+1 ) ) IF( K.LE.KBOT-4 ) TST1
= TST1 + ABS( H( K+4, K+1 ) ) END IF IF( ABS( H( K+1, K ) ).LE.MAX( SMLNUM,
ULP*TST1 ) ) THEN H12 = MAX( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) ) H21 =
MIN( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) ) H11 = MAX( ABS( H( K+1, K+1 )
), ABS( H( K, K )-H( K+1, K+1 ) ) ) H22 = MIN( ABS( H( K+1, K+1 ) ), ABS( H(
K, K )-H( K+1, K+1 ) ) ) SCL = H11 + H12 TST2 = H22*( H11 / SCL ) IF(
TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE. MAX( SMLNUM, ULP*TST2 ) )H( K+1, K )
= ZERO END IF END IF 130 CONTINUE MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )
DO 140 M = MTOP, MEND K = KRCOL + 3*( M-1 ) REFSUM = V( 1, M )*V( 3, M )*H(
K+4, K+3 ) H( K+4, K+1 ) = -REFSUM H( K+4, K+2 ) = -REFSUM*V( 2, M ) H( K+4,
K+3 ) = H( K+4, K+3 ) - REFSUM*V( 3, M ) 140 CONTINUE 150 CONTINUE IF( ACCUM
) THEN IF( WANTT ) THEN JTOP = 1 JBOT = N ELSE JTOP = KTOP JBOT = KBOT END
IF IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR. ( NDCOL.GT.KBOT ) .OR. (
NS.LE.2 ) ) THEN K1 = MAX( 1, KTOP-INCOL ) NU = ( KDU-MAX( 0, NDCOL-KBOT ) )
- K1 + 1 DO 160 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH JLEN = MIN( NH,
JBOT-JCOL+1 ) CALL DGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ), LDU, H(
INCOL+K1, JCOL ), LDH, ZERO, WH, LDWH ) CALL DLACPY( 'ALL', NU, JLEN, WH,
LDWH, H( INCOL+K1, JCOL ), LDH ) 160 CONTINUE DO 170 JROW = JTOP, MAX( KTOP,
INCOL ) - 1, NV JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW ) CALL DGEMM( 'N',
'N', JLEN, NU, NU, ONE, H( JROW, INCOL+K1 ), LDH, U( K1, K1 ), LDU, ZERO,
WV, LDWV ) CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV, H( JROW, INCOL+K1 ), LDH
) 170 CONTINUE IF( WANTZ ) THEN DO 180 JROW = ILOZ, IHIZ, NV JLEN = MIN( NV,
IHIZ-JROW+1 ) CALL DGEMM( 'N', 'N', JLEN, NU, NU, ONE, Z( JROW, INCOL+K1 ),
LDZ, U( K1, K1 ), LDU, ZERO, WV, LDWV ) CALL DLACPY( 'ALL', JLEN, NU, WV,
LDWV, Z( JROW, INCOL+K1 ), LDZ ) 180 CONTINUE END IF ELSE I2 = ( KDU+1 ) / 2
I4 = KDU J2 = I4 - I2 J4 = KDU KZS = ( J4-J2 ) - ( NS+1 ) KNZ = NS + 1 DO
190 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH JLEN = MIN( NH, JBOT-JCOL+1 )
CALL DLACPY( 'ALL', KNZ, JLEN, H( INCOL+1+J2, JCOL ), LDH, WH( KZS+1, 1 ),
LDWH ) CALL DLASET( 'ALL', KZS, JLEN, ZERO, ZERO, WH, LDWH ) CALL DTRMM(
'L', 'U', 'C', 'N', KNZ, JLEN, ONE, U( J2+1, 1+KZS ), LDU, WH( KZS+1, 1 ),
LDWH ) CALL DGEMM( 'C', 'N', I2, JLEN, J2, ONE, U, LDU, H( INCOL+1, JCOL ),
LDH, ONE, WH, LDWH ) CALL DLACPY( 'ALL', J2, JLEN, H( INCOL+1, JCOL ), LDH,
WH( I2+1, 1 ), LDWH ) CALL DTRMM( 'L', 'L', 'C', 'N', J2, JLEN, ONE, U( 1,
I2+1 ), LDU, WH( I2+1, 1 ), LDWH ) CALL DGEMM( 'C', 'N', I4-I2, JLEN, J4-J2,
ONE, U( J2+1, I2+1 ), LDU, H( INCOL+1+J2, JCOL ), LDH, ONE, WH( I2+1, 1 ),
LDWH ) CALL DLACPY( 'ALL', KDU, JLEN, WH, LDWH, H( INCOL+1, JCOL ), LDH )
190 CONTINUE DO 200 JROW = JTOP, MAX( INCOL, KTOP ) - 1, NV JLEN = MIN( NV,
MAX( INCOL, KTOP )-JROW ) CALL DLACPY( 'ALL', JLEN, KNZ, H( JROW, INCOL+1+J2
), LDH, WV( 1, 1+KZS ), LDWV ) CALL DLASET( 'ALL', JLEN, KZS, ZERO, ZERO,
WV, LDWV ) CALL DTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE, U( J2+1, 1+KZS ),
LDU, WV( 1, 1+KZS ), LDWV ) CALL DGEMM( 'N', 'N', JLEN, I2, J2, ONE, H(
JROW, INCOL+1 ), LDH, U, LDU, ONE, WV, LDWV ) CALL DLACPY( 'ALL', JLEN, J2,
H( JROW, INCOL+1 ), LDH, WV( 1, 1+I2 ), LDWV ) CALL DTRMM( 'R', 'L', 'N',
'N', JLEN, I4-I2, ONE, U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV ) CALL DGEMM(
'N', 'N', JLEN, I4-I2, J4-J2, ONE, H( JROW, INCOL+1+J2 ), LDH, U( J2+1, I2+1
), LDU, ONE, WV( 1, 1+I2 ), LDWV ) CALL DLACPY( 'ALL', JLEN, KDU, WV, LDWV,
H( JROW, INCOL+1 ), LDH ) 200 CONTINUE IF( WANTZ ) THEN DO 210 JROW = ILOZ,
IHIZ, NV JLEN = MIN( NV, IHIZ-JROW+1 ) CALL DLACPY( 'ALL', JLEN, KNZ, Z(
JROW, INCOL+1+J2 ), LDZ, WV( 1, 1+KZS ), LDWV ) CALL DLASET( 'ALL', JLEN,
KZS, ZERO, ZERO, WV, LDWV ) CALL DTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ), LDWV ) CALL DGEMM( 'N', 'N', JLEN,
I2, J2, ONE, Z( JROW, INCOL+1 ), LDZ, U, LDU, ONE, WV, LDWV ) CALL DLACPY(
'ALL', JLEN, J2, Z( JROW, INCOL+1 ), LDZ, WV( 1, 1+I2 ), LDWV ) CALL DTRMM(
'R', 'L', 'N', 'N', JLEN, I4-I2, ONE, U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV
) CALL DGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE, Z( JROW, INCOL+1+J2 ), LDZ,
U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ), LDWV ) CALL DLACPY( 'ALL', JLEN,
KDU, WV, LDWV, Z( JROW, INCOL+1 ), LDZ ) 210 CONTINUE END IF END IF END IF
220 CONTINUE END