1CLAQR3(1) LAPACK auxiliary routine (version 3.2) CLAQR3(1)
2
7 SUBROUTINE CLAQR3( 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 H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
17 WORK( * ), WV( LDWV, * ), Z( LDZ, * )
18 COMPLEX ZERO, ONE
20 PARAMETER ( ZERO = ( 0.0e0, 0.0e0 ), ONE = ( 1.0e0, 0.0e0 ) )
22 REAL RZERO, RONE
24 PARAMETER ( RZERO = 0.0e0, RONE = 1.0e0 )
26 COMPLEX BETA, CDUM, S, TAU
28 REAL 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 REAL SLAMCH
35 INTEGER ILAENV
37 EXTERNAL SLAMCH, ILAENV
39 EXTERNAL CCOPY, CGEHRD, CGEMM, CLACPY, CLAHQR, CLAQR4, CLARF,
41 CLARFG, CLASET, CTREXC, CUNMHR, SLABAD
42 INTRINSIC ABS, AIMAG, CMPLX, CONJG, INT, MAX, MIN, REAL
44 REAL CABS1
46 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
48 JW = MIN( NW, KBOT-KTOP+1 )
50 IF( JW.LE.2 ) THEN
52 LWKOPT = 1
54 ELSE
56 CALL CGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
58 LWK1 = INT( WORK( 1 ) )
60 CALL CUNMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V,
62 LDV, WORK, -1, INFO )
63 LWK2 = INT( WORK( 1 ) )
65 CALL CLAQR4( .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 ) = CMPLX( 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 = SLAMCH( 'SAFE MINIMUM' )
94 SAFMAX = RONE / SAFMIN
96 CALL SLABAD( SAFMIN, SAFMAX )
98 ULP = SLAMCH( 'PRECISION' )
100 SMLNUM = SAFMIN*( REAL( 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 CLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT
143 )
144 CALL CCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ),
146 LDT+1 )
147 CALL CLASET( 'A', JW, JW, ZERO, ONE, V, LDV )
149 NMIN = ILAENV( 12, 'CLAQR3', 'SV', JW, 1, JW, LWORK )
151 IF( JW.GT.NMIN ) THEN
153 CALL CLAQR4( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP
155 ), 1, JW, V, LDV, WORK, LWORK, INFQR )
156 ELSE
158 CALL CLAHQR( .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 CTREXC( '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 CTREXC( '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 CCOPY( NS, V, LDV, WORK, 1 )
226 DO 50 I = 1, NS
228 WORK( I ) = CONJG( WORK( I ) )
230 50 CONTINUE
232 BETA = WORK( 1 )
234 CALL CLARFG( NS, BETA, WORK( 2 ), 1, TAU )
236 WORK( 1 ) = ONE
238 CALL CLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT
240 )
241 CALL CLARF( 'L', NS, JW, WORK, 1, CONJG( TAU ), T, LDT,
243 WORK( JW+1 ) )
244 CALL CLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, WORK( JW+1
246 ) )
247 CALL CLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, WORK( JW+1
249 ) )
250 CALL CGEHRD( 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*CONJG( V( 1, 1
257 ) )
258 CALL CLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH
260 )
261 CALL CCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
263 LDH+1 )
264 IF( NS.GT.1 .AND. S.NE.ZERO ) CALL CUNMHR( '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 CGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ),
284 LDH, V, LDV, ZERO, WV, LDWV )
285 CALL CLACPY( '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 CGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, H( KWTOP,
298 KCOL ), LDH, ZERO, T, LDT )
299 CALL CLACPY( '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 CGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ),
314 LDZ, V, LDV, ZERO, WV, LDWV )
315 CALL CLACPY( '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 ) = CMPLX( LWKOPT, 0 )
330 END
332
334 LAPACK auxiliary routine (versionNo3v.e2m)ber 2008 CLAQR3(1)