SYNOPSIS¶
- SUBROUTINE
DLAQR3(
- 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, LWK3, LWKOPT, NMIN LOGICAL BULGE,
SORTED DOUBLE PRECISION DLAMCH INTEGER ILAENV EXTERNAL DLAMCH, ILAENV
EXTERNAL DCOPY, DGEHRD, DGEMM, DLABAD, DLACPY, DLAHQR, DLANV2, DLAQR4,
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 ) ) CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT,
SR, SI, 1, JW, V, LDV, WORK, -1, INFQR ) LWK3 = INT( WORK( 1 ) ) LWKOPT =
MAX( JW+MAX( LWK1, LWK2 ), LWK3 ) 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 ) NMIN = ILAENV( 12, 'DLAQR3', 'SV',
JW, 1, JW, LWORK ) IF( JW.GT.NMIN ) THEN CALL DLAQR4( .true., .true., JW, 1,
JW, T, LDT, SR( KWTOP ), SI( KWTOP ), 1, JW, V, LDV, WORK, LWORK, INFQR )
ELSE CALL DLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ), SI( KWTOP
), 1, JW, V, LDV, INFQR ) END IF 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