LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
dlaqr3.f
Go to the documentation of this file.
1 *> \brief \b DLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLAQR3 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr3.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr3.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr3.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
22 * IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T,
23 * LDT, NV, WV, LDWV, WORK, LWORK )
24 *
25 * .. Scalar Arguments ..
26 * INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
27 * $ LDZ, LWORK, N, ND, NH, NS, NV, NW
28 * LOGICAL WANTT, WANTZ
29 * ..
30 * .. Array Arguments ..
31 * DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
32 * $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
33 * $ Z( LDZ, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> Aggressive early deflation:
43 *>
44 *> DLAQR3 accepts as input an upper Hessenberg matrix
45 *> H and performs an orthogonal similarity transformation
46 *> designed to detect and deflate fully converged eigenvalues from
47 *> a trailing principal submatrix. On output H has been over-
48 *> written by a new Hessenberg matrix that is a perturbation of
49 *> an orthogonal similarity transformation of H. It is to be
50 *> hoped that the final version of H has many zero subdiagonal
51 *> entries.
52 *> \endverbatim
53 *
54 * Arguments:
55 * ==========
56 *
57 *> \param[in] WANTT
58 *> \verbatim
59 *> WANTT is LOGICAL
60 *> If .TRUE., then the Hessenberg matrix H is fully updated
61 *> so that the quasi-triangular Schur factor may be
62 *> computed (in cooperation with the calling subroutine).
63 *> If .FALSE., then only enough of H is updated to preserve
64 *> the eigenvalues.
65 *> \endverbatim
66 *>
67 *> \param[in] WANTZ
68 *> \verbatim
69 *> WANTZ is LOGICAL
70 *> If .TRUE., then the orthogonal matrix Z is updated so
71 *> so that the orthogonal Schur factor may be computed
72 *> (in cooperation with the calling subroutine).
73 *> If .FALSE., then Z is not referenced.
74 *> \endverbatim
75 *>
76 *> \param[in] N
77 *> \verbatim
78 *> N is INTEGER
79 *> The order of the matrix H and (if WANTZ is .TRUE.) the
80 *> order of the orthogonal matrix Z.
81 *> \endverbatim
82 *>
83 *> \param[in] KTOP
84 *> \verbatim
85 *> KTOP is INTEGER
86 *> It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
87 *> KBOT and KTOP together determine an isolated block
88 *> along the diagonal of the Hessenberg matrix.
89 *> \endverbatim
90 *>
91 *> \param[in] KBOT
92 *> \verbatim
93 *> KBOT is INTEGER
94 *> It is assumed without a check that either
95 *> KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together
96 *> determine an isolated block along the diagonal of the
97 *> Hessenberg matrix.
98 *> \endverbatim
99 *>
100 *> \param[in] NW
101 *> \verbatim
102 *> NW is INTEGER
103 *> Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
104 *> \endverbatim
105 *>
106 *> \param[in,out] H
107 *> \verbatim
108 *> H is DOUBLE PRECISION array, dimension (LDH,N)
109 *> On input the initial N-by-N section of H stores the
110 *> Hessenberg matrix undergoing aggressive early deflation.
111 *> On output H has been transformed by an orthogonal
112 *> similarity transformation, perturbed, and the returned
113 *> to Hessenberg form that (it is to be hoped) has some
114 *> zero subdiagonal entries.
115 *> \endverbatim
116 *>
117 *> \param[in] LDH
118 *> \verbatim
119 *> LDH is INTEGER
120 *> Leading dimension of H just as declared in the calling
121 *> subroutine. N <= LDH
122 *> \endverbatim
123 *>
124 *> \param[in] ILOZ
125 *> \verbatim
126 *> ILOZ is INTEGER
127 *> \endverbatim
128 *>
129 *> \param[in] IHIZ
130 *> \verbatim
131 *> IHIZ is INTEGER
132 *> Specify the rows of Z to which transformations must be
133 *> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
134 *> \endverbatim
135 *>
136 *> \param[in,out] Z
137 *> \verbatim
138 *> Z is DOUBLE PRECISION array, dimension (LDZ,N)
139 *> IF WANTZ is .TRUE., then on output, the orthogonal
140 *> similarity transformation mentioned above has been
141 *> accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.
142 *> If WANTZ is .FALSE., then Z is unreferenced.
143 *> \endverbatim
144 *>
145 *> \param[in] LDZ
146 *> \verbatim
147 *> LDZ is INTEGER
148 *> The leading dimension of Z just as declared in the
149 *> calling subroutine. 1 <= LDZ.
150 *> \endverbatim
151 *>
152 *> \param[out] NS
153 *> \verbatim
154 *> NS is INTEGER
155 *> The number of unconverged (ie approximate) eigenvalues
156 *> returned in SR and SI that may be used as shifts by the
157 *> calling subroutine.
158 *> \endverbatim
159 *>
160 *> \param[out] ND
161 *> \verbatim
162 *> ND is INTEGER
163 *> The number of converged eigenvalues uncovered by this
164 *> subroutine.
165 *> \endverbatim
166 *>
167 *> \param[out] SR
168 *> \verbatim
169 *> SR is DOUBLE PRECISION array, dimension (KBOT)
170 *> \endverbatim
171 *>
172 *> \param[out] SI
173 *> \verbatim
174 *> SI is DOUBLE PRECISION array, dimension (KBOT)
175 *> On output, the real and imaginary parts of approximate
176 *> eigenvalues that may be used for shifts are stored in
177 *> SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
178 *> SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
179 *> The real and imaginary parts of converged eigenvalues
180 *> are stored in SR(KBOT-ND+1) through SR(KBOT) and
181 *> SI(KBOT-ND+1) through SI(KBOT), respectively.
182 *> \endverbatim
183 *>
184 *> \param[out] V
185 *> \verbatim
186 *> V is DOUBLE PRECISION array, dimension (LDV,NW)
187 *> An NW-by-NW work array.
188 *> \endverbatim
189 *>
190 *> \param[in] LDV
191 *> \verbatim
192 *> LDV is INTEGER
193 *> The leading dimension of V just as declared in the
194 *> calling subroutine. NW <= LDV
195 *> \endverbatim
196 *>
197 *> \param[in] NH
198 *> \verbatim
199 *> NH is INTEGER
200 *> The number of columns of T. NH >= NW.
201 *> \endverbatim
202 *>
203 *> \param[out] T
204 *> \verbatim
205 *> T is DOUBLE PRECISION array, dimension (LDT,NW)
206 *> \endverbatim
207 *>
208 *> \param[in] LDT
209 *> \verbatim
210 *> LDT is INTEGER
211 *> The leading dimension of T just as declared in the
212 *> calling subroutine. NW <= LDT
213 *> \endverbatim
214 *>
215 *> \param[in] NV
216 *> \verbatim
217 *> NV is INTEGER
218 *> The number of rows of work array WV available for
219 *> workspace. NV >= NW.
220 *> \endverbatim
221 *>
222 *> \param[out] WV
223 *> \verbatim
224 *> WV is DOUBLE PRECISION array, dimension (LDWV,NW)
225 *> \endverbatim
226 *>
227 *> \param[in] LDWV
228 *> \verbatim
229 *> LDWV is INTEGER
230 *> The leading dimension of W just as declared in the
231 *> calling subroutine. NW <= LDV
232 *> \endverbatim
233 *>
234 *> \param[out] WORK
235 *> \verbatim
236 *> WORK is DOUBLE PRECISION array, dimension (LWORK)
237 *> On exit, WORK(1) is set to an estimate of the optimal value
238 *> of LWORK for the given values of N, NW, KTOP and KBOT.
239 *> \endverbatim
240 *>
241 *> \param[in] LWORK
242 *> \verbatim
243 *> LWORK is INTEGER
244 *> The dimension of the work array WORK. LWORK = 2*NW
245 *> suffices, but greater efficiency may result from larger
246 *> values of LWORK.
247 *>
248 *> If LWORK = -1, then a workspace query is assumed; DLAQR3
249 *> only estimates the optimal workspace size for the given
250 *> values of N, NW, KTOP and KBOT. The estimate is returned
251 *> in WORK(1). No error message related to LWORK is issued
252 *> by XERBLA. Neither H nor Z are accessed.
253 *> \endverbatim
254 *
255 * Authors:
256 * ========
257 *
258 *> \author Univ. of Tennessee
259 *> \author Univ. of California Berkeley
260 *> \author Univ. of Colorado Denver
261 *> \author NAG Ltd.
262 *
263 *> \ingroup doubleOTHERauxiliary
264 *
265 *> \par Contributors:
266 * ==================
267 *>
268 *> Karen Braman and Ralph Byers, Department of Mathematics,
269 *> University of Kansas, USA
270 *>
271 * =====================================================================
272  SUBROUTINE dlaqr3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
273  $ IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T,
274  $ LDT, NV, WV, LDWV, WORK, LWORK )
275 *
276 * -- LAPACK auxiliary routine --
277 * -- LAPACK is a software package provided by Univ. of Tennessee, --
278 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
279 *
280 * .. Scalar Arguments ..
281  INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
282  $ LDZ, LWORK, N, ND, NH, NS, NV, NW
283  LOGICAL WANTT, WANTZ
284 * ..
285 * .. Array Arguments ..
286  DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
287  $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
288  $ z( ldz, * )
289 * ..
290 *
291 * ================================================================
292 * .. Parameters ..
293  DOUBLE PRECISION ZERO, ONE
294  PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
295 * ..
296 * .. Local Scalars ..
297  DOUBLE PRECISION AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S,
298  $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
299  INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
300  $ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2, LWK3,
301  $ lwkopt, nmin
302  LOGICAL BULGE, SORTED
303 * ..
304 * .. External Functions ..
305  DOUBLE PRECISION DLAMCH
306  INTEGER ILAENV
307  EXTERNAL dlamch, ilaenv
308 * ..
309 * .. External Subroutines ..
310  EXTERNAL dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr,
312  $ dtrexc
313 * ..
314 * .. Intrinsic Functions ..
315  INTRINSIC abs, dble, int, max, min, sqrt
316 * ..
317 * .. Executable Statements ..
318 *
319 * ==== Estimate optimal workspace. ====
320 *
321  jw = min( nw, kbot-ktop+1 )
322  IF( jw.LE.2 ) THEN
323  lwkopt = 1
324  ELSE
325 *
326 * ==== Workspace query call to DGEHRD ====
327 *
328  CALL dgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
329  lwk1 = int( work( 1 ) )
330 *
331 * ==== Workspace query call to DORMHR ====
332 *
333  CALL dormhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
334  $ work, -1, info )
335  lwk2 = int( work( 1 ) )
336 *
337 * ==== Workspace query call to DLAQR4 ====
338 *
339  CALL dlaqr4( .true., .true., jw, 1, jw, t, ldt, sr, si, 1, jw,
340  $ v, ldv, work, -1, infqr )
341  lwk3 = int( work( 1 ) )
342 *
343 * ==== Optimal workspace ====
344 *
345  lwkopt = max( jw+max( lwk1, lwk2 ), lwk3 )
346  END IF
347 *
348 * ==== Quick return in case of workspace query. ====
349 *
350  IF( lwork.EQ.-1 ) THEN
351  work( 1 ) = dble( lwkopt )
352  RETURN
353  END IF
354 *
355 * ==== Nothing to do ...
356 * ... for an empty active block ... ====
357  ns = 0
358  nd = 0
359  work( 1 ) = one
360  IF( ktop.GT.kbot )
361  $ RETURN
362 * ... nor for an empty deflation window. ====
363  IF( nw.LT.1 )
364  $ RETURN
365 *
366 * ==== Machine constants ====
367 *
368  safmin = dlamch( 'SAFE MINIMUM' )
369  safmax = one / safmin
370  CALL dlabad( safmin, safmax )
371  ulp = dlamch( 'PRECISION' )
372  smlnum = safmin*( dble( n ) / ulp )
373 *
374 * ==== Setup deflation window ====
375 *
376  jw = min( nw, kbot-ktop+1 )
377  kwtop = kbot - jw + 1
378  IF( kwtop.EQ.ktop ) THEN
379  s = zero
380  ELSE
381  s = h( kwtop, kwtop-1 )
382  END IF
383 *
384  IF( kbot.EQ.kwtop ) THEN
385 *
386 * ==== 1-by-1 deflation window: not much to do ====
387 *
388  sr( kwtop ) = h( kwtop, kwtop )
389  si( kwtop ) = zero
390  ns = 1
391  nd = 0
392  IF( abs( s ).LE.max( smlnum, ulp*abs( h( kwtop, kwtop ) ) ) )
393  $ THEN
394  ns = 0
395  nd = 1
396  IF( kwtop.GT.ktop )
397  $ h( kwtop, kwtop-1 ) = zero
398  END IF
399  work( 1 ) = one
400  RETURN
401  END IF
402 *
403 * ==== Convert to spike-triangular form. (In case of a
404 * . rare QR failure, this routine continues to do
405 * . aggressive early deflation using that part of
406 * . the deflation window that converged using INFQR
407 * . here and there to keep track.) ====
408 *
409  CALL dlacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
410  CALL dcopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
411 *
412  CALL dlaset( 'A', jw, jw, zero, one, v, ldv )
413  nmin = ilaenv( 12, 'DLAQR3', 'SV', jw, 1, jw, lwork )
414  IF( jw.GT.nmin ) THEN
415  CALL dlaqr4( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
416  $ si( kwtop ), 1, jw, v, ldv, work, lwork, infqr )
417  ELSE
418  CALL dlahqr( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
419  $ si( kwtop ), 1, jw, v, ldv, infqr )
420  END IF
421 *
422 * ==== DTREXC needs a clean margin near the diagonal ====
423 *
424  DO 10 j = 1, jw - 3
425  t( j+2, j ) = zero
426  t( j+3, j ) = zero
427  10 CONTINUE
428  IF( jw.GT.2 )
429  $ t( jw, jw-2 ) = zero
430 *
431 * ==== Deflation detection loop ====
432 *
433  ns = jw
434  ilst = infqr + 1
435  20 CONTINUE
436  IF( ilst.LE.ns ) THEN
437  IF( ns.EQ.1 ) THEN
438  bulge = .false.
439  ELSE
440  bulge = t( ns, ns-1 ).NE.zero
441  END IF
442 *
443 * ==== Small spike tip test for deflation ====
444 *
445  IF( .NOT. bulge ) THEN
446 *
447 * ==== Real eigenvalue ====
448 *
449  foo = abs( t( ns, ns ) )
450  IF( foo.EQ.zero )
451  $ foo = abs( s )
452  IF( abs( s*v( 1, ns ) ).LE.max( smlnum, ulp*foo ) ) THEN
453 *
454 * ==== Deflatable ====
455 *
456  ns = ns - 1
457  ELSE
458 *
459 * ==== Undeflatable. Move it up out of the way.
460 * . (DTREXC can not fail in this case.) ====
461 *
462  ifst = ns
463  CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
464  $ info )
465  ilst = ilst + 1
466  END IF
467  ELSE
468 *
469 * ==== Complex conjugate pair ====
470 *
471  foo = abs( t( ns, ns ) ) + sqrt( abs( t( ns, ns-1 ) ) )*
472  $ sqrt( abs( t( ns-1, ns ) ) )
473  IF( foo.EQ.zero )
474  $ foo = abs( s )
475  IF( max( abs( s*v( 1, ns ) ), abs( s*v( 1, ns-1 ) ) ).LE.
476  $ max( smlnum, ulp*foo ) ) THEN
477 *
478 * ==== Deflatable ====
479 *
480  ns = ns - 2
481  ELSE
482 *
483 * ==== Undeflatable. Move them up out of the way.
484 * . Fortunately, DTREXC does the right thing with
485 * . ILST in case of a rare exchange failure. ====
486 *
487  ifst = ns
488  CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
489  $ info )
490  ilst = ilst + 2
491  END IF
492  END IF
493 *
494 * ==== End deflation detection loop ====
495 *
496  GO TO 20
497  END IF
498 *
499 * ==== Return to Hessenberg form ====
500 *
501  IF( ns.EQ.0 )
502  $ s = zero
503 *
504  IF( ns.LT.jw ) THEN
505 *
506 * ==== sorting diagonal blocks of T improves accuracy for
507 * . graded matrices. Bubble sort deals well with
508 * . exchange failures. ====
509 *
510  sorted = .false.
511  i = ns + 1
512  30 CONTINUE
513  IF( sorted )
514  $ GO TO 50
515  sorted = .true.
516 *
517  kend = i - 1
518  i = infqr + 1
519  IF( i.EQ.ns ) THEN
520  k = i + 1
521  ELSE IF( t( i+1, i ).EQ.zero ) THEN
522  k = i + 1
523  ELSE
524  k = i + 2
525  END IF
526  40 CONTINUE
527  IF( k.LE.kend ) THEN
528  IF( k.EQ.i+1 ) THEN
529  evi = abs( t( i, i ) )
530  ELSE
531  evi = abs( t( i, i ) ) + sqrt( abs( t( i+1, i ) ) )*
532  $ sqrt( abs( t( i, i+1 ) ) )
533  END IF
534 *
535  IF( k.EQ.kend ) THEN
536  evk = abs( t( k, k ) )
537  ELSE IF( t( k+1, k ).EQ.zero ) THEN
538  evk = abs( t( k, k ) )
539  ELSE
540  evk = abs( t( k, k ) ) + sqrt( abs( t( k+1, k ) ) )*
541  $ sqrt( abs( t( k, k+1 ) ) )
542  END IF
543 *
544  IF( evi.GE.evk ) THEN
545  i = k
546  ELSE
547  sorted = .false.
548  ifst = i
549  ilst = k
550  CALL dtrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
551  $ info )
552  IF( info.EQ.0 ) THEN
553  i = ilst
554  ELSE
555  i = k
556  END IF
557  END IF
558  IF( i.EQ.kend ) THEN
559  k = i + 1
560  ELSE IF( t( i+1, i ).EQ.zero ) THEN
561  k = i + 1
562  ELSE
563  k = i + 2
564  END IF
565  GO TO 40
566  END IF
567  GO TO 30
568  50 CONTINUE
569  END IF
570 *
571 * ==== Restore shift/eigenvalue array from T ====
572 *
573  i = jw
574  60 CONTINUE
575  IF( i.GE.infqr+1 ) THEN
576  IF( i.EQ.infqr+1 ) THEN
577  sr( kwtop+i-1 ) = t( i, i )
578  si( kwtop+i-1 ) = zero
579  i = i - 1
580  ELSE IF( t( i, i-1 ).EQ.zero ) THEN
581  sr( kwtop+i-1 ) = t( i, i )
582  si( kwtop+i-1 ) = zero
583  i = i - 1
584  ELSE
585  aa = t( i-1, i-1 )
586  cc = t( i, i-1 )
587  bb = t( i-1, i )
588  dd = t( i, i )
589  CALL dlanv2( aa, bb, cc, dd, sr( kwtop+i-2 ),
590  $ si( kwtop+i-2 ), sr( kwtop+i-1 ),
591  $ si( kwtop+i-1 ), cs, sn )
592  i = i - 2
593  END IF
594  GO TO 60
595  END IF
596 *
597  IF( ns.LT.jw .OR. s.EQ.zero ) THEN
598  IF( ns.GT.1 .AND. s.NE.zero ) THEN
599 *
600 * ==== Reflect spike back into lower triangle ====
601 *
602  CALL dcopy( ns, v, ldv, work, 1 )
603  beta = work( 1 )
604  CALL dlarfg( ns, beta, work( 2 ), 1, tau )
605  work( 1 ) = one
606 *
607  CALL dlaset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
608 *
609  CALL dlarf( 'L', ns, jw, work, 1, tau, t, ldt,
610  $ work( jw+1 ) )
611  CALL dlarf( 'R', ns, ns, work, 1, tau, t, ldt,
612  $ work( jw+1 ) )
613  CALL dlarf( 'R', jw, ns, work, 1, tau, v, ldv,
614  $ work( jw+1 ) )
615 *
616  CALL dgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
617  $ lwork-jw, info )
618  END IF
619 *
620 * ==== Copy updated reduced window into place ====
621 *
622  IF( kwtop.GT.1 )
623  $ h( kwtop, kwtop-1 ) = s*v( 1, 1 )
624  CALL dlacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
625  CALL dcopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
626  $ ldh+1 )
627 *
628 * ==== Accumulate orthogonal matrix in order update
629 * . H and Z, if requested. ====
630 *
631  IF( ns.GT.1 .AND. s.NE.zero )
632  $ CALL dormhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
633  $ work( jw+1 ), lwork-jw, info )
634 *
635 * ==== Update vertical slab in H ====
636 *
637  IF( wantt ) THEN
638  ltop = 1
639  ELSE
640  ltop = ktop
641  END IF
642  DO 70 krow = ltop, kwtop - 1, nv
643  kln = min( nv, kwtop-krow )
644  CALL dgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
645  $ ldh, v, ldv, zero, wv, ldwv )
646  CALL dlacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
647  70 CONTINUE
648 *
649 * ==== Update horizontal slab in H ====
650 *
651  IF( wantt ) THEN
652  DO 80 kcol = kbot + 1, n, nh
653  kln = min( nh, n-kcol+1 )
654  CALL dgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
655  $ h( kwtop, kcol ), ldh, zero, t, ldt )
656  CALL dlacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
657  $ ldh )
658  80 CONTINUE
659  END IF
660 *
661 * ==== Update vertical slab in Z ====
662 *
663  IF( wantz ) THEN
664  DO 90 krow = iloz, ihiz, nv
665  kln = min( nv, ihiz-krow+1 )
666  CALL dgemm( 'N', 'N', kln, jw, jw, one, z( krow, kwtop ),
667  $ ldz, v, ldv, zero, wv, ldwv )
668  CALL dlacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
669  $ ldz )
670  90 CONTINUE
671  END IF
672  END IF
673 *
674 * ==== Return the number of deflations ... ====
675 *
676  nd = jw - ns
677 *
678 * ==== ... and the number of shifts. (Subtracting
679 * . INFQR from the spike length takes care
680 * . of the case of a rare QR failure while
681 * . calculating eigenvalues of the deflation
682 * . window.) ====
683 *
684  ns = ns - infqr
685 *
686 * ==== Return optimal workspace. ====
687 *
688  work( 1 ) = dble( lwkopt )
689 *
690 * ==== End of DLAQR3 ====
691 *
692  END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:167
subroutine dlanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition: dlanv2.f:127
subroutine dlaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
DLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: dlaqr4.f:263
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:124
subroutine dlaqr3(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
DLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate...
Definition: dlaqr3.f:275
subroutine dlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition: dlahqr.f:207
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:106
subroutine dtrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
DTREXC
Definition: dtrexc.f:148
subroutine dormhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMHR
Definition: dormhr.f:178