← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/lapack/double/zhseqr.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/zhseqr.f
 
1
*> \brief \b ZHSEQR
 
2
*
 
3
*  =========== DOCUMENTATION ===========
 
4
*
 
5
* Online html documentation available at 
 
6
*            http://www.netlib.org/lapack/explore-html/ 
 
7
*
 
8
*> \htmlonly
 
9
*> Download ZHSEQR + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhseqr.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhseqr.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhseqr.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE ZHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ,
 
22
*                          WORK, LWORK, INFO )
 
23
* 
 
24
*       .. Scalar Arguments ..
 
25
*       INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
 
26
*       CHARACTER          COMPZ, JOB
 
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
*>    ZHSEQR computes the eigenvalues of a Hessenberg matrix H
 
39
*>    and, optionally, the matrices T and Z from the Schur decomposition
 
40
*>    H = Z T Z**H, where T is an upper triangular matrix (the
 
41
*>    Schur form), and Z is the unitary matrix of Schur vectors.
 
42
*>
 
43
*>    Optionally Z may be postmultiplied into an input unitary
 
44
*>    matrix Q so that this routine can give the Schur factorization
 
45
*>    of a matrix A which has been reduced to the Hessenberg form H
 
46
*>    by the unitary matrix Q:  A = Q*H*Q**H = (QZ)*H*(QZ)**H.
 
47
*> \endverbatim
 
48
*
 
49
*  Arguments:
 
50
*  ==========
 
51
*
 
52
*> \param[in] JOB
 
53
*> \verbatim
 
54
*>          JOB is CHARACTER*1
 
55
*>           = 'E':  compute eigenvalues only;
 
56
*>           = 'S':  compute eigenvalues and the Schur form T.
 
57
*> \endverbatim
 
58
*>
 
59
*> \param[in] COMPZ
 
60
*> \verbatim
 
61
*>          COMPZ is CHARACTER*1
 
62
*>           = 'N':  no Schur vectors are computed;
 
63
*>           = 'I':  Z is initialized to the unit matrix and the matrix Z
 
64
*>                   of Schur vectors of H is returned;
 
65
*>           = 'V':  Z must contain an unitary matrix Q on entry, and
 
66
*>                   the product Q*Z is returned.
 
67
*> \endverbatim
 
68
*>
 
69
*> \param[in] N
 
70
*> \verbatim
 
71
*>          N is INTEGER
 
72
*>           The order of the matrix H.  N .GE. 0.
 
73
*> \endverbatim
 
74
*>
 
75
*> \param[in] ILO
 
76
*> \verbatim
 
77
*>          ILO is INTEGER
 
78
*> \endverbatim
 
79
*>
 
80
*> \param[in] IHI
 
81
*> \verbatim
 
82
*>          IHI is INTEGER
 
83
*>
 
84
*>           It is assumed that H is already upper triangular in rows
 
85
*>           and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
 
86
*>           set by a previous call to ZGEBAL, and then passed to ZGEHRD
 
87
*>           when the matrix output by ZGEBAL is reduced to Hessenberg
 
88
*>           form. Otherwise ILO and IHI should be set to 1 and N
 
89
*>           respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
 
90
*>           If N = 0, then ILO = 1 and IHI = 0.
 
91
*> \endverbatim
 
92
*>
 
93
*> \param[in,out] H
 
94
*> \verbatim
 
95
*>          H is COMPLEX*16 array, dimension (LDH,N)
 
96
*>           On entry, the upper Hessenberg matrix H.
 
97
*>           On exit, if INFO = 0 and JOB = 'S', H contains the upper
 
98
*>           triangular matrix T from the Schur decomposition (the
 
99
*>           Schur form). If INFO = 0 and JOB = 'E', the contents of
 
100
*>           H are unspecified on exit.  (The output value of H when
 
101
*>           INFO.GT.0 is given under the description of INFO below.)
 
102
*>
 
103
*>           Unlike earlier versions of ZHSEQR, this subroutine may
 
104
*>           explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
 
105
*>           or j = IHI+1, IHI+2, ... N.
 
106
*> \endverbatim
 
107
*>
 
108
*> \param[in] LDH
 
109
*> \verbatim
 
110
*>          LDH is INTEGER
 
111
*>           The leading dimension of the array H. LDH .GE. max(1,N).
 
112
*> \endverbatim
 
113
*>
 
114
*> \param[out] W
 
115
*> \verbatim
 
116
*>          W is COMPLEX*16 array, dimension (N)
 
117
*>           The computed eigenvalues. If JOB = 'S', the eigenvalues are
 
118
*>           stored in the same order as on the diagonal of the Schur
 
119
*>           form returned in H, with W(i) = H(i,i).
 
120
*> \endverbatim
 
121
*>
 
122
*> \param[in,out] Z
 
123
*> \verbatim
 
124
*>          Z is COMPLEX*16 array, dimension (LDZ,N)
 
125
*>           If COMPZ = 'N', Z is not referenced.
 
126
*>           If COMPZ = 'I', on entry Z need not be set and on exit,
 
127
*>           if INFO = 0, Z contains the unitary matrix Z of the Schur
 
128
*>           vectors of H.  If COMPZ = 'V', on entry Z must contain an
 
129
*>           N-by-N matrix Q, which is assumed to be equal to the unit
 
130
*>           matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
 
