1CLAQR2(1) LAPACK auxiliary routine (version 3.1) CLAQR2(1)
2
7 SUBROUTINE CLAQR2( 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, LWKOPT
33 REAL SLAMCH
35 EXTERNAL SLAMCH
37 EXTERNAL CCOPY, CGEHRD, CGEMM, CLACPY, CLAHQR, CLARF, CLARFG,
39 CLASET, CTREXC, CUNGHR, SLABAD
40 INTRINSIC ABS, AIMAG, CMPLX, CONJG, INT, MAX, MIN, REAL
42 REAL CABS1
44 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
46 JW = MIN( NW, KBOT-KTOP+1 )
48 IF( JW.LE.2 ) THEN
50 LWKOPT = 1
52 ELSE
54 CALL CGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
56 LWK1 = INT( WORK( 1 ) )
58 CALL CUNGHR( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
60 LWK2 = INT( WORK( 1 ) )
62 LWKOPT = JW + MAX( LWK1, LWK2 )
64 END IF
66 IF( LWORK.EQ.-1 ) THEN
68 WORK( 1 ) = CMPLX( LWKOPT, 0 )
70 RETURN
72 END IF
74 NS = 0
76 ND = 0
78 IF( KTOP.GT.KBOT ) RETURN
80 IF( NW.LT.1 ) RETURN
82 SAFMIN = SLAMCH( 'SAFE MINIMUM' )
84 SAFMAX = RONE / SAFMIN
86 CALL SLABAD( SAFMIN, SAFMAX )
88 ULP = SLAMCH( 'PRECISION' )
90 SMLNUM = SAFMIN*( REAL( N ) / ULP )
92 JW = MIN( NW, KBOT-KTOP+1 )
94 KWTOP = KBOT - JW + 1
96 IF( KWTOP.EQ.KTOP ) THEN
98 S = ZERO
100 ELSE
102 S = H( KWTOP, KWTOP-1 )
104 END IF
106 IF( KBOT.EQ.KWTOP ) THEN
108 SH( KWTOP ) = H( KWTOP, KWTOP )
110 NS = 1
112 ND = 0
114 IF( CABS1( S ).LE.MAX( SMLNUM, ULP*CABS1( H( KWTOP,
116 KWTOP ) ) ) ) THEN
117 NS = 0
119 ND = 1
121 IF( KWTOP.GT.KTOP ) H( KWTOP, KWTOP-1 ) = ZERO
123 END IF
125 RETURN
127 END IF
129 CALL CLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT
131 )
132 CALL CCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ),
134 LDT+1 )
135 CALL CLASET( 'A', JW, JW, ZERO, ONE, V, LDV )
137 CALL CLAHQR( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP
139 ), 1, JW, V, LDV, INFQR )
140 NS = JW
142 ILST = INFQR + 1
144 DO 10 KNT = INFQR + 1, JW
146 FOO = CABS1( T( NS, NS ) )
148 IF( FOO.EQ.RZERO ) FOO = CABS1( S )
150 IF( CABS1( S )*CABS1( V( 1, NS ) ).LE.MAX( SMLNUM,
152 ULP*FOO ) ) THEN
153 NS = NS - 1
155 ELSE
157 IFST = NS
159 CALL CTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, INFO )
161 ILST = ILST + 1
163 END IF
165 10 CONTINUE
167 IF( NS.EQ.0 ) S = ZERO
169 IF( NS.LT.JW ) THEN
171 DO 30 I = INFQR + 1, NS
173 IFST = I
175 DO 20 J = I + 1, NS
177 IF( CABS1( T( J, J ) ).GT.CABS1( T( IFST, IFST ) ) )
179 IFST = J
180 20 CONTINUE
182 ILST = I
184 IF( IFST.NE.ILST ) CALL CTREXC( 'V', JW, T, LDT, V, LDV,
186 IFST, ILST, INFO )
187 30 CONTINUE
189 END IF
191 DO 40 I = INFQR + 1, JW
193 SH( KWTOP+I-1 ) = T( I, I )
195 40 CONTINUE
197 IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
199 IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
201 CALL CCOPY( NS, V, LDV, WORK, 1 )
203 DO 50 I = 1, NS
205 WORK( I ) = CONJG( WORK( I ) )
207 50 CONTINUE
209 BETA = WORK( 1 )
211 CALL CLARFG( NS, BETA, WORK( 2 ), 1, TAU )
213 WORK( 1 ) = ONE
215 CALL CLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT
217 )
218 CALL CLARF( 'L', NS, JW, WORK, 1, CONJG( TAU ), T, LDT,
220 WORK( JW+1 ) )
221 CALL CLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, WORK( JW+1
223 ) )
224 CALL CLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, WORK( JW+1
226 ) )
227 CALL CGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
229 LWORK-JW, INFO )
230 END IF
232 IF( KWTOP.GT.1 ) H( KWTOP, KWTOP-1 ) = S*CONJG( V( 1, 1
234 ) )
235 CALL CLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH
237 )
238 CALL CCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
240 LDH+1 )
241 IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
243 CALL CUNGHR( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
245 LWORK-JW, INFO )
246 CALL CGEMM( 'N', 'N', JW, NS, NS, ONE, V, LDV, T, LDT,
248 ZERO, WV, LDWV )
249 CALL CLACPY( 'A', JW, NS, WV, LDWV, V, LDV )
251 END IF
253 IF( WANTT ) THEN
255 LTOP = 1
257 ELSE
259 LTOP = KTOP
261 END IF
263 DO 60 KROW = LTOP, KWTOP - 1, NV
265 KLN = MIN( NV, KWTOP-KROW )
267 CALL CGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ),
269 LDH, V, LDV, ZERO, WV, LDWV )
270 CALL CLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ),
272 LDH )
273 60 CONTINUE
275 IF( WANTT ) THEN
277 DO 70 KCOL = KBOT + 1, N, NH
279 KLN = MIN( NH, N-KCOL+1 )
281 CALL CGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, H( KWTOP,
283 KCOL ), LDH, ZERO, T, LDT )
284 CALL CLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ), LDH
286 )
287 70 CONTINUE
289 END IF
291 IF( WANTZ ) THEN
293 DO 80 KROW = ILOZ, IHIZ, NV
295 KLN = MIN( NV, IHIZ-KROW+1 )
297 CALL CGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ),
299 LDZ, V, LDV, ZERO, WV, LDWV )
300 CALL CLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ),
302 LDZ )
303 80 CONTINUE
305 END IF
307 END IF
309 ND = JW - NS
311 NS = NS - INFQR
313 WORK( 1 ) = CMPLX( LWKOPT, 0 )
315 END
317
319 LAPACK auxiliary routine (versionNo3v.e1m)ber 2006 CLAQR2(1)