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

PURPOSE

506 LAPACK auxiliary routine (versionNo3v.e1m)ber 2006                       DLAQR3(1)