131
*>           if INFO = 0, Z contains Q*Z.
 
132
*>           Normally Q is the unitary matrix generated by ZUNGHR
 
133
*>           after the call to ZGEHRD which formed the Hessenberg matrix
 
134
*>           H. (The output value of Z when INFO.GT.0 is given under
 
135
*>           the description of INFO below.)
 
136
*> \endverbatim
 
137
*>
 
138
*> \param[in] LDZ
 
139
*> \verbatim
 
140
*>          LDZ is INTEGER
 
141
*>           The leading dimension of the array Z.  if COMPZ = 'I' or
 
142
*>           COMPZ = 'V', then LDZ.GE.MAX(1,N).  Otherwize, LDZ.GE.1.
 
143
*> \endverbatim
 
144
*>
 
145
*> \param[out] WORK
 
146
*> \verbatim
 
147
*>          WORK is COMPLEX*16 array, dimension (LWORK)
 
148
*>           On exit, if INFO = 0, WORK(1) returns an estimate of
 
149
*>           the optimal value for LWORK.
 
150
*> \endverbatim
 
151
*>
 
152
*> \param[in] LWORK
 
153
*> \verbatim
 
154
*>          LWORK is INTEGER
 
155
*>           The dimension of the array WORK.  LWORK .GE. max(1,N)
 
156
*>           is sufficient and delivers very good and sometimes
 
157
*>           optimal performance.  However, LWORK as large as 11*N
 
158
*>           may be required for optimal performance.  A workspace
 
159
*>           query is recommended to determine the optimal workspace
 
160
*>           size.
 
161
*>
 
162
*>           If LWORK = -1, then ZHSEQR does a workspace query.
 
163
*>           In this case, ZHSEQR checks the input parameters and
 
164
*>           estimates the optimal workspace size for the given
 
165
*>           values of N, ILO and IHI.  The estimate is returned
 
166
*>           in WORK(1).  No error message related to LWORK is
 
167
*>           issued by XERBLA.  Neither H nor Z are accessed.
 
168
*> \endverbatim
 
169
*>
 
170
*> \param[out] INFO
 
171
*> \verbatim
 
172
*>          INFO is INTEGER
 
173
*>             =  0:  successful exit
 
174
*>           .LT. 0:  if INFO = -i, the i-th argument had an illegal
 
175
*>                    value
 
176
*>           .GT. 0:  if INFO = i, ZHSEQR failed to compute all of
 
177
*>                the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
 
178
*>                and WI contain those eigenvalues which have been
 
179
*>                successfully computed.  (Failures are rare.)
 
180
*>
 
181
*>                If INFO .GT. 0 and JOB = 'E', then on exit, the
 
182
*>                remaining unconverged eigenvalues are the eigen-
 
183
*>                values of the upper Hessenberg matrix rows and
 
184
*>                columns ILO through INFO of the final, output
 
185
*>                value of H.
 
186
*>
 
187
*>                If INFO .GT. 0 and JOB   = 'S', then on exit
 
188
*>
 
189
*>           (*)  (initial value of H)*U  = U*(final value of H)
 
190
*>
 
191
*>                where U is a unitary matrix.  The final
 
192
*>                value of  H is upper Hessenberg and triangular in
 
193
*>                rows and columns INFO+1 through IHI.
 
194
*>
 
195
*>                If INFO .GT. 0 and COMPZ = 'V', then on exit
 
196
*>
 
197
*>                  (final value of Z)  =  (initial value of Z)*U
 
198
*>
 
199
*>                where U is the unitary matrix in (*) (regard-
 
200
*>                less of the value of JOB.)
 
201
*>
 
202
*>                If INFO .GT. 0 and COMPZ = 'I', then on exit
 
203
*>                      (final value of Z)  = U
 
204
*>                where U is the unitary matrix in (*) (regard-
 
205
*>                less of the value of JOB.)
 
206
*>
 
207
*>                If INFO .GT. 0 and COMPZ = 'N', then Z is not
 
208
*>                accessed.
 
209
*> \endverbatim
 
210
*
 
211
*  Authors:
 
212
*  ========
 
213
*
 
214
*> \author Univ. of Tennessee 
 
215
*> \author Univ. of California Berkeley 
 
216
*> \author Univ. of Colorado Denver 
 
217
*> \author NAG Ltd. 
 
218
*
 
219
*> \date November 2011
 
220
*
 
221
*> \ingroup complex16OTHERcomputational
 
222
*
 
223
*> \par Contributors:
 
224
*  ==================
 
225
*>
 
226
*>       Karen Braman and Ralph Byers, Department of Mathematics,
 
227
*>       University of Kansas, USA
 
228
*
 
229
*> \par Further Details:
 
230
*  =====================
 
231
*>
 
232
*> \verbatim
 
233
*>
 
234
*>             Default values supplied by
 
