SYNOPSIS¶
- SUBROUTINE
DLAQR2(
- WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR,
SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK )
INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV, LDZ, LWORK,
N, ND, NH, NS, NV, NW LOGICAL WANTT, WANTZ DOUBLE PRECISION H( LDH, * ), SI(
* ), SR( * ), T( LDT, * ), V( LDV, * ), WORK( * ), WV( LDWV, * ), Z( LDZ, *
) DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 ) DOUBLE
PRECISION AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S, SAFMAX, SAFMIN,
SMLNUM, SN, TAU, ULP INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2, LWKOPT LOGICAL BULGE, SORTED
DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH EXTERNAL DCOPY, DGEHRD, DGEMM,
DLABAD, DLACPY, DLAHQR, DLANV2, DLARF, DLARFG, DLASET, DORMHR, DTREXC
INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT JW = MIN( NW, KBOT-KTOP+1 ) IF(
JW.LE.2 ) THEN LWKOPT = 1 ELSE CALL DGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK,
-1, INFO ) LWK1 = INT( WORK( 1 ) ) CALL DORMHR( 'R', 'N', JW, JW, 1, JW-1,
T, LDT, WORK, V, LDV, WORK, -1, INFO ) LWK2 = INT( WORK( 1 ) ) LWKOPT = JW +
MAX( LWK1, LWK2 ) END IF IF( LWORK.EQ.-1 ) THEN WORK( 1 ) = DBLE( LWKOPT )
RETURN END IF NS = 0 ND = 0 WORK( 1 ) = ONE IF( KTOP.GT.KBOT ) RETURN IF(
NW.LT.1 ) RETURN SAFMIN = DLAMCH( 'SAFE MINIMUM' ) SAFMAX = ONE / SAFMIN
CALL DLABAD( SAFMIN, SAFMAX ) ULP = DLAMCH( 'PRECISION' ) SMLNUM = SAFMIN*(
DBLE( N ) / ULP ) JW = MIN( NW, KBOT-KTOP+1 ) KWTOP = KBOT - JW + 1 IF(
KWTOP.EQ.KTOP ) THEN S = ZERO ELSE S = H( KWTOP, KWTOP-1 ) END IF IF(
KBOT.EQ.KWTOP ) THEN SR( KWTOP ) = H( KWTOP, KWTOP ) SI( KWTOP ) = ZERO NS =
1 ND = 0 IF( ABS( S ).LE.MAX( SMLNUM, ULP*ABS( H( KWTOP, KWTOP ) ) ) ) THEN
NS = 0 ND = 1 IF( KWTOP.GT.KTOP ) H( KWTOP, KWTOP-1 ) = ZERO END IF WORK( 1
) = ONE RETURN END IF CALL DLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T,
LDT ) CALL DCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 ) CALL
DLASET( 'A', JW, JW, ZERO, ONE, V, LDV ) CALL DLAHQR( .true., .true., JW, 1,
JW, T, LDT, SR( KWTOP ), SI( KWTOP ), 1, JW, V, LDV, INFQR ) DO 10 J = 1, JW
- 3 T( J+2, J ) = ZERO T( J+3, J ) = ZERO 10 CONTINUE IF( JW.GT.2 ) T( JW,
JW-2 ) = ZERO NS = JW ILST = INFQR + 1 20 CONTINUE IF( ILST.LE.NS ) THEN IF(
NS.EQ.1 ) THEN BULGE = .FALSE. ELSE BULGE = T( NS, NS-1 ).NE.ZERO END IF IF(
FOO = ABS( T( NS, NS ) ) IF( FOO.EQ.ZERO ) FOO = ABS( S ) IF( ABS( S*V( 1,
NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN NS = NS - 1 ELSE IFST = NS CALL
DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK, INFO ) ILST = ILST + 1
END IF ELSE FOO = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )* SQRT(
ABS( T( NS-1, NS ) ) ) IF( FOO.EQ.ZERO ) FOO = ABS( S ) IF( MAX( ABS( S*V(
1, NS ) ), ABS( S*V( 1, NS-1 ) ) ).LE. MAX( SMLNUM, ULP*FOO ) ) THEN NS = NS
- 2 ELSE IFST = NS CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
INFO ) ILST = ILST + 2 END IF END IF GO TO 20 END IF IF( NS.EQ.0 ) S = ZERO
IF( NS.LT.JW ) THEN SORTED = .false. I = NS + 1 30 CONTINUE IF( SORTED ) GO
TO 50 SORTED = .true. KEND = I - 1 I = INFQR + 1 IF( I.EQ.NS ) THEN K = I +
1 ELSE IF( T( I+1, I ).EQ.ZERO ) THEN K = I + 1 ELSE K = I + 2 END IF 40
CONTINUE IF( K.LE.KEND ) THEN IF( K.EQ.I+1 ) THEN EVI = ABS( T( I, I ) )
ELSE EVI = ABS( T( I, I ) ) + SQRT( ABS( T( I+1, I ) ) )* SQRT( ABS( T( I,
I+1 ) ) ) END IF IF( K.EQ.KEND ) THEN EVK = ABS( T( K, K ) ) ELSE IF( T(
K+1, K ).EQ.ZERO ) THEN EVK = ABS( T( K, K ) ) ELSE EVK = ABS( T( K, K ) ) +
SQRT( ABS( T( K+1, K ) ) )* SQRT( ABS( T( K, K+1 ) ) ) END IF IF( EVI.GE.EVK
) THEN I = K ELSE SORTED = .false. IFST = I ILST = K CALL DTREXC( 'V', JW,
T, LDT, V, LDV, IFST, ILST, WORK, INFO ) IF( INFO.EQ.0 ) THEN I = ILST ELSE
I = K END IF END IF IF( I.EQ.KEND ) THEN K = I + 1 ELSE IF( T( I+1, I
).EQ.ZERO ) THEN K = I + 1 ELSE K = I + 2 END IF GO TO 40 END IF GO TO 30 50
CONTINUE END IF I = JW 60 CONTINUE IF( I.GE.INFQR+1 ) THEN IF( I.EQ.INFQR+1
) THEN SR( KWTOP+I-1 ) = T( I, I ) SI( KWTOP+I-1 ) = ZERO I = I - 1 ELSE IF(
T( I, I-1 ).EQ.ZERO ) THEN SR( KWTOP+I-1 ) = T( I, I ) SI( KWTOP+I-1 ) =
ZERO I = I - 1 ELSE AA = T( I-1, I-1 ) CC = T( I, I-1 ) BB = T( I-1, I ) DD
= T( I, I ) CALL DLANV2( AA, BB, CC, DD, SR( KWTOP+I-2 ), SI( KWTOP+I-2 ),
SR( KWTOP+I-1 ), SI( KWTOP+I-1 ), CS, SN ) I = I - 2 END IF GO TO 60 END IF
IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN IF( NS.GT.1 .AND. S.NE.ZERO ) THEN CALL
DCOPY( NS, V, LDV, WORK, 1 ) BETA = WORK( 1 ) CALL DLARFG( NS, BETA, WORK( 2
), 1, TAU ) WORK( 1 ) = ONE CALL DLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3,
1 ), LDT ) CALL DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT, WORK( JW+1 ) )
CALL DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, WORK( JW+1 ) ) CALL DLARF(
'R', JW, NS, WORK, 1, TAU, V, LDV, WORK( JW+1 ) ) CALL DGEHRD( JW, 1, NS, T,
LDT, WORK, WORK( JW+1 ), LWORK-JW, INFO ) END IF IF( KWTOP.GT.1 ) H( KWTOP,
KWTOP-1 ) = S*V( 1, 1 ) CALL DLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ),
LDH ) CALL DCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ), LDH+1 ) IF(
NS.GT.1 .AND. S.NE.ZERO ) CALL DORMHR( 'R', 'N', JW, NS, 1, NS, T, LDT,
WORK, V, LDV, WORK( JW+1 ), LWORK-JW, INFO ) IF( WANTT ) THEN LTOP = 1 ELSE
LTOP = KTOP END IF DO 70 KROW = LTOP, KWTOP - 1, NV KLN = MIN( NV,
KWTOP-KROW ) CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ), LDH,
V, LDV, ZERO, WV, LDWV ) CALL DLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP
), LDH ) 70 CONTINUE IF( WANTT ) THEN DO 80 KCOL = KBOT + 1, N, NH KLN =
MIN( NH, N-KCOL+1 ) CALL DGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, H(
KWTOP, KCOL ), LDH, ZERO, T, LDT ) CALL DLACPY( 'A', JW, KLN, T, LDT, H(
KWTOP, KCOL ), LDH ) 80 CONTINUE END IF IF( WANTZ ) THEN DO 90 KROW = ILOZ,
IHIZ, NV KLN = MIN( NV, IHIZ-KROW+1 ) CALL DGEMM( 'N', 'N', KLN, JW, JW,
ONE, Z( KROW, KWTOP ), LDZ, V, LDV, ZERO, WV, LDWV ) CALL DLACPY( 'A', KLN,
JW, WV, LDWV, Z( KROW, KWTOP ), LDZ ) 90 CONTINUE END IF END IF ND = JW - NS
NS = NS - INFQR WORK( 1 ) = DBLE( LWKOPT ) END