← Back to branch summary

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

« back to all changes in this revision

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