235
*>             ILAENV(ISPEC,'ZHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
 
236
*>             It is suggested that these defaults be adjusted in order
 
237
*>             to attain best performance in each particular
 
238
*>             computational environment.
 
239
*>
 
240
*>            ISPEC=12: The ZLAHQR vs ZLAQR0 crossover point.
 
241
*>                      Default: 75. (Must be at least 11.)
 
242
*>
 
243
*>            ISPEC=13: Recommended deflation window size.
 
244
*>                      This depends on ILO, IHI and NS.  NS is the
 
245
*>                      number of simultaneous shifts returned
 
246
*>                      by ILAENV(ISPEC=15).  (See ISPEC=15 below.)
 
247
*>                      The default for (IHI-ILO+1).LE.500 is NS.
 
248
*>                      The default for (IHI-ILO+1).GT.500 is 3*NS/2.
 
249
*>
 
250
*>            ISPEC=14: Nibble crossover point. (See IPARMQ for
 
251
*>                      details.)  Default: 14% of deflation window
 
252
*>                      size.
 
253
*>
 
254
*>            ISPEC=15: Number of simultaneous shifts in a multishift
 
255
*>                      QR iteration.
 
256
*>
 
257
*>                      If IHI-ILO+1 is ...
 
258
*>
 
259
*>                      greater than      ...but less    ... the
 
260
*>                      or equal to ...      than        default is
 
261
*>
 
262
*>                           1               30          NS =   2(+)
 
263
*>                          30               60          NS =   4(+)
 
264
*>                          60              150          NS =  10(+)
 
265
*>                         150              590          NS =  **
 
266
*>                         590             3000          NS =  64
 
267
*>                        3000             6000          NS = 128
 
268
*>                        6000             infinity      NS = 256
 
269
*>
 
270
*>                  (+)  By default some or all matrices of this order
 
271
*>                       are passed to the implicit double shift routine
 
272
*>                       ZLAHQR and this parameter is ignored.  See
 
273
*>                       ISPEC=12 above and comments in IPARMQ for
 
274
*>                       details.
 
275
*>
 
276
*>                 (**)  The asterisks (**) indicate an ad-hoc
 
277
*>                       function of N increasing from 10 to 64.
 
278
*>
 
279
*>            ISPEC=16: Select structured matrix multiply.
 
280
*>                      If the number of simultaneous shifts (specified
 
281
*>                      by ISPEC=15) is less than 14, then the default
 
282
*>                      for ISPEC=16 is 0.  Otherwise the default for
 
283
*>                      ISPEC=16 is 2.
 
284
*> \endverbatim
 
285
*
 
286
*> \par References:
 
287
*  ================
 
288
*>
 
289
*>       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
 
290
*>       Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
 
291
*>       Performance, SIAM Journal of Matrix Analysis, volume 23, pages
 
292
*>       929--947, 2002.
 
293
*> \n
 
294
*>       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
 
295
*>       Algorithm Part II: Aggressive Early Deflation, SIAM Journal
 
296
*>       of Matrix Analysis, volume 23, pages 948--973, 2002.
 
297
*
 
298
*  =====================================================================
1
299
      SUBROUTINE ZHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ,
2
300
     $                   WORK, LWORK, INFO )
3
 
C$Id: zhseqr.f 19697 2010-10-29 16:57:34Z d3y133 $                          
4
301
*
5
 
*  -- LAPACK routine (version 3.0) --
6
 
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
7
 
*     Courant Institute, Argonne National Lab, and Rice University
8
 
*     June 30, 1999
 
302
*  -- LAPACK computational routine (version 3.4.0) --
 
303
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
304
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
305
*     November 2011
9
306
*
10
307
*     .. Scalar Arguments ..
 
308
      INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
11
309
      CHARACTER          COMPZ, JOB
12
 
      INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
13
310
*     ..
14
311
*     .. Array Arguments ..
15
312
      COMPLEX*16         H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
16
313
*     ..
17
314
*
18
 
*  Purpose
19
 
*  =======
20
 
*
21
 
*  ZHSEQR computes the eigenvalues of a complex upper Hessenberg
22
 
*  matrix H, and, optionally, the matrices T and Z from the Schur
23
 
*  decomposition H = Z T Z**H, where T is an upper triangular matrix
24
 
*  (the Schur form), and Z is the unitary matrix of Schur vectors.
25
 
*
26
 
*  Optionally Z may be postmultiplied into an input unitary matrix Q,
27
 
*  so that this routine can give the Schur factorization of a matrix A
28
 
*  which has been reduced to the Hessenberg form H by the unitary
29
 
*  matrix Q:  A = Q*H*Q**H = (QZ)*T*(QZ)**H.
30
 
*
31
 
*  Arguments
32
 
*  =========
33
 
*
34
 
*  JOB     (input) CHARACTER*1
35
 
*          = 'E': compute eigenvalues only;
36
 
*          = 'S': compute eigenvalues and the Schur form T.
37
 
*
38
 
*  COMPZ   (input) CHARACTER*1
39
 
*          = 'N': no Schur vectors are computed;
40
 
*          = 'I': Z is initialized to the unit matrix and the matrix Z
41
 
*                 of Schur vectors of H is returned;
42
 
*          = 'V': Z must contain an unitary matrix Q on entry, and
43
 
*                 the product Q*Z is returned.
44
 
*
45
 
*  N       (input) INTEGER
46
 
*          The order of the matrix H.  N >= 0.
47
 
*
48
 
*  ILO     (input) INTEGER
49
 
*  IHI     (input) INTEGER
50
 
*          It is assumed that H is already upper triangular in rows
51
 
*          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
52
 
*          set by a previous call to ZGEBAL, and then passed to CGEHRD
53
 
*          when the matrix output by ZGEBAL is reduced to Hessenberg
54
 
