Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / zlaqr4.f
1 *> \brief \b ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition.
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 *            http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZLAQR4 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqr4.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqr4.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqr4.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE ZLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
22 *                          IHIZ, Z, LDZ, WORK, LWORK, INFO )
23 *
24 *       .. Scalar Arguments ..
25 *       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
26 *       LOGICAL            WANTT, WANTZ
27 *       ..
28 *       .. Array Arguments ..
29 *       COMPLEX*16         H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
30 *       ..
31 *
32 *
33 *> \par Purpose:
34 *  =============
35 *>
36 *> \verbatim
37 *>
38 *>    ZLAQR4 implements one level of recursion for ZLAQR0.
39 *>    It is a complete implementation of the small bulge multi-shift
40 *>    QR algorithm.  It may be called by ZLAQR0 and, for large enough
41 *>    deflation window size, it may be called by ZLAQR3.  This
42 *>    subroutine is identical to ZLAQR0 except that it calls ZLAQR2
43 *>    instead of ZLAQR3.
44 *>
45 *>    ZLAQR4 computes the eigenvalues of a Hessenberg matrix H
46 *>    and, optionally, the matrices T and Z from the Schur decomposition
47 *>    H = Z T Z**H, where T is an upper triangular matrix (the
48 *>    Schur form), and Z is the unitary matrix of Schur vectors.
49 *>
50 *>    Optionally Z may be postmultiplied into an input unitary
51 *>    matrix Q so that this routine can give the Schur factorization
52 *>    of a matrix A which has been reduced to the Hessenberg form H
53 *>    by the unitary matrix Q:  A = Q*H*Q**H = (QZ)*H*(QZ)**H.
54 *> \endverbatim
55 *
56 *  Arguments:
57 *  ==========
58 *
59 *> \param[in] WANTT
60 *> \verbatim
61 *>          WANTT is LOGICAL
62 *>          = .TRUE. : the full Schur form T is required;
63 *>          = .FALSE.: only eigenvalues are required.
64 *> \endverbatim
65 *>
66 *> \param[in] WANTZ
67 *> \verbatim
68 *>          WANTZ is LOGICAL
69 *>          = .TRUE. : the matrix of Schur vectors Z is required;
70 *>          = .FALSE.: Schur vectors are not required.
71 *> \endverbatim
72 *>
73 *> \param[in] N
74 *> \verbatim
75 *>          N is INTEGER
76 *>           The order of the matrix H.  N .GE. 0.
77 *> \endverbatim
78 *>
79 *> \param[in] ILO
80 *> \verbatim
81 *>          ILO is INTEGER
82 *> \endverbatim
83 *>
84 *> \param[in] IHI
85 *> \verbatim
86 *>          IHI is INTEGER
87 *>           It is assumed that H is already upper triangular in rows
88 *>           and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
89 *>           H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
90 *>           previous call to ZGEBAL, and then passed to ZGEHRD when the
91 *>           matrix output by ZGEBAL is reduced to Hessenberg form.
92 *>           Otherwise, ILO and IHI should be set to 1 and N,
93 *>           respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
94 *>           If N = 0, then ILO = 1 and IHI = 0.
95 *> \endverbatim
96 *>
97 *> \param[in,out] H
98 *> \verbatim
99 *>          H is COMPLEX*16 array, dimension (LDH,N)
100 *>           On entry, the upper Hessenberg matrix H.
101 *>           On exit, if INFO = 0 and WANTT is .TRUE., then H
102 *>           contains the upper triangular matrix T from the Schur
103 *>           decomposition (the Schur form). If INFO = 0 and WANT is
104 *>           .FALSE., then the contents of H are unspecified on exit.
105 *>           (The output value of H when INFO.GT.0 is given under the
106 *>           description of INFO below.)
107 *>
108 *>           This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
109 *>           j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
110 *> \endverbatim
111 *>
112 *> \param[in] LDH
113 *> \verbatim
114 *>          LDH is INTEGER
115 *>           The leading dimension of the array H. LDH .GE. max(1,N).
116 *> \endverbatim
117 *>
118 *> \param[out] W
119 *> \verbatim
120 *>          W is COMPLEX*16 array, dimension (N)
121 *>           The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored
122 *>           in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are
123 *>           stored in the same order as on the diagonal of the Schur
124 *>           form returned in H, with W(i) = H(i,i).
125 *> \endverbatim
126 *>
127 *> \param[in] ILOZ
128 *> \verbatim
129 *>          ILOZ is INTEGER
130 *> \endverbatim
131 *>
132 *> \param[in] IHIZ
133 *> \verbatim
134 *>          IHIZ is INTEGER
135 *>           Specify the rows of Z to which transformations must be
136 *>           applied if WANTZ is .TRUE..
137 *>           1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
138 *> \endverbatim
139 *>
140 *> \param[in,out] Z
141 *> \verbatim
142 *>          Z is COMPLEX*16 array, dimension (LDZ,IHI)
143 *>           If WANTZ is .FALSE., then Z is not referenced.
144 *>           If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
145 *>           replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
146 *>           orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
147 *>           (The output value of Z when INFO.GT.0 is given under
148 *>           the description of INFO below.)
149 *> \endverbatim
150 *>
151 *> \param[in] LDZ
152 *> \verbatim
153 *>          LDZ is INTEGER
154 *>           The leading dimension of the array Z.  if WANTZ is .TRUE.
155 *>           then LDZ.GE.MAX(1,IHIZ).  Otherwize, LDZ.GE.1.
156 *> \endverbatim
157 *>
158 *> \param[out] WORK
159 *> \verbatim
160 *>          WORK is COMPLEX*16 array, dimension LWORK
161 *>           On exit, if LWORK = -1, WORK(1) returns an estimate of
162 *>           the optimal value for LWORK.
163 *> \endverbatim
164 *>
165 *> \param[in] LWORK
166 *> \verbatim
167 *>          LWORK is INTEGER
168 *>           The dimension of the array WORK.  LWORK .GE. max(1,N)
169 *>           is sufficient, but LWORK typically as large as 6*N may
170 *>           be required for optimal performance.  A workspace query
171 *>           to determine the optimal workspace size is recommended.
172 *>
173 *>           If LWORK = -1, then ZLAQR4 does a workspace query.
174 *>           In this case, ZLAQR4 checks the input parameters and
175 *>           estimates the optimal workspace size for the given
176 *>           values of N, ILO and IHI.  The estimate is returned
177 *>           in WORK(1).  No error message related to LWORK is
178 *>           issued by XERBLA.  Neither H nor Z are accessed.
179 *> \endverbatim
180 *>
181 *> \param[out] INFO
182 *> \verbatim
183 *>          INFO is INTEGER
184 *>             =  0:  successful exit
185 *>           .GT. 0:  if INFO = i, ZLAQR4 failed to compute all of
186 *>                the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
187 *>                and WI contain those eigenvalues which have been
188 *>                successfully computed.  (Failures are rare.)
189 *>
190 *>                If INFO .GT. 0 and WANT is .FALSE., then on exit,
191 *>                the remaining unconverged eigenvalues are the eigen-
192 *>                values of the upper Hessenberg matrix rows and
193 *>                columns ILO through INFO of the final, output
194 *>                value of H.
195 *>
196 *>                If INFO .GT. 0 and WANTT is .TRUE., then on exit
197 *>
198 *>           (*)  (initial value of H)*U  = U*(final value of H)
199 *>
200 *>                where U is a unitary matrix.  The final
201 *>                value of  H is upper Hessenberg and triangular in
202 *>                rows and columns INFO+1 through IHI.
203 *>
204 *>                If INFO .GT. 0 and WANTZ is .TRUE., then on exit
205 *>
206 *>                  (final value of Z(ILO:IHI,ILOZ:IHIZ)
207 *>                   =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
208 *>
209 *>                where U is the unitary matrix in (*) (regard-
210 *>                less of the value of WANTT.)
211 *>
212 *>                If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
213 *>                accessed.
214 *> \endverbatim
215 *
216 *  Authors:
217 *  ========
218 *
219 *> \author Univ. of Tennessee
220 *> \author Univ. of California Berkeley
221 *> \author Univ. of Colorado Denver
222 *> \author NAG Ltd.
223 *
224 *> \date September 2012
225 *
226 *> \ingroup complex16OTHERauxiliary
227 *
228 *> \par Contributors:
229 *  ==================
230 *>
231 *>       Karen Braman and Ralph Byers, Department of Mathematics,
232 *>       University of Kansas, USA
233 *
234 *> \par References:
235 *  ================
236 *>
237 *>       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
238 *>       Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
239 *>       Performance, SIAM Journal of Matrix Analysis, volume 23, pages
240 *>       929--947, 2002.
241 *> \n
242 *>       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
243 *>       Algorithm Part II: Aggressive Early Deflation, SIAM Journal
244 *>       of Matrix Analysis, volume 23, pages 948--973, 2002.
245 *>
246 *  =====================================================================
247       SUBROUTINE ZLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
248      $                   IHIZ, Z, LDZ, WORK, LWORK, INFO )
249 *
250 *  -- LAPACK auxiliary routine (version 3.4.2) --
251 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
252 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253 *     September 2012
254 *
255 *     .. Scalar Arguments ..
256       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
257       LOGICAL            WANTT, WANTZ
258 *     ..
259 *     .. Array Arguments ..
260       COMPLEX*16         H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
261 *     ..
262 *
263 *  ================================================================
264 *
265 *     .. Parameters ..
266 *
267 *     ==== Matrices of order NTINY or smaller must be processed by
268 *     .    ZLAHQR because of insufficient subdiagonal scratch space.
269 *     .    (This is a hard limit.) ====
270       INTEGER            NTINY
271       PARAMETER          ( NTINY = 11 )
272 *
273 *     ==== Exceptional deflation windows:  try to cure rare
274 *     .    slow convergence by varying the size of the
275 *     .    deflation window after KEXNW iterations. ====
276       INTEGER            KEXNW
277       PARAMETER          ( KEXNW = 5 )
278 *
279 *     ==== Exceptional shifts: try to cure rare slow convergence
280 *     .    with ad-hoc exceptional shifts every KEXSH iterations.
281 *     .    ====
282       INTEGER            KEXSH
283       PARAMETER          ( KEXSH = 6 )
284 *
285 *     ==== The constant WILK1 is used to form the exceptional
286 *     .    shifts. ====
287       DOUBLE PRECISION   WILK1
288       PARAMETER          ( WILK1 = 0.75d0 )
289       COMPLEX*16         ZERO, ONE
290       PARAMETER          ( ZERO = ( 0.0d0, 0.0d0 ),
291      $                   ONE = ( 1.0d0, 0.0d0 ) )
292       DOUBLE PRECISION   TWO
293       PARAMETER          ( TWO = 2.0d0 )
294 *     ..
295 *     .. Local Scalars ..
296       COMPLEX*16         AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
297       DOUBLE PRECISION   S
298       INTEGER            I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
299      $                   KT, KTOP, KU, KV, KWH, KWTOP, KWV, LD, LS,
300      $                   LWKOPT, NDEC, NDFL, NH, NHO, NIBBLE, NMIN, NS,
301      $                   NSMAX, NSR, NVE, NW, NWMAX, NWR, NWUPBD
302       LOGICAL            SORTED
303       CHARACTER          JBCMPZ*2
304 *     ..
305 *     .. External Functions ..
306       INTEGER            ILAENV
307       EXTERNAL           ILAENV
308 *     ..
309 *     .. Local Arrays ..
310       COMPLEX*16         ZDUM( 1, 1 )
311 *     ..
312 *     .. External Subroutines ..
313       EXTERNAL           ZLACPY, ZLAHQR, ZLAQR2, ZLAQR5
314 *     ..
315 *     .. Intrinsic Functions ..
316       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, INT, MAX, MIN, MOD,
317      $                   SQRT
318 *     ..
319 *     .. Statement Functions ..
320       DOUBLE PRECISION   CABS1
321 *     ..
322 *     .. Statement Function definitions ..
323       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
324 *     ..
325 *     .. Executable Statements ..
326       INFO = 0
327 *
328 *     ==== Quick return for N = 0: nothing to do. ====
329 *
330       IF( N.EQ.0 ) THEN
331          WORK( 1 ) = ONE
332          RETURN
333       END IF
334 *
335       IF( N.LE.NTINY ) THEN
336 *
337 *        ==== Tiny matrices must use ZLAHQR. ====
338 *
339          LWKOPT = 1
340          IF( LWORK.NE.-1 )
341      $      CALL ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
342      $                   IHIZ, Z, LDZ, INFO )
343       ELSE
344 *
345 *        ==== Use small bulge multi-shift QR with aggressive early
346 *        .    deflation on larger-than-tiny matrices. ====
347 *
348 *        ==== Hope for the best. ====
349 *
350          INFO = 0
351 *
352 *        ==== Set up job flags for ILAENV. ====
353 *
354          IF( WANTT ) THEN
355             JBCMPZ( 1: 1 ) = 'S'
356          ELSE
357             JBCMPZ( 1: 1 ) = 'E'
358          END IF
359          IF( WANTZ ) THEN
360             JBCMPZ( 2: 2 ) = 'V'
361          ELSE
362             JBCMPZ( 2: 2 ) = 'N'
363          END IF
364 *
365 *        ==== NWR = recommended deflation window size.  At this
366 *        .    point,  N .GT. NTINY = 11, so there is enough
367 *        .    subdiagonal workspace for NWR.GE.2 as required.
368 *        .    (In fact, there is enough subdiagonal space for
369 *        .    NWR.GE.3.) ====
370 *
371          NWR = ILAENV( 13, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK )
372          NWR = MAX( 2, NWR )
373          NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR )
374 *
375 *        ==== NSR = recommended number of simultaneous shifts.
376 *        .    At this point N .GT. NTINY = 11, so there is at
377 *        .    enough subdiagonal workspace for NSR to be even
378 *        .    and greater than or equal to two as required. ====
379 *
380          NSR = ILAENV( 15, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK )
381          NSR = MIN( NSR, ( N+6 ) / 9, IHI-ILO )
382          NSR = MAX( 2, NSR-MOD( NSR, 2 ) )
383 *
384 *        ==== Estimate optimal workspace ====
385 *
386 *        ==== Workspace query call to ZLAQR2 ====
387 *
388          CALL ZLAQR2( WANTT, WANTZ, N, ILO, IHI, NWR+1, H, LDH, ILOZ,
389      $                IHIZ, Z, LDZ, LS, LD, W, H, LDH, N, H, LDH, N, H,
390      $                LDH, WORK, -1 )
391 *
392 *        ==== Optimal workspace = MAX(ZLAQR5, ZLAQR2) ====
393 *
394          LWKOPT = MAX( 3*NSR / 2, INT( WORK( 1 ) ) )
395 *
396 *        ==== Quick return in case of workspace query. ====
397 *
398          IF( LWORK.EQ.-1 ) THEN
399             WORK( 1 ) = DCMPLX( LWKOPT, 0 )
400             RETURN
401          END IF
402 *
403 *        ==== ZLAHQR/ZLAQR0 crossover point ====
404 *
405          NMIN = ILAENV( 12, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK )
406          NMIN = MAX( NTINY, NMIN )
407 *
408 *        ==== Nibble crossover point ====
409 *
410          NIBBLE = ILAENV( 14, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK )
411          NIBBLE = MAX( 0, NIBBLE )
412 *
413 *        ==== Accumulate reflections during ttswp?  Use block
414 *        .    2-by-2 structure during matrix-matrix multiply? ====
415 *
416          KACC22 = ILAENV( 16, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK )
417          KACC22 = MAX( 0, KACC22 )
418          KACC22 = MIN( 2, KACC22 )
419 *
420 *        ==== NWMAX = the largest possible deflation window for
421 *        .    which there is sufficient workspace. ====
422 *
423          NWMAX = MIN( ( N-1 ) / 3, LWORK / 2 )
424          NW = NWMAX
425 *
426 *        ==== NSMAX = the Largest number of simultaneous shifts
427 *        .    for which there is sufficient workspace. ====
428 *
429          NSMAX = MIN( ( N+6 ) / 9, 2*LWORK / 3 )
430          NSMAX = NSMAX - MOD( NSMAX, 2 )
431 *
432 *        ==== NDFL: an iteration count restarted at deflation. ====
433 *
434          NDFL = 1
435 *
436 *        ==== ITMAX = iteration limit ====
437 *
438          ITMAX = MAX( 30, 2*KEXSH )*MAX( 10, ( IHI-ILO+1 ) )
439 *
440 *        ==== Last row and column in the active block ====
441 *
442          KBOT = IHI
443 *
444 *        ==== Main Loop ====
445 *
446          DO 70 IT = 1, ITMAX
447 *
448 *           ==== Done when KBOT falls below ILO ====
449 *
450             IF( KBOT.LT.ILO )
451      $         GO TO 80
452 *
453 *           ==== Locate active block ====
454 *
455             DO 10 K = KBOT, ILO + 1, -1
456                IF( H( K, K-1 ).EQ.ZERO )
457      $            GO TO 20
458    10       CONTINUE
459             K = ILO
460    20       CONTINUE
461             KTOP = K
462 *
463 *           ==== Select deflation window size:
464 *           .    Typical Case:
465 *           .      If possible and advisable, nibble the entire
466 *           .      active block.  If not, use size MIN(NWR,NWMAX)
467 *           .      or MIN(NWR+1,NWMAX) depending upon which has
468 *           .      the smaller corresponding subdiagonal entry
469 *           .      (a heuristic).
470 *           .
471 *           .    Exceptional Case:
472 *           .      If there have been no deflations in KEXNW or
473 *           .      more iterations, then vary the deflation window
474 *           .      size.   At first, because, larger windows are,
475 *           .      in general, more powerful than smaller ones,
476 *           .      rapidly increase the window to the maximum possible.
477 *           .      Then, gradually reduce the window size. ====
478 *
479             NH = KBOT - KTOP + 1
480             NWUPBD = MIN( NH, NWMAX )
481             IF( NDFL.LT.KEXNW ) THEN
482                NW = MIN( NWUPBD, NWR )
483             ELSE
484                NW = MIN( NWUPBD, 2*NW )
485             END IF
486             IF( NW.LT.NWMAX ) THEN
487                IF( NW.GE.NH-1 ) THEN
488                   NW = NH
489                ELSE
490                   KWTOP = KBOT - NW + 1
491                   IF( CABS1( H( KWTOP, KWTOP-1 ) ).GT.
492      $                CABS1( H( KWTOP-1, KWTOP-2 ) ) )NW = NW + 1
493                END IF
494             END IF
495             IF( NDFL.LT.KEXNW ) THEN
496                NDEC = -1
497             ELSE IF( NDEC.GE.0 .OR. NW.GE.NWUPBD ) THEN
498                NDEC = NDEC + 1
499                IF( NW-NDEC.LT.2 )
500      $            NDEC = 0
501                NW = NW - NDEC
502             END IF
503 *
504 *           ==== Aggressive early deflation:
505 *           .    split workspace under the subdiagonal into
506 *           .      - an nw-by-nw work array V in the lower
507 *           .        left-hand-corner,
508 *           .      - an NW-by-at-least-NW-but-more-is-better
509 *           .        (NW-by-NHO) horizontal work array along
510 *           .        the bottom edge,
511 *           .      - an at-least-NW-but-more-is-better (NHV-by-NW)
512 *           .        vertical work array along the left-hand-edge.
513 *           .        ====
514 *
515             KV = N - NW + 1
516             KT = NW + 1
517             NHO = ( N-NW-1 ) - KT + 1
518             KWV = NW + 2
519             NVE = ( N-NW ) - KWV + 1
520 *
521 *           ==== Aggressive early deflation ====
522 *
523             CALL ZLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
524      $                   IHIZ, Z, LDZ, LS, LD, W, H( KV, 1 ), LDH, NHO,
525      $                   H( KV, KT ), LDH, NVE, H( KWV, 1 ), LDH, WORK,
526      $                   LWORK )
527 *
528 *           ==== Adjust KBOT accounting for new deflations. ====
529 *
530             KBOT = KBOT - LD
531 *
532 *           ==== KS points to the shifts. ====
533 *
534             KS = KBOT - LS + 1
535 *
536 *           ==== Skip an expensive QR sweep if there is a (partly
537 *           .    heuristic) reason to expect that many eigenvalues
538 *           .    will deflate without it.  Here, the QR sweep is
539 *           .    skipped if many eigenvalues have just been deflated
540 *           .    or if the remaining active block is small.
541 *
542             IF( ( LD.EQ.0 ) .OR. ( ( 100*LD.LE.NW*NIBBLE ) .AND. ( KBOT-
543      $          KTOP+1.GT.MIN( NMIN, NWMAX ) ) ) ) THEN
544 *
545 *              ==== NS = nominal number of simultaneous shifts.
546 *              .    This may be lowered (slightly) if ZLAQR2
547 *              .    did not provide that many shifts. ====
548 *
549                NS = MIN( NSMAX, NSR, MAX( 2, KBOT-KTOP ) )
550                NS = NS - MOD( NS, 2 )
551 *
552 *              ==== If there have been no deflations
553 *              .    in a multiple of KEXSH iterations,
554 *              .    then try exceptional shifts.
555 *              .    Otherwise use shifts provided by
556 *              .    ZLAQR2 above or from the eigenvalues
557 *              .    of a trailing principal submatrix. ====
558 *
559                IF( MOD( NDFL, KEXSH ).EQ.0 ) THEN
560                   KS = KBOT - NS + 1
561                   DO 30 I = KBOT, KS + 1, -2
562                      W( I ) = H( I, I ) + WILK1*CABS1( H( I, I-1 ) )
563                      W( I-1 ) = W( I )
564    30             CONTINUE
565                ELSE
566 *
567 *                 ==== Got NS/2 or fewer shifts? Use ZLAHQR
568 *                 .    on a trailing principal submatrix to
569 *                 .    get more. (Since NS.LE.NSMAX.LE.(N+6)/9,
570 *                 .    there is enough space below the subdiagonal
571 *                 .    to fit an NS-by-NS scratch array.) ====
572 *
573                   IF( KBOT-KS+1.LE.NS / 2 ) THEN
574                      KS = KBOT - NS + 1
575                      KT = N - NS + 1
576                      CALL ZLACPY( 'A', NS, NS, H( KS, KS ), LDH,
577      $                            H( KT, 1 ), LDH )
578                      CALL ZLAHQR( .false., .false., NS, 1, NS,
579      $                            H( KT, 1 ), LDH, W( KS ), 1, 1, ZDUM,
580      $                            1, INF )
581                      KS = KS + INF
582 *
583 *                    ==== In case of a rare QR failure use
584 *                    .    eigenvalues of the trailing 2-by-2
585 *                    .    principal submatrix.  Scale to avoid
586 *                    .    overflows, underflows and subnormals.
587 *                    .    (The scale factor S can not be zero,
588 *                    .    because H(KBOT,KBOT-1) is nonzero.) ====
589 *
590                      IF( KS.GE.KBOT ) THEN
591                         S = CABS1( H( KBOT-1, KBOT-1 ) ) +
592      $                      CABS1( H( KBOT, KBOT-1 ) ) +
593      $                      CABS1( H( KBOT-1, KBOT ) ) +
594      $                      CABS1( H( KBOT, KBOT ) )
595                         AA = H( KBOT-1, KBOT-1 ) / S
596                         CC = H( KBOT, KBOT-1 ) / S
597                         BB = H( KBOT-1, KBOT ) / S
598                         DD = H( KBOT, KBOT ) / S
599                         TR2 = ( AA+DD ) / TWO
600                         DET = ( AA-TR2 )*( DD-TR2 ) - BB*CC
601                         RTDISC = SQRT( -DET )
602                         W( KBOT-1 ) = ( TR2+RTDISC )*S
603                         W( KBOT ) = ( TR2-RTDISC )*S
604 *
605                         KS = KBOT - 1
606                      END IF
607                   END IF
608 *
609                   IF( KBOT-KS+1.GT.NS ) THEN
610 *
611 *                    ==== Sort the shifts (Helps a little) ====
612 *
613                      SORTED = .false.
614                      DO 50 K = KBOT, KS + 1, -1
615                         IF( SORTED )
616      $                     GO TO 60
617                         SORTED = .true.
618                         DO 40 I = KS, K - 1
619                            IF( CABS1( W( I ) ).LT.CABS1( W( I+1 ) ) )
620      $                          THEN
621                               SORTED = .false.
622                               SWAP = W( I )
623                               W( I ) = W( I+1 )
624                               W( I+1 ) = SWAP
625                            END IF
626    40                   CONTINUE
627    50                CONTINUE
628    60                CONTINUE
629                   END IF
630                END IF
631 *
632 *              ==== If there are only two shifts, then use
633 *              .    only one.  ====
634 *
635                IF( KBOT-KS+1.EQ.2 ) THEN
636                   IF( CABS1( W( KBOT )-H( KBOT, KBOT ) ).LT.
637      $                CABS1( W( KBOT-1 )-H( KBOT, KBOT ) ) ) THEN
638                      W( KBOT-1 ) = W( KBOT )
639                   ELSE
640                      W( KBOT ) = W( KBOT-1 )
641                   END IF
642                END IF
643 *
644 *              ==== Use up to NS of the the smallest magnatiude
645 *              .    shifts.  If there aren't NS shifts available,
646 *              .    then use them all, possibly dropping one to
647 *              .    make the number of shifts even. ====
648 *
649                NS = MIN( NS, KBOT-KS+1 )
650                NS = NS - MOD( NS, 2 )
651                KS = KBOT - NS + 1
652 *
653 *              ==== Small-bulge multi-shift QR sweep:
654 *              .    split workspace under the subdiagonal into
655 *              .    - a KDU-by-KDU work array U in the lower
656 *              .      left-hand-corner,
657 *              .    - a KDU-by-at-least-KDU-but-more-is-better
658 *              .      (KDU-by-NHo) horizontal work array WH along
659 *              .      the bottom edge,
660 *              .    - and an at-least-KDU-but-more-is-better-by-KDU
661 *              .      (NVE-by-KDU) vertical work WV arrow along
662 *              .      the left-hand-edge. ====
663 *
664                KDU = 3*NS - 3
665                KU = N - KDU + 1
666                KWH = KDU + 1
667                NHO = ( N-KDU+1-4 ) - ( KDU+1 ) + 1
668                KWV = KDU + 4
669                NVE = N - KDU - KWV + 1
670 *
671 *              ==== Small-bulge multi-shift QR sweep ====
672 *
673                CALL ZLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NS,
674      $                      W( KS ), H, LDH, ILOZ, IHIZ, Z, LDZ, WORK,
675      $                      3, H( KU, 1 ), LDH, NVE, H( KWV, 1 ), LDH,
676      $                      NHO, H( KU, KWH ), LDH )
677             END IF
678 *
679 *           ==== Note progress (or the lack of it). ====
680 *
681             IF( LD.GT.0 ) THEN
682                NDFL = 1
683             ELSE
684                NDFL = NDFL + 1
685             END IF
686 *
687 *           ==== End of main loop ====
688    70    CONTINUE
689 *
690 *        ==== Iteration limit exceeded.  Set INFO to show where
691 *        .    the problem occurred and exit. ====
692 *
693          INFO = KBOT
694    80    CONTINUE
695       END IF
696 *
697 *     ==== Return the optimal value of LWORK. ====
698 *
699       WORK( 1 ) = DCMPLX( LWKOPT, 0 )
700 *
701 *     ==== End of ZLAQR4 ====
702 *
703       END