← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to src/lapack/double/zlaqr4.f

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lapack/double/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