*          form. Otherwise ILO and IHI should be set to 1 and N
55
 
*          respectively.
56
 
*          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
57
 
*
58
 
*  H       (input/output) COMPLEX*16 array, dimension (LDH,N)
59
 
*          On entry, the upper Hessenberg matrix H.
60
 
*          On exit, if JOB = 'S', H contains the upper triangular matrix
61
 
*          T from the Schur decomposition (the Schur form). If
62
 
*          JOB = 'E', the contents of H are unspecified on exit.
63
 
*
64
 
*  LDH     (input) INTEGER
65
 
*          The leading dimension of the array H. LDH >= max(1,N).
66
 
*
67
 
*  W       (output) COMPLEX*16 array, dimension (N)
68
 
*          The computed eigenvalues. If JOB = 'S', the eigenvalues are
69
 
*          stored in the same order as on the diagonal of the Schur form
70
 
*          returned in H, with W(i) = H(i,i).
71
 
*
72
 
*  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)
73
 
*          If COMPZ = 'N': Z is not referenced.
74
 
*          If COMPZ = 'I': on entry, Z need not be set, and on exit, Z
75
 
*          contains the unitary matrix Z of the Schur vectors of H.
76
 
*          If COMPZ = 'V': on entry Z must contain an N-by-N matrix Q,
77
 
*          which is assumed to be equal to the unit matrix except for
78
 
*          the submatrix Z(ILO:IHI,ILO:IHI); on exit Z contains Q*Z.
79
 
*          Normally Q is the unitary matrix generated by ZUNGHR after
80
 
*          the call to ZGEHRD which formed the Hessenberg matrix H.
81
 
*
82
 
*  LDZ     (input) INTEGER
83
 
*          The leading dimension of the array Z.
84
 
*          LDZ >= max(1,N) if COMPZ = 'I' or 'V'; LDZ >= 1 otherwise.
85
 
*
86
 
*  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
87
 
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
88
 
*
89
 
*  LWORK   (input) INTEGER
90
 
*          The dimension of the array WORK.  LWORK >= max(1,N).
91
 
*
92
 
*          If LWORK = -1, then a workspace query is assumed; the routine
93
 
*          only calculates the optimal size of the WORK array, returns
94
 
*          this value as the first entry of the WORK array, and no error
95
 
*          message related to LWORK is issued by XERBLA.
96
 
*
97
 
*  INFO    (output) INTEGER
98
 
*          = 0:  successful exit
99
 
*          < 0:  if INFO = -i, the i-th argument had an illegal value
100
 
*          > 0:  if INFO = i, ZHSEQR failed to compute all the
101
 
*                eigenvalues in a total of 30*(IHI-ILO+1) iterations;
102
 
*                elements 1:ilo-1 and i+1:n of W contain those
103
 
*                eigenvalues which have been successfully computed.
104
 
*
105
315
*  =====================================================================
106
316
*
107
317
*     .. Parameters ..
 
318
*
 
319
*     ==== Matrices of order NTINY or smaller must be processed by
 
320
*     .    ZLAHQR because of insufficient subdiagonal scratch space.
 
321
*     .    (This is a hard limit.) ====
 
322
      INTEGER            NTINY
 
323
      PARAMETER          ( NTINY = 11 )
 
324
*
 
325
*     ==== NL allocates some local workspace to help small matrices
 
326
*     .    through a rare ZLAHQR failure.  NL .GT. NTINY = 11 is
 
327
*     .    required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
 
328
*     .    mended.  (The default value of NMIN is 75.)  Using NL = 49
 
329
*     .    allows up to six simultaneous shifts and a 16-by-16
 
330
*     .    deflation window.  ====
 
331
      INTEGER            NL
 
332
      PARAMETER          ( NL = 49 )
108
333
      COMPLEX*16         ZERO, ONE
109
 
      PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
110
 
     $                   ONE = ( 1.0D+0, 0.0D+0 ) )
111
 
      DOUBLE PRECISION   RZERO, RONE, CONST
112
 
      PARAMETER          ( RZERO = 0.0D+0, RONE = 1.0D+0,
113
 
     $                   CONST = 1.5D+0 )
114
 
      INTEGER            NSMAX, LDS
115
 
      PARAMETER          ( NSMAX = 15, LDS = NSMAX )
 
334
      PARAMETER          ( ZERO = ( 0.0d0, 0.0d0 ),
 
335
     $                   ONE = ( 1.0d0, 0.0d0 ) )
 
336
      DOUBLE PRECISION   RZERO
 
337
      PARAMETER          ( RZERO = 0.0d0 )
 
338
*     ..
 
339
*     .. Local Arrays ..
 
340
      COMPLEX*16         HL( NL, NL ), WORKL( NL )
116
341
*     ..
117
342
*     .. Local Scalars ..
 
343
      INTEGER            KBOT, NMIN
118
344
      LOGICAL            INITZ, LQUERY, WANTT, WANTZ
119
 
      INTEGER            I, I1, I2, IERR, II, ITEMP, ITN, ITS, J, K, L,
120
 
     $                   MAXB, NH, NR, NS, NV
121
 
      DOUBLE PRECISION   OVFL, RTEMP, SMLNUM, TST1, ULP, UNFL
122
 
      COMPLEX*16         CDUM, TAU, TEMP
123
 
*     ..
124
 
*     .. Local Arrays ..
125
 
      DOUBLE PRECISION   RWORK( 1 )
126
 
      COMPLEX*16         S( LDS, NSMAX ), V( NSMAX+1 ), VV( NSMAX+1 )
127
345
*     ..
128
346
*     .. External Functions ..
 
347
      INTEGER            ILAENV
129
348
      LOGICAL            LSAME
130
 
      INTEGER            ILAENV, IZAMAX
131
 
      DOUBLE PRECISION   DLAMCH, DLAPY2, ZLANHS
132
 
      EXTERNAL           LSAME, ILAENV, IZAMAX, DLAMCH, DLAPY2, ZLANHS
 
349
      EXTERNAL           ILAENV, LSAME
133
350
*     ..
134
351
*     .. External Subroutines ..
135
 
      EXTERNAL           XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLACPY, ZLAHQR,
136
 
     $                   ZLARFG, ZLARFX, ZLASET, ZSCAL
 
352
      EXTERNAL           XERBLA, ZCOPY, ZLACPY, ZLAHQR, ZLAQR0, ZLASET
137
353
*     ..
138
354
*     .. Intrinsic Functions ..
139
 
      INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX, MIN
140
 
*     ..
141
 
*     .. Statement Functions ..
142
 
      DOUBLE PRECISION   CABS1
143
 
*     ..
144
 
*     .. Statement Function definitions ..
145
 
      CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
 
355
      INTRINSIC          DBLE, DCMPLX, MAX, MIN
146
356
*     ..
147
357
*     .. Executable Statements ..
148
358
*
149
 
*     Decode and test the input parameters
 
359
*     ==== Decode and check the input parameters. ====
150
360
*
151
361
      WANTT = LSAME( JOB, 'S' )
152
362
      INITZ = LSAME( COMPZ, 'I' )
153
363
      WANTZ = INITZ .OR. LSAME( COMPZ, 'V' )
 
364
      WORK( 1 ) = DCMPLX( DBLE( MAX( 1, N ) ), RZERO )
 
365
      LQUERY = LWORK.EQ.-1
154
366
*
155
367
      INFO = 0
156
 
      WORK( 1 ) = MAX( 1, N )
157
 
      LQUERY = ( LWORK.EQ.-1 )
158
368
      IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN
159
369
         INFO = -1
160
370
      ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
167
377
         INFO = -5
168
378
      ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
169
379
         INFO = -7
170
 
      ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. LDZ.LT.MAX( 1, N ) ) THEN
 
380
      ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
171
381
         INFO = -10
172
382
      ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
173
383
         INFO = -12
174
384
      END IF
 
385
*
175
386
      IF( INFO.NE.0 ) THEN
 
387
*
 
388
*        ==== Quick return in case of invalid argument. ====
 
389
*
176
390
         CALL XERBLA( 'ZHSEQR', -INFO )
177
391
         RETURN
 
392
*
 
393
      ELSE IF( N.EQ.0 ) THEN
 
394
*
 
395
*        ==== Quick return in case N = 0; nothing to do. ====
 
396
*
 
397
         RETURN
 
398
*
178
399
      ELSE IF( LQUERY ) THEN
179
 
         RETURN
180
 
      END IF
181
 
*
182
 
*     Initialize Z, if necessary
183
 
*
184
 
      IF( INITZ )
185
 
     $   CALL ZLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
186
 
*
187
 
*     Store the eigenvalues isolated by ZGEBAL.
188
 
*
189
 
      DO 10 I = 1, ILO - 1
190
 
         W( I ) = H( I, I )
191
 
   10 CONTINUE
192
 
      DO 20 I = IHI + 1, N
193
 
         W( I ) = H( I, I )
194
 
   20 CONTINUE
195
 
*
196
 
*     Quick return if possible.
197
 
*
198
 
      IF( N.EQ.0 )
199
 
     $   RETURN
200
 
      IF( ILO.EQ.IHI ) THEN
201
 
         W( ILO ) = H( ILO, ILO )
202
 
         RETURN
203
 
      END IF
204
 
*
205
 
*     Set rows and columns ILO to IHI to zero below the first
206
 
*     subdiagonal.
207
 
*
208
 
      DO 40 J = ILO, IHI - 2
209
 
         DO 30 I = J + 2, N
210
 
            H( I, J ) = ZERO
211
 
   30    CONTINUE
212
 
   40 CONTINUE
213
 
      NH = IHI - ILO + 1
214
 
*
215
 
*     I1 and I2 are the indices of the first row and last column of H
216
 
*     to which transformations must be applied. If eigenvalues only are
217
 
*     being computed, I1 and I2 are re-set inside the main loop.
218
 
*
219
 
      IF( WANTT ) THEN
220
 
         I1 = 1
221
 
         I2 = N
 
400
*
 
401
*        ==== Quick return in case of a workspace query ====
 
402
*
 
403
         CALL ZLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILO, IHI, Z,
 
404
     $                LDZ, WORK, LWORK, INFO )
 
405
*        ==== Ensure reported workspace size is backward-compatible with
 
406
*        .    previous LAPACK versions. ====
 
407
         WORK( 1 ) = DCMPLX( MAX( DBLE( WORK( 1 ) ), DBLE( MAX( 1,
 
408
     $               N ) ) ), RZERO )
 
409
         RETURN
 
410
*
222
411
      ELSE
223
 
         I1 = ILO
224
 
         I2 = IHI
225
 
      END IF
226
 
*
227
 
*     Ensure that the subdiagonal elements are real.
228
 
*
229
 
      DO 50 I = ILO + 1, IHI
230
 
         TEMP = H( I, I-1 )
231
 
         IF( DIMAG( TEMP ).NE.RZERO ) THEN
232
 
            RTEMP = DLAPY2( DBLE( TEMP ), DIMAG( TEMP ) )
233
 
            H( I, I-1 ) = RTEMP
234
 
            TEMP = TEMP / RTEMP
235
 
            IF( I2.GT.I )
236
 
     $         CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH )
237
 
            CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 )
238
 
            IF( I.LT.IHI )
239
 
     $         H( I+1, I ) = TEMP*H( I+1, I )
240
 
            IF( WANTZ )
241
 
     $         CALL ZSCAL( NH, TEMP, Z( ILO, I ), 1 )
242
 
         END IF
243
 
   50 CONTINUE
244
 
*
245
 
*     Determine the order of the multi-shift QR algorithm to be used.
246
 
*
247
 
      NS = ILAENV( 4, 'ZHSEQR', JOB // COMPZ, N, ILO, IHI, -1 )
248
 
      MAXB = ILAENV( 8, 'ZHSEQR', JOB // COMPZ, N, ILO, IHI, -1 )
249
 
      IF( NS.LE.1 .OR. NS.GT.NH .OR. MAXB.GE.NH ) THEN
250
 
*
251
 
*        Use the standard double-shift algorithm
252
 
*
253
 
         CALL ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILO, IHI, Z,
254
 
     $                LDZ, INFO )
255
 
         RETURN
256
 
      END IF
257
 
      MAXB = MAX( 2, MAXB )
258
 
      NS = MIN( NS, MAXB, NSMAX )
259
 
*
260
 
*     Now 1 < NS <= MAXB < NH.
261
 
*
262
 
*     Set machine-dependent constants for the stopping criterion.
263
 
*     If norm(H) <= sqrt(OVFL), overflow should not occur.
264
 
*
265
 
      UNFL = DLAMCH( 'Safe minimum' )
266
 
      OVFL = RONE / UNFL
267
 
      CALL DLABAD( UNFL, OVFL )
268
 
      ULP = DLAMCH( 'Precision' )
269
 
      SMLNUM = UNFL*( NH / ULP )
270
 
*
271
 
*     ITN is the total number of multiple-shift QR iterations allowed.
272
 
*
273
 
      ITN = 30*NH
274
 
*
275
 
*     The main loop begins here. I is the loop index and decreases from
276
 
*     IHI to ILO in steps of at most MAXB. Each iteration of the loop
277
 
*     works with the active submatrix in rows and columns L to I.
278
 
*     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
279
 
*     H(L,L-1) is negligible so that the matrix splits.
280
 
*
281
 
      I = IHI
282
 
   60 CONTINUE
283
 
      IF( I.LT.ILO )
284
 
     $   GO TO 180
285
 
*
286
 
*     Perform multiple-shift QR iterations on rows and columns ILO to I
287
 
*     until a submatrix of order at most MAXB splits off at the bottom
288
 
*     because a subdiagonal element has become negligible.
289
 
*
290
 
      L = ILO
291
 
      DO 160 ITS = 0, ITN
292
 
*
293
 
*        Look for a single small subdiagonal element.
294
 
*
295
 
         DO 70 K = I, L + 1, -1
296
 
            TST1 = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
297
 
            IF( TST1.EQ.RZERO )
298
 
     $         TST1 = ZLANHS( '1', I-L+1, H( L, L ), LDH, RWORK )
299
 
            IF( ABS( DBLE( H( K, K-1 ) ) ).LE.MAX( ULP*TST1, SMLNUM ) )
300
 
     $         GO TO 80
301
 
   70    CONTINUE
302
 
   80    CONTINUE
303
 
         L = K
304
 
         IF( L.GT.ILO ) THEN
305
 
*
306
 
*           H(L,L-1) is negligible.
307
 
*
308
 
            H( L, L-1 ) = ZERO
309
 
         END IF
310
 
*
311
 
*        Exit from loop if a submatrix of order <= MAXB has split off.
312
 
*
313
 
         IF( L.GE.I-MAXB+1 )
314
 
     $      GO TO 170
315
 
*
316
 
*        Now the active submatrix is in rows and columns L to I. If
317
 
*        eigenvalues only are being computed, only the active submatrix
318
 
*        need be transformed.
319
 
*
320
 
         IF( .NOT.WANTT ) THEN
321
 
            I1 = L
322
 
            I2 = I
323
 
         END IF
324
 
*
325
 
         IF( ITS.EQ.20 .OR. ITS.EQ.30 ) THEN
326
 
*
327
 
*           Exceptional shifts.
328
 
*
329
 
            DO 90 II = I - NS + 1, I
330
 
               W( II ) = CONST*( ABS( DBLE( H( II, II-1 ) ) )+
331
 
     $                   ABS( DBLE( H( II, II ) ) ) )
332
 
   90       CONTINUE
 
412
*
 
413
*        ==== copy eigenvalues isolated by ZGEBAL ====
 
414
*
 
415
         IF( ILO.GT.1 )
 
416
     $      CALL ZCOPY( ILO-1, H, LDH+1, W, 1 )
 
417
         IF( IHI.LT.N )
 
418
     $      CALL ZCOPY( N-IHI, H( IHI+1, IHI+1 ), LDH+1, W( IHI+1 ), 1 )
 
419
*
 
420
*        ==== Initialize Z, if requested ====
 
421
*
 
422
         IF( INITZ )
 
423
     $      CALL ZLASET( 'A', N, N, ZERO, ONE, Z, LDZ )
 
424
*
 
425
*        ==== Quick return if possible ====
 
426
*
 
427
         IF( ILO.EQ.IHI ) THEN
 
428
            W( ILO ) = H( ILO, ILO )
 
429
            RETURN
 
430
         END IF
 
431
*
 
432
*        ==== ZLAHQR/ZLAQR0 crossover point ====
 
433
*
 
434
         NMIN = ILAENV( 12, 'ZHSEQR', JOB( : 1 ) // COMPZ( : 1 ), N,
 
435
     $          ILO, IHI, LWORK )
 
436
         NMIN = MAX( NTINY, NMIN )
 
437
*
 
438
*        ==== ZLAQR0 for big matrices; ZLAHQR for small ones ====
 
439
*
 
440
         IF( N.GT.NMIN ) THEN
 
441
            CALL ZLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILO, IHI,
 
442
     $                   Z, LDZ, WORK, LWORK, INFO )
333
443
         ELSE
334
444
*
335
 
*           Use eigenvalues of trailing submatrix of order NS as shifts.
336
 
*
337
 
            CALL ZLACPY( 'Full', NS, NS, H( I-NS+1, I-NS+1 ), LDH, S,
338
 
     $                   LDS )
339
 
            CALL ZLAHQR( .FALSE., .FALSE., NS, 1, NS, S, LDS,
340
 
     $                   W( I-NS+1 ), 1, NS, Z, LDZ, IERR )
341
 
            IF( IERR.GT.0 ) THEN
342
 
*
343
 
*              If ZLAHQR failed to compute all NS eigenvalues, use the
344
 
*              unconverged diagonal elements as the remaining shifts.
345
 
*
346
 
               DO 100 II = 1, IERR
347
 
                  W( I-NS+II ) = S( II, II )
348
 
  100          CONTINUE
349
 
            END IF
350
 
         END IF
351
 
*
352
 
*        Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns))
353
 
*        where G is the Hessenberg submatrix H(L:I,L:I) and w is
354
 
*        the vector of shifts (stored in W). The result is
355
 
*        stored in the local array V.
356
 
*
357
 
         V( 1 ) = ONE
358
 
         DO 110 II = 2, NS + 1
359
 
            V( II ) = ZERO
360
 
  110    CONTINUE
361
 
         NV = 1
362
 
         DO 130 J = I - NS + 1, I
363
 
            CALL ZCOPY( NV+1, V, 1, VV, 1 )
364
 
            CALL ZGEMV( 'No transpose', NV+1, NV, ONE, H( L, L ), LDH,
365
 
     $                  VV, 1, -W( J ), V, 1 )
366
 
            NV = NV + 1
367
 
*
368
 
*           Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero,
369
 
*           reset it to the unit vector.
370
 
*
371
 
            ITEMP = IZAMAX( NV, V, 1 )
372
 
            RTEMP = CABS1( V( ITEMP ) )
373
 
            IF( RTEMP.EQ.RZERO ) THEN
374
 
               V( 1 ) = ONE
375
 
               DO 120 II = 2, NV
376
 
                  V( II ) = ZERO
377
 
  120          CONTINUE
378
 
            ELSE
379
 
               RTEMP = MAX( RTEMP, SMLNUM )
380
 
               CALL ZDSCAL( NV, RONE / RTEMP, V, 1 )
381
 
            END IF
382
 
  130    CONTINUE
383
 
*
384
 
*        Multiple-shift QR step
385
 
*
386
 
         DO 150 K = L, I - 1
387
 
*
388
 
*           The first iteration of this loop determines a reflection G
389
 
*           from the vector V and applies it from left and right to H,
390
 
*           thus creating a nonzero bulge below the subdiagonal.
391
 
*
392
 
*           Each subsequent iteration determines a reflection G to
393
 
*           restore the Hessenberg form in the (K-1)th column, and thus
394
 
*           chases the bulge one step toward the bottom of the active
395
 
*           submatrix. NR is the order of G.
396
 
*
397
 
            NR = MIN( NS+1, I-K+1 )
398
 
            IF( K.GT.L )
399
 
     $         CALL ZCOPY( NR, H( K, K-1 ), 1, V, 1 )
400
 
            CALL ZLARFG( NR, V( 1 ), V( 2 ), 1, TAU )
401
 
            IF( K.GT.L ) THEN
402
 
               H( K, K-1 ) = V( 1 )
403
 
               DO 140 II = K + 1, I
404
 
                  H( II, K-1 ) = ZERO
405
 
  140          CONTINUE
406
 
            END IF
407
 
            V( 1 ) = ONE
408
 
*
409
 
*           Apply G' from the left to transform the rows of the matrix
410
 
*           in columns K to I2.
411
 
*
412
 
            CALL ZLARFX( 'Left', NR, I2-K+1, V, DCONJG( TAU ),
413
 
     $                   H( K, K ), LDH, WORK )
414
 
*
415
 
*           Apply G from the right to transform the columns of the
416
 
*           matrix in rows I1 to min(K+NR,I).
417
 
*
418
 
            CALL ZLARFX( 'Right', MIN( K+NR, I )-I1+1, NR, V, TAU,
419
 
     $                   H( I1, K ), LDH, WORK )
420
 
*
421
 
            IF( WANTZ ) THEN
422
 
*
423
 
*              Accumulate transformations in the matrix Z
424
 
*
425
 
               CALL ZLARFX( 'Right', NH, NR, V, TAU, Z( ILO, K ), LDZ,
426
 
     $                      WORK )
427
 
            END IF
428
 
  150    CONTINUE
429
 
*
430
 
*        Ensure that H(I,I-1) is real.
431
 
*
432
 
         TEMP = H( I, I-1 )
433
 
         IF( DIMAG( TEMP ).NE.RZERO ) THEN
434
 
            RTEMP = DLAPY2( DBLE( TEMP ), DIMAG( TEMP ) )
435
 
            H( I, I-1 ) = RTEMP
436
 
            TEMP = TEMP / RTEMP
437
 
            IF( I2.GT.I )
438
 
     $         CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH )
439
 
            CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 )
440
 
            IF( WANTZ ) THEN
441
 
               CALL ZSCAL( NH, TEMP, Z( ILO, I ), 1 )
442
 
            END IF
443
 
         END IF
444
 
*
445
 
  160 CONTINUE
446
 
*
447
 
*     Failure to converge in remaining number of iterations
448
 
*
449
 
      INFO = I
450
 
      RETURN
451
 
*
452
 
  170 CONTINUE
453
 
*
454
 
*     A submatrix of order <= MAXB in rows and columns L to I has split
455
 
*     off. Use the double-shift QR algorithm to handle it.
456
 
*
457
 
      CALL ZLAHQR( WANTT, WANTZ, N, L, I, H, LDH, W, ILO, IHI, Z, LDZ,
458
 
     $             INFO )
459
 
      IF( INFO.GT.0 )
460
 
     $   RETURN
461
 
*
462
 
*     Decrement number of remaining iterations, and return to start of
463
 
*     the main loop with a new value of I.
464
 
*
465
 
      ITN = ITN - ITS
466
 
      I = L - 1
467
 
      GO TO 60
468
 
*
469
 
  180 CONTINUE
470
 
      WORK( 1 ) = MAX( 1, N )
471
 
      RETURN
472
 
*
473
 
*     End of ZHSEQR
 
445
*           ==== Small matrix ====
 
446
*
 
447
            CALL ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILO, IHI,
 
448
     $                   Z, LDZ, INFO )
 
449
*
 
450
            IF( INFO.GT.0 ) THEN
 
451
*
 
452
*              ==== A rare ZLAHQR failure!  ZLAQR0 sometimes succeeds
 
453
*              .    when ZLAHQR fails. ====
 
454
*
 
455
               KBOT = INFO
 
456
*
 
457
               IF( N.GE.NL ) THEN
 
458
*
 
459
*                 ==== Larger matrices have enough subdiagonal scratch
 
460
*                 .    space to call ZLAQR0 directly. ====
 
461
*
 
462
                  CALL ZLAQR0( WANTT, WANTZ, N, ILO, KBOT, H, LDH, W,
 
463
     $                         ILO, IHI, Z, LDZ, WORK, LWORK, INFO )
 
464
*
 
465
               ELSE
 
466
*
 
467
*                 ==== Tiny matrices don't have enough subdiagonal
 
468
*                 .    scratch space to benefit from ZLAQR0.  Hence,
 
469
*                 .    tiny matrices must be copied into a larger
 
470
*                 .    array before calling ZLAQR0. ====
 
471
*
 
472
                  CALL ZLACPY( 'A', N, N, H, LDH, HL, NL )
 
473
                  HL( N+1, N ) = ZERO
 
474
                  CALL ZLASET( 'A', NL, NL-N, ZERO, ZERO, HL( 1, N+1 ),
 
475
     $                         NL )
 
476
                  CALL ZLAQR0( WANTT, WANTZ, NL, ILO, KBOT, HL, NL, W,
 
477
     $                         ILO, IHI, Z, LDZ, WORKL, NL, INFO )
 
478
                  IF( WANTT .OR. INFO.NE.0 )
 
479
     $               CALL ZLACPY( 'A', N, N, HL, NL, H, LDH )
 
480
               END IF
 
481
            END IF
 
482
         END IF
 
483
*
 
484
*        ==== Clear out the trash, if necessary. ====
 
485
*
 
486
         IF( ( WANTT .OR. INFO.NE.0 ) .AND. N.GT.2 )
 
487
     $      CALL ZLASET( 'L', N-2, N-2, ZERO, ZERO, H( 3, 1 ), LDH )
 
488
*
 
489
*        ==== Ensure reported workspace size is backward-compatible with
 
490
*        .    previous LAPACK versions. ====
 
491
*
 
492
         WORK( 1 ) = DCMPLX( MAX( DBLE( MAX( 1, N ) ),
 
493
     $               DBLE( WORK( 1 ) ) ), RZERO )
 
494
      END IF
 
495
*
 
496
*     ==== End of ZHSEQR ====
474
497
*
475
498
      END