1ZLAQR3(1) LAPACK auxiliary routine (version 3.2) ZLAQR3(1)
2
7 SUBROUTINE ZLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ,
8 Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV,
9 LDWV, WORK, LWORK )
10 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV, LDZ,
12 LWORK, N, ND, NH, NS, NV, NW
13 LOGICAL WANTT, WANTZ
15 COMPLEX*16 H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
17 WORK( * ), WV( LDWV, * ), Z( LDZ, * )
18 COMPLEX*16 ZERO, ONE
20 PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ), ONE = ( 1.0d0, 0.0d0 ) )
22 DOUBLE PRECISION RZERO, RONE
24 PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0 )
26 COMPLEX*16 BETA, CDUM, S, TAU
28 DOUBLE PRECISION FOO, SAFMAX, SAFMIN, SMLNUM, ULP
30 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN, KNT,
32 KROW, KWTOP, LTOP, LWK1, LWK2, LWK3, LWKOPT, NMIN
33 DOUBLE PRECISION DLAMCH
35 INTEGER ILAENV
37 EXTERNAL DLAMCH, ILAENV
39 EXTERNAL DLABAD, ZCOPY, ZGEHRD, ZGEMM, ZLACPY, ZLAHQR,
41 ZLAQR4, ZLARF, ZLARFG, ZLASET, ZTREXC, ZUNMHR
42 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, MAX, MIN
44 DOUBLE PRECISION CABS1
46 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
48 JW = MIN( NW, KBOT-KTOP+1 )
50 IF( JW.LE.2 ) THEN
52 LWKOPT = 1
54 ELSE
56 CALL ZGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
58 LWK1 = INT( WORK( 1 ) )
60 CALL ZUNMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V,
62 LDV, WORK, -1, INFO )
63 LWK2 = INT( WORK( 1 ) )
65 CALL ZLAQR4( .true., .true., JW, 1, JW, T, LDT, SH, 1,
67 JW, V, LDV, WORK, -1, INFQR )
68 LWK3 = INT( WORK( 1 ) )
70 LWKOPT = MAX( JW+MAX( LWK1, LWK2 ), LWK3 )
72 END IF
74 IF( LWORK.EQ.-1 ) THEN
76 WORK( 1 ) = DCMPLX( LWKOPT, 0 )
78 RETURN
80 END IF
82 NS = 0
84 ND = 0
86 WORK( 1 ) = ONE
88 IF( KTOP.GT.KBOT ) RETURN
90 IF( NW.LT.1 ) RETURN
92 SAFMIN = DLAMCH( 'SAFE MINIMUM' )
94 SAFMAX = RONE / SAFMIN
96 CALL DLABAD( SAFMIN, SAFMAX )
98 ULP = DLAMCH( 'PRECISION' )
100 SMLNUM = SAFMIN*( DBLE( N ) / ULP )
102 JW = MIN( NW, KBOT-KTOP+1 )
104 KWTOP = KBOT - JW + 1
106 IF( KWTOP.EQ.KTOP ) THEN
108 S = ZERO
110 ELSE
112 S = H( KWTOP, KWTOP-1 )
114 END IF
116 IF( KBOT.EQ.KWTOP ) THEN
118 SH( KWTOP ) = H( KWTOP, KWTOP )
120 NS = 1
122 ND = 0
124 IF( CABS1( S ).LE.MAX( SMLNUM, ULP*CABS1( H( KWTOP,
126 KWTOP ) ) ) ) THEN
127 NS = 0
129 ND = 1
131 IF( KWTOP.GT.KTOP ) H( KWTOP, KWTOP-1 ) = ZERO
133 END IF
135 WORK( 1 ) = ONE
137 RETURN
139 END IF
141 CALL ZLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT
143 )
144 CALL ZCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ),
146 LDT+1 )
147 CALL ZLASET( 'A', JW, JW, ZERO, ONE, V, LDV )
149 NMIN = ILAENV( 12, 'ZLAQR3', 'SV', JW, 1, JW, LWORK )
151 IF( JW.GT.NMIN ) THEN
153 CALL ZLAQR4( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP
155 ), 1, JW, V, LDV, WORK, LWORK, INFQR )
156 ELSE
158 CALL ZLAHQR( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP
160 ), 1, JW, V, LDV, INFQR )
161 END IF
163 NS = JW
165 ILST = INFQR + 1
167 DO 10 KNT = INFQR + 1, JW
169 FOO = CABS1( T( NS, NS ) )
171 IF( FOO.EQ.RZERO ) FOO = CABS1( S )
173 IF( CABS1( S )*CABS1( V( 1, NS ) ).LE.MAX( SMLNUM,
175 ULP*FOO ) ) THEN
176 NS = NS - 1
178 ELSE
180 IFST = NS
182 CALL ZTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, INFO )
184 ILST = ILST + 1
186 END IF
188 10 CONTINUE
190 IF( NS.EQ.0 ) S = ZERO
192 IF( NS.LT.JW ) THEN
194 DO 30 I = INFQR + 1, NS
196 IFST = I
198 DO 20 J = I + 1, NS
200 IF( CABS1( T( J, J ) ).GT.CABS1( T( IFST, IFST ) ) )
202 IFST = J
203 20 CONTINUE
205 ILST = I
207 IF( IFST.NE.ILST ) CALL ZTREXC( 'V', JW, T, LDT, V, LDV,
209 IFST, ILST, INFO )
210 30 CONTINUE
212 END IF
214 DO 40 I = INFQR + 1, JW
216 SH( KWTOP+I-1 ) = T( I, I )
218 40 CONTINUE
220 IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
222 IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
224 CALL ZCOPY( NS, V, LDV, WORK, 1 )
226 DO 50 I = 1, NS
228 WORK( I ) = DCONJG( WORK( I ) )
230 50 CONTINUE
232 BETA = WORK( 1 )
234 CALL ZLARFG( NS, BETA, WORK( 2 ), 1, TAU )
236 WORK( 1 ) = ONE
238 CALL ZLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT
240 )
241 CALL ZLARF( 'L', NS, JW, WORK, 1, DCONJG( TAU ), T, LDT,
243 WORK( JW+1 ) )
244 CALL ZLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, WORK( JW+1
246 ) )
247 CALL ZLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, WORK( JW+1
249 ) )
250 CALL ZGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
252 LWORK-JW, INFO )
253 END IF
255 IF( KWTOP.GT.1 ) H( KWTOP, KWTOP-1 ) = S*DCONJG( V( 1, 1
257 ) )
258 CALL ZLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH
260 )
261 CALL ZCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
263 LDH+1 )
264 IF( NS.GT.1 .AND. S.NE.ZERO ) CALL ZUNMHR( 'R', 'N', JW,
266 NS, 1, NS, T, LDT, WORK, V, LDV, WORK( JW+1 ),
267 LWORK-JW, INFO )
268 IF( WANTT ) THEN
270 LTOP = 1
272 ELSE
274 LTOP = KTOP
276 END IF
278 DO 60 KROW = LTOP, KWTOP - 1, NV
280 KLN = MIN( NV, KWTOP-KROW )
282 CALL ZGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ),
284 LDH, V, LDV, ZERO, WV, LDWV )
285 CALL ZLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ),
287 LDH )
288 60 CONTINUE
290 IF( WANTT ) THEN
292 DO 70 KCOL = KBOT + 1, N, NH
294 KLN = MIN( NH, N-KCOL+1 )
296 CALL ZGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, H( KWTOP,
298 KCOL ), LDH, ZERO, T, LDT )
299 CALL ZLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ), LDH
301 )
302 70 CONTINUE
304 END IF
306 IF( WANTZ ) THEN
308 DO 80 KROW = ILOZ, IHIZ, NV
310 KLN = MIN( NV, IHIZ-KROW+1 )
312 CALL ZGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ),
314 LDZ, V, LDV, ZERO, WV, LDWV )
315 CALL ZLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ),
317 LDZ )
318 80 CONTINUE
320 END IF
322 END IF
324 ND = JW - NS
326 NS = NS - INFQR
328 WORK( 1 ) = DCMPLX( LWKOPT, 0 )
330 END
332
334 LAPACK auxiliary routine (versionNo3v.e2m)ber 2008 ZLAQR3(1)