1DLAQR3(1)           LAPACK auxiliary routine (version 3.2)           DLAQR3(1)
2
3
4

NAME

SYNOPSIS

7       SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ,
8                          Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV,  WV,
9                          LDWV, WORK, LWORK )
10
11           INTEGER        IHIZ,  ILOZ,  KBOT,  KTOP, LDH, LDT, LDV, LDWV, LDZ,
12                          LWORK, N, ND, NH, NS, NV, NW
13
14           LOGICAL        WANTT, WANTZ
15
16           DOUBLE         PRECISION H( LDH, * ), SI( * ), SR( * ), T(  LDT,  *
17                          ),  V( LDV, * ), WORK( * ), WV( LDWV, * ), Z( LDZ, *
18                          )
19
20           DOUBLE         PRECISION ZERO, ONE
21
22           PARAMETER      ( ZERO = 0.0d0, ONE = 1.0d0 )
23
24           DOUBLE         PRECISION AA, BB, BETA, CC, CS, DD, EVI,  EVK,  FOO,
25                          S, SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
26
27           INTEGER        I,  IFST,  ILST,  INFO, INFQR, J, JW, K, KCOL, KEND,
28                          KLN, KROW, KWTOP, LTOP, LWK1,  LWK2,  LWK3,  LWKOPT,
29                          NMIN
30
31           LOGICAL        BULGE, SORTED
32
33           DOUBLE         PRECISION DLAMCH
34
35           INTEGER        ILAENV
36
37           EXTERNAL       DLAMCH, ILAENV
38
39           EXTERNAL       DCOPY,   DGEHRD,   DGEMM,  DLABAD,  DLACPY,  DLAHQR,
40                          DLANV2,  DLAQR4,  DLARF,  DLARFG,  DLASET,   DORMHR,
41                          DTREXC
42
43           INTRINSIC      ABS, DBLE, INT, MAX, MIN, SQRT
44
45           JW             = MIN( NW, KBOT-KTOP+1 )
46
47           IF(            JW.LE.2 ) THEN
48
49           LWKOPT         = 1
50
51           ELSE
52
53           CALL           DGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
54
55           LWK1           = INT( WORK( 1 ) )
56
57           CALL           DORMHR(  'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V,
58                          LDV, WORK, -1, INFO )
59
60           LWK2           = INT( WORK( 1 ) )
61
62           CALL           DLAQR4( .true., .true., JW, 1, JW, T, LDT,  SR,  SI,
63                          1, JW, V, LDV, WORK, -1, INFQR )
64
65           LWK3           = INT( WORK( 1 ) )
66
67           LWKOPT         = MAX( JW+MAX( LWK1, LWK2 ), LWK3 )
68
69           END            IF
70
71           IF(            LWORK.EQ.-1 ) THEN
72
73           WORK(          1 ) = DBLE( LWKOPT )
74
75           RETURN
76
77           END            IF
78
79           NS             = 0
80
81           ND             = 0
82
83           WORK(          1 ) = ONE
84
85           IF(            KTOP.GT.KBOT ) RETURN
86
87           IF(            NW.LT.1 ) RETURN
88
89           SAFMIN         = DLAMCH( 'SAFE MINIMUM' )
90
91           SAFMAX         = ONE / SAFMIN
92
93           CALL           DLABAD( SAFMIN, SAFMAX )
94
95           ULP            = DLAMCH( 'PRECISION' )
96
97           SMLNUM         = SAFMIN*( DBLE( N ) / ULP )
98
99           JW             = MIN( NW, KBOT-KTOP+1 )
100
101           KWTOP          = KBOT - JW + 1
102
103           IF(            KWTOP.EQ.KTOP ) THEN
104
105           S              = ZERO
106
107           ELSE
108
109           S              = H( KWTOP, KWTOP-1 )
110
111           END            IF
112
113           IF(            KBOT.EQ.KWTOP ) THEN
114
115           SR(            KWTOP ) = H( KWTOP, KWTOP )
116
117           SI(            KWTOP ) = ZERO
118
119           NS             = 1
120
121           ND             = 0
122
123           IF(            ABS(  S ).LE.MAX( SMLNUM, ULP*ABS( H( KWTOP, KWTOP )
124                          ) ) ) THEN
125
126           NS             = 0
127
128           ND             = 1
129
130           IF(            KWTOP.GT.KTOP ) H( KWTOP, KWTOP-1 ) = ZERO
131
132           END            IF
133
134           WORK(          1 ) = ONE
135
136           RETURN
137
138           END            IF
139
140           CALL           DLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T,  LDT
141                          )
142
143           CALL           DCOPY(  JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ),
144                          LDT+1 )
145
146           CALL           DLASET( 'A', JW, JW, ZERO, ONE, V, LDV )
147
148           NMIN           = ILAENV( 12, 'DLAQR3', 'SV', JW, 1, JW, LWORK )
149
150           IF(            JW.GT.NMIN ) THEN
151
152           CALL           DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP
153                          ), SI( KWTOP ), 1, JW, V, LDV, WORK, LWORK, INFQR )
154
155           ELSE
156
157           CALL           DLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP
158                          ), SI( KWTOP ), 1, JW, V, LDV, INFQR )
159
160           END            IF
161
162           DO             10 J = 1, JW - 3
163
164           T(             J+2, J ) = ZERO
165
166           T(             J+3, J ) = ZERO
167
168           10             CONTINUE
169
170           IF(            JW.GT.2 ) T( JW, JW-2 ) = ZERO
171
172           NS             = JW
173
174           ILST           = INFQR + 1
175
176           20             CONTINUE
177
178           IF(            ILST.LE.NS ) THEN
179
180           IF(            NS.EQ.1 ) THEN
181
182           BULGE          = .FALSE.
183
184           ELSE
185
186           BULGE          = T( NS, NS-1 ).NE.ZERO
187
188           END            IF
189
190           IF(
191
192           FOO            = ABS( T( NS, NS ) )
193
194           IF(            FOO.EQ.ZERO ) FOO = ABS( S )
195
196           IF(            ABS( S*V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN
197
198           NS             = NS - 1
199
200           ELSE
201
202           IFST           = NS
203
204           CALL           DTREXC( 'V', JW, T, LDT, V, LDV, IFST,  ILST,  WORK,
205                          INFO )
206
207           ILST           = ILST + 1
208
209           END            IF
210
211           ELSE
212
213           FOO            = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )*
214                          SQRT( ABS( T( NS-1, NS ) ) )
215
216           IF(            FOO.EQ.ZERO ) FOO = ABS( S )
217
218           IF(            MAX( ABS( S*V( 1, NS ) ), ABS(  S*V(  1,  NS-1  )  )
219                          ).LE.  MAX( SMLNUM, ULP*FOO ) ) THEN
220
221           NS             = NS - 2
222
223           ELSE
224
225           IFST           = NS
226
227           CALL           DTREXC(  'V',  JW, T, LDT, V, LDV, IFST, ILST, WORK,
228                          INFO )
229
230           ILST           = ILST + 2
231
232           END            IF
233
234           END            IF
235
236           GO             TO 20
237
238           END            IF
239
240           IF(            NS.EQ.0 ) S = ZERO
241
242           IF(            NS.LT.JW ) THEN
243
244           SORTED         = .false.
245
246           I              = NS + 1
247
248           30             CONTINUE
249
250           IF(            SORTED ) GO TO 50
251
252           SORTED         = .true.
253
254           KEND           = I - 1
255
256           I              = INFQR + 1
257
258           IF(            I.EQ.NS ) THEN
259
260           K              = I + 1
261
262           ELSE           IF( T( I+1, I ).EQ.ZERO ) THEN
263
264           K              = I + 1
265
266           ELSE
267
268           K              = I + 2
269
270           END            IF
271
272           40             CONTINUE
273
274           IF(            K.LE.KEND ) THEN
275
276           IF(            K.EQ.I+1 ) THEN
277
278           EVI            = ABS( T( I, I ) )
279
280           ELSE
281
282           EVI            = ABS( T( I, I ) ) + SQRT( ABS( T( I+1,  I  )  )  )*
283                          SQRT( ABS( T( I, I+1 ) ) )
284
285           END            IF
286
287           IF(            K.EQ.KEND ) THEN
288
289           EVK            = ABS( T( K, K ) )
290
291           ELSE           IF( T( K+1, K ).EQ.ZERO ) THEN
292
293           EVK            = ABS( T( K, K ) )
294
295           ELSE
296
297           EVK            =  ABS(  T(  K,  K ) ) + SQRT( ABS( T( K+1, K ) ) )*
298                          SQRT( ABS( T( K, K+1 ) ) )
299
300           END            IF
301
302           IF(            EVI.GE.EVK ) THEN
303
304           I              = K
305
306           ELSE
307
308           SORTED         = .false.
309
310           IFST           = I
311
312           ILST           = K
313
314           CALL           DTREXC( 'V', JW, T, LDT, V, LDV, IFST,  ILST,  WORK,
315                          INFO )
316
317           IF(            INFO.EQ.0 ) THEN
318
319           I              = ILST
320
321           ELSE
322
323           I              = K
324
325           END            IF
326
327           END            IF
328
329           IF(            I.EQ.KEND ) THEN
330
331           K              = I + 1
332
333           ELSE           IF( T( I+1, I ).EQ.ZERO ) THEN
334
335           K              = I + 1
336
337           ELSE
338
339           K              = I + 2
340
341           END            IF
342
343           GO             TO 40
344
345           END            IF
346
347           GO             TO 30
348
349           50             CONTINUE
350
351           END            IF
352
353           I              = JW
354
355           60             CONTINUE
356
357           IF(            I.GE.INFQR+1 ) THEN
358
359           IF(            I.EQ.INFQR+1 ) THEN
360
361           SR(            KWTOP+I-1 ) = T( I, I )
362
363           SI(            KWTOP+I-1 ) = ZERO
364
365           I              = I - 1
366
367           ELSE           IF( T( I, I-1 ).EQ.ZERO ) THEN
368
369           SR(            KWTOP+I-1 ) = T( I, I )
370
371           SI(            KWTOP+I-1 ) = ZERO
372
373           I              = I - 1
374
375           ELSE
376
377           AA             = T( I-1, I-1 )
378
379           CC             = T( I, I-1 )
380
381           BB             = T( I-1, I )
382
383           DD             = T( I, I )
384
385           CALL           DLANV2(  AA,  BB,  CC,  DD,  SR(  KWTOP+I-2  ),  SI(
386                          KWTOP+I-2 ), SR( KWTOP+I-1 ), SI( KWTOP+I-1  ),  CS,
387                          SN )
388
389           I              = I - 2
390
391           END            IF
392
393           GO             TO 60
394
395           END            IF
396
397           IF(            NS.LT.JW .OR. S.EQ.ZERO ) THEN
398
399           IF(            NS.GT.1 .AND. S.NE.ZERO ) THEN
400
401           CALL           DCOPY( NS, V, LDV, WORK, 1 )
402
403           BETA           = WORK( 1 )
404
405           CALL           DLARFG( NS, BETA, WORK( 2 ), 1, TAU )
406
407           WORK(          1 ) = ONE
408
409           CALL           DLASET(  'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT
410                          )
411
412           CALL           DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT, WORK( JW+1
413                          ) )
414
415           CALL           DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, WORK( JW+1
416                          ) )
417
418           CALL           DLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, WORK( JW+1
419                          ) )
420
421           CALL           DGEHRD(  JW,  1,  NS,  T,  LDT,  WORK, WORK( JW+1 ),
422                          LWORK-JW, INFO )
423
424           END            IF
425
426           IF(            KWTOP.GT.1 ) H( KWTOP, KWTOP-1 ) = S*V( 1, 1 )
427
428           CALL           DLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ),  LDH
429                          )
430
431           CALL           DCOPY(  JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
432                          LDH+1 )
433
434           IF(            NS.GT.1 .AND. S.NE.ZERO ) CALL DORMHR( 'R', 'N', JW,
435                          NS,  1,  NS,  T,  LDT,  WORK,  V, LDV, WORK( JW+1 ),
436                          LWORK-JW, INFO )
437
438           IF(            WANTT ) THEN
439
440           LTOP           = 1
441
442           ELSE
443
444           LTOP           = KTOP
445
446           END            IF
447
448           DO             70 KROW = LTOP, KWTOP - 1, NV
449
450           KLN            = MIN( NV, KWTOP-KROW )
451
452           CALL           DGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ),
453                          LDH, V, LDV, ZERO, WV, LDWV )
454
455           CALL           DLACPY(  'A',  KLN,  JW, WV, LDWV, H( KROW, KWTOP ),
456                          LDH )
457
458           70             CONTINUE
459
460           IF(            WANTT ) THEN
461
462           DO             80 KCOL = KBOT + 1, N, NH
463
464           KLN            = MIN( NH, N-KCOL+1 )
465
466           CALL           DGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, H( KWTOP,
467                          KCOL ), LDH, ZERO, T, LDT )
468
469           CALL           DLACPY(  'A', JW, KLN, T, LDT, H( KWTOP, KCOL ), LDH
470                          )
471
472           80             CONTINUE
473
474           END            IF
475
476           IF(            WANTZ ) THEN
477
478           DO             90 KROW = ILOZ, IHIZ, NV
479
480           KLN            = MIN( NV, IHIZ-KROW+1 )
481
482           CALL           DGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ),
483                          LDZ, V, LDV, ZERO, WV, LDWV )
484
485           CALL           DLACPY(  'A',  KLN,  JW, WV, LDWV, Z( KROW, KWTOP ),
486                          LDZ )
487
488           90             CONTINUE
489
490           END            IF
491
492           END            IF
493
494           ND             = JW - NS
495
496           NS             = NS - INFQR
497
498           WORK(          1 ) = DBLE( LWKOPT )
499
500           END
501

PURPOSE

503 LAPACK auxiliary routine (versionNo3v.e2m)ber 2008                       DLAQR3(1)