LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zlahqr.f
Go to the documentation of this file.
1 *> \brief \b ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZLAHQR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahqr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahqr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahqr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
22 * IHIZ, Z, LDZ, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
26 * LOGICAL WANTT, WANTZ
27 * ..
28 * .. Array Arguments ..
29 * COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZLAHQR is an auxiliary routine called by CHSEQR to update the
39 *> eigenvalues and Schur decomposition already computed by CHSEQR, by
40 *> dealing with the Hessenberg submatrix in rows and columns ILO to
41 *> IHI.
42 *> \endverbatim
43 *
44 * Arguments:
45 * ==========
46 *
47 *> \param[in] WANTT
48 *> \verbatim
49 *> WANTT is LOGICAL
50 *> = .TRUE. : the full Schur form T is required;
51 *> = .FALSE.: only eigenvalues are required.
52 *> \endverbatim
53 *>
54 *> \param[in] WANTZ
55 *> \verbatim
56 *> WANTZ is LOGICAL
57 *> = .TRUE. : the matrix of Schur vectors Z is required;
58 *> = .FALSE.: Schur vectors are not required.
59 *> \endverbatim
60 *>
61 *> \param[in] N
62 *> \verbatim
63 *> N is INTEGER
64 *> The order of the matrix H. N >= 0.
65 *> \endverbatim
66 *>
67 *> \param[in] ILO
68 *> \verbatim
69 *> ILO is INTEGER
70 *> \endverbatim
71 *>
72 *> \param[in] IHI
73 *> \verbatim
74 *> IHI is INTEGER
75 *> It is assumed that H is already upper triangular in rows and
76 *> columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
77 *> ZLAHQR works primarily with the Hessenberg submatrix in rows
78 *> and columns ILO to IHI, but applies transformations to all of
79 *> H if WANTT is .TRUE..
80 *> 1 <= ILO <= max(1,IHI); IHI <= N.
81 *> \endverbatim
82 *>
83 *> \param[in,out] H
84 *> \verbatim
85 *> H is COMPLEX*16 array, dimension (LDH,N)
86 *> On entry, the upper Hessenberg matrix H.
87 *> On exit, if INFO is zero and if WANTT is .TRUE., then H
88 *> is upper triangular in rows and columns ILO:IHI. If INFO
89 *> is zero and if WANTT is .FALSE., then the contents of H
90 *> are unspecified on exit. The output state of H in case
91 *> INF is positive is below under the description of INFO.
92 *> \endverbatim
93 *>
94 *> \param[in] LDH
95 *> \verbatim
96 *> LDH is INTEGER
97 *> The leading dimension of the array H. LDH >= max(1,N).
98 *> \endverbatim
99 *>
100 *> \param[out] W
101 *> \verbatim
102 *> W is COMPLEX*16 array, dimension (N)
103 *> The computed eigenvalues ILO to IHI are stored in the
104 *> corresponding elements of W. If WANTT is .TRUE., the
105 *> eigenvalues are stored in the same order as on the diagonal
106 *> of the Schur form returned in H, with W(i) = H(i,i).
107 *> \endverbatim
108 *>
109 *> \param[in] ILOZ
110 *> \verbatim
111 *> ILOZ is INTEGER
112 *> \endverbatim
113 *>
114 *> \param[in] IHIZ
115 *> \verbatim
116 *> IHIZ is INTEGER
117 *> Specify the rows of Z to which transformations must be
118 *> applied if WANTZ is .TRUE..
119 *> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
120 *> \endverbatim
121 *>
122 *> \param[in,out] Z
123 *> \verbatim
124 *> Z is COMPLEX*16 array, dimension (LDZ,N)
125 *> If WANTZ is .TRUE., on entry Z must contain the current
126 *> matrix Z of transformations accumulated by CHSEQR, and on
127 *> exit Z has been updated; transformations are applied only to
128 *> the submatrix Z(ILOZ:IHIZ,ILO:IHI).
129 *> If WANTZ is .FALSE., Z is not referenced.
130 *> \endverbatim
131 *>
132 *> \param[in] LDZ
133 *> \verbatim
134 *> LDZ is INTEGER
135 *> The leading dimension of the array Z. LDZ >= max(1,N).
136 *> \endverbatim
137 *>
138 *> \param[out] INFO
139 *> \verbatim
140 *> INFO is INTEGER
141 *> = 0: successful exit
142 *> .GT. 0: if INFO = i, ZLAHQR failed to compute all the
143 *> eigenvalues ILO to IHI in a total of 30 iterations
144 *> per eigenvalue; elements i+1:ihi of W contain
145 *> those eigenvalues which have been successfully
146 *> computed.
147 *>
148 *> If INFO .GT. 0 and WANTT is .FALSE., then on exit,
149 *> the remaining unconverged eigenvalues are the
150 *> eigenvalues of the upper Hessenberg matrix
151 *> rows and columns ILO thorugh INFO of the final,
152 *> output value of H.
153 *>
154 *> If INFO .GT. 0 and WANTT is .TRUE., then on exit
155 *> (*) (initial value of H)*U = U*(final value of H)
156 *> where U is an orthognal matrix. The final
157 *> value of H is upper Hessenberg and triangular in
158 *> rows and columns INFO+1 through IHI.
159 *>
160 *> If INFO .GT. 0 and WANTZ is .TRUE., then on exit
161 *> (final value of Z) = (initial value of Z)*U
162 *> where U is the orthogonal matrix in (*)
163 *> (regardless of the value of WANTT.)
164 *> \endverbatim
165 *
166 * Authors:
167 * ========
168 *
169 *> \author Univ. of Tennessee
170 *> \author Univ. of California Berkeley
171 *> \author Univ. of Colorado Denver
172 *> \author NAG Ltd.
173 *
174 *> \date September 2012
175 *
176 *> \ingroup complex16OTHERauxiliary
177 *
178 *> \par Contributors:
179 * ==================
180 *>
181 *> \verbatim
182 *>
183 *> 02-96 Based on modifications by
184 *> David Day, Sandia National Laboratory, USA
185 *>
186 *> 12-04 Further modifications by
187 *> Ralph Byers, University of Kansas, USA
188 *> This is a modified version of ZLAHQR from LAPACK version 3.0.
189 *> It is (1) more robust against overflow and underflow and
190 *> (2) adopts the more conservative Ahues & Tisseur stopping
191 *> criterion (LAWN 122, 1997).
192 *> \endverbatim
193 *>
194 * =====================================================================
195  SUBROUTINE zlahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
196  $ ihiz, z, ldz, info )
197 *
198 * -- LAPACK auxiliary routine (version 3.4.2) --
199 * -- LAPACK is a software package provided by Univ. of Tennessee, --
200 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201 * September 2012
202 *
203 * .. Scalar Arguments ..
204  INTEGER ihi, ihiz, ilo, iloz, info, ldh, ldz, n
205  LOGICAL wantt, wantz
206 * ..
207 * .. Array Arguments ..
208  COMPLEX*16 h( ldh, * ), w( * ), z( ldz, * )
209 * ..
210 *
211 * =========================================================
212 *
213 * .. Parameters ..
214  INTEGER itmax
215  parameter( itmax = 30 )
216  COMPLEX*16 zero, one
217  parameter( zero = ( 0.0d0, 0.0d0 ),
218  $ one = ( 1.0d0, 0.0d0 ) )
219  DOUBLE PRECISION rzero, rone, half
220  parameter( rzero = 0.0d0, rone = 1.0d0, half = 0.5d0 )
221  DOUBLE PRECISION dat1
222  parameter( dat1 = 3.0d0 / 4.0d0 )
223 * ..
224 * .. Local Scalars ..
225  COMPLEX*16 cdum, h11, h11s, h22, sc, sum, t, t1, temp, u,
226  $ v2, x, y
227  DOUBLE PRECISION aa, ab, ba, bb, h10, h21, rtemp, s, safmax,
228  $ safmin, smlnum, sx, t2, tst, ulp
229  INTEGER i, i1, i2, its, j, jhi, jlo, k, l, m, nh, nz
230 * ..
231 * .. Local Arrays ..
232  COMPLEX*16 v( 2 )
233 * ..
234 * .. External Functions ..
235  COMPLEX*16 zladiv
236  DOUBLE PRECISION dlamch
237  EXTERNAL zladiv, dlamch
238 * ..
239 * .. External Subroutines ..
240  EXTERNAL dlabad, zcopy, zlarfg, zscal
241 * ..
242 * .. Statement Functions ..
243  DOUBLE PRECISION cabs1
244 * ..
245 * .. Intrinsic Functions ..
246  INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
247 * ..
248 * .. Statement Function definitions ..
249  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
250 * ..
251 * .. Executable Statements ..
252 *
253  info = 0
254 *
255 * Quick return if possible
256 *
257  IF( n.EQ.0 )
258  $ return
259  IF( ilo.EQ.ihi ) THEN
260  w( ilo ) = h( ilo, ilo )
261  return
262  END IF
263 *
264 * ==== clear out the trash ====
265  DO 10 j = ilo, ihi - 3
266  h( j+2, j ) = zero
267  h( j+3, j ) = zero
268  10 continue
269  IF( ilo.LE.ihi-2 )
270  $ h( ihi, ihi-2 ) = zero
271 * ==== ensure that subdiagonal entries are real ====
272  IF( wantt ) THEN
273  jlo = 1
274  jhi = n
275  ELSE
276  jlo = ilo
277  jhi = ihi
278  END IF
279  DO 20 i = ilo + 1, ihi
280  IF( dimag( h( i, i-1 ) ).NE.rzero ) THEN
281 * ==== The following redundant normalization
282 * . avoids problems with both gradual and
283 * . sudden underflow in ABS(H(I,I-1)) ====
284  sc = h( i, i-1 ) / cabs1( h( i, i-1 ) )
285  sc = dconjg( sc ) / abs( sc )
286  h( i, i-1 ) = abs( h( i, i-1 ) )
287  CALL zscal( jhi-i+1, sc, h( i, i ), ldh )
288  CALL zscal( min( jhi, i+1 )-jlo+1, dconjg( sc ),
289  $ h( jlo, i ), 1 )
290  IF( wantz )
291  $ CALL zscal( ihiz-iloz+1, dconjg( sc ), z( iloz, i ), 1 )
292  END IF
293  20 continue
294 *
295  nh = ihi - ilo + 1
296  nz = ihiz - iloz + 1
297 *
298 * Set machine-dependent constants for the stopping criterion.
299 *
300  safmin = dlamch( 'SAFE MINIMUM' )
301  safmax = rone / safmin
302  CALL dlabad( safmin, safmax )
303  ulp = dlamch( 'PRECISION' )
304  smlnum = safmin*( dble( nh ) / ulp )
305 *
306 * I1 and I2 are the indices of the first row and last column of H
307 * to which transformations must be applied. If eigenvalues only are
308 * being computed, I1 and I2 are set inside the main loop.
309 *
310  IF( wantt ) THEN
311  i1 = 1
312  i2 = n
313  END IF
314 *
315 * The main loop begins here. I is the loop index and decreases from
316 * IHI to ILO in steps of 1. Each iteration of the loop works
317 * with the active submatrix in rows and columns L to I.
318 * Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
319 * H(L,L-1) is negligible so that the matrix splits.
320 *
321  i = ihi
322  30 continue
323  IF( i.LT.ilo )
324  $ go to 150
325 *
326 * Perform QR iterations on rows and columns ILO to I until a
327 * submatrix of order 1 splits off at the bottom because a
328 * subdiagonal element has become negligible.
329 *
330  l = ilo
331  DO 130 its = 0, itmax
332 *
333 * Look for a single small subdiagonal element.
334 *
335  DO 40 k = i, l + 1, -1
336  IF( cabs1( h( k, k-1 ) ).LE.smlnum )
337  $ go to 50
338  tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
339  IF( tst.EQ.zero ) THEN
340  IF( k-2.GE.ilo )
341  $ tst = tst + abs( dble( h( k-1, k-2 ) ) )
342  IF( k+1.LE.ihi )
343  $ tst = tst + abs( dble( h( k+1, k ) ) )
344  END IF
345 * ==== The following is a conservative small subdiagonal
346 * . deflation criterion due to Ahues & Tisseur (LAWN 122,
347 * . 1997). It has better mathematical foundation and
348 * . improves accuracy in some examples. ====
349  IF( abs( dble( h( k, k-1 ) ) ).LE.ulp*tst ) THEN
350  ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
351  ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
352  aa = max( cabs1( h( k, k ) ),
353  $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
354  bb = min( cabs1( h( k, k ) ),
355  $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
356  s = aa + ab
357  IF( ba*( ab / s ).LE.max( smlnum,
358  $ ulp*( bb*( aa / s ) ) ) )go to 50
359  END IF
360  40 continue
361  50 continue
362  l = k
363  IF( l.GT.ilo ) THEN
364 *
365 * H(L,L-1) is negligible
366 *
367  h( l, l-1 ) = zero
368  END IF
369 *
370 * Exit from loop if a submatrix of order 1 has split off.
371 *
372  IF( l.GE.i )
373  $ go to 140
374 *
375 * Now the active submatrix is in rows and columns L to I. If
376 * eigenvalues only are being computed, only the active submatrix
377 * need be transformed.
378 *
379  IF( .NOT.wantt ) THEN
380  i1 = l
381  i2 = i
382  END IF
383 *
384  IF( its.EQ.10 ) THEN
385 *
386 * Exceptional shift.
387 *
388  s = dat1*abs( dble( h( l+1, l ) ) )
389  t = s + h( l, l )
390  ELSE IF( its.EQ.20 ) THEN
391 *
392 * Exceptional shift.
393 *
394  s = dat1*abs( dble( h( i, i-1 ) ) )
395  t = s + h( i, i )
396  ELSE
397 *
398 * Wilkinson's shift.
399 *
400  t = h( i, i )
401  u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
402  s = cabs1( u )
403  IF( s.NE.rzero ) THEN
404  x = half*( h( i-1, i-1 )-t )
405  sx = cabs1( x )
406  s = max( s, cabs1( x ) )
407  y = s*sqrt( ( x / s )**2+( u / s )**2 )
408  IF( sx.GT.rzero ) THEN
409  IF( dble( x / sx )*dble( y )+dimag( x / sx )*
410  $ dimag( y ).LT.rzero )y = -y
411  END IF
412  t = t - u*zladiv( u, ( x+y ) )
413  END IF
414  END IF
415 *
416 * Look for two consecutive small subdiagonal elements.
417 *
418  DO 60 m = i - 1, l + 1, -1
419 *
420 * Determine the effect of starting the single-shift QR
421 * iteration at row M, and see if this would make H(M,M-1)
422 * negligible.
423 *
424  h11 = h( m, m )
425  h22 = h( m+1, m+1 )
426  h11s = h11 - t
427  h21 = dble( h( m+1, m ) )
428  s = cabs1( h11s ) + abs( h21 )
429  h11s = h11s / s
430  h21 = h21 / s
431  v( 1 ) = h11s
432  v( 2 ) = h21
433  h10 = dble( h( m, m-1 ) )
434  IF( abs( h10 )*abs( h21 ).LE.ulp*
435  $ ( cabs1( h11s )*( cabs1( h11 )+cabs1( h22 ) ) ) )
436  $ go to 70
437  60 continue
438  h11 = h( l, l )
439  h22 = h( l+1, l+1 )
440  h11s = h11 - t
441  h21 = dble( h( l+1, l ) )
442  s = cabs1( h11s ) + abs( h21 )
443  h11s = h11s / s
444  h21 = h21 / s
445  v( 1 ) = h11s
446  v( 2 ) = h21
447  70 continue
448 *
449 * Single-shift QR step
450 *
451  DO 120 k = m, i - 1
452 *
453 * The first iteration of this loop determines a reflection G
454 * from the vector V and applies it from left and right to H,
455 * thus creating a nonzero bulge below the subdiagonal.
456 *
457 * Each subsequent iteration determines a reflection G to
458 * restore the Hessenberg form in the (K-1)th column, and thus
459 * chases the bulge one step toward the bottom of the active
460 * submatrix.
461 *
462 * V(2) is always real before the call to ZLARFG, and hence
463 * after the call T2 ( = T1*V(2) ) is also real.
464 *
465  IF( k.GT.m )
466  $ CALL zcopy( 2, h( k, k-1 ), 1, v, 1 )
467  CALL zlarfg( 2, v( 1 ), v( 2 ), 1, t1 )
468  IF( k.GT.m ) THEN
469  h( k, k-1 ) = v( 1 )
470  h( k+1, k-1 ) = zero
471  END IF
472  v2 = v( 2 )
473  t2 = dble( t1*v2 )
474 *
475 * Apply G from the left to transform the rows of the matrix
476 * in columns K to I2.
477 *
478  DO 80 j = k, i2
479  sum = dconjg( t1 )*h( k, j ) + t2*h( k+1, j )
480  h( k, j ) = h( k, j ) - sum
481  h( k+1, j ) = h( k+1, j ) - sum*v2
482  80 continue
483 *
484 * Apply G from the right to transform the columns of the
485 * matrix in rows I1 to min(K+2,I).
486 *
487  DO 90 j = i1, min( k+2, i )
488  sum = t1*h( j, k ) + t2*h( j, k+1 )
489  h( j, k ) = h( j, k ) - sum
490  h( j, k+1 ) = h( j, k+1 ) - sum*dconjg( v2 )
491  90 continue
492 *
493  IF( wantz ) THEN
494 *
495 * Accumulate transformations in the matrix Z
496 *
497  DO 100 j = iloz, ihiz
498  sum = t1*z( j, k ) + t2*z( j, k+1 )
499  z( j, k ) = z( j, k ) - sum
500  z( j, k+1 ) = z( j, k+1 ) - sum*dconjg( v2 )
501  100 continue
502  END IF
503 *
504  IF( k.EQ.m .AND. m.GT.l ) THEN
505 *
506 * If the QR step was started at row M > L because two
507 * consecutive small subdiagonals were found, then extra
508 * scaling must be performed to ensure that H(M,M-1) remains
509 * real.
510 *
511  temp = one - t1
512  temp = temp / abs( temp )
513  h( m+1, m ) = h( m+1, m )*dconjg( temp )
514  IF( m+2.LE.i )
515  $ h( m+2, m+1 ) = h( m+2, m+1 )*temp
516  DO 110 j = m, i
517  IF( j.NE.m+1 ) THEN
518  IF( i2.GT.j )
519  $ CALL zscal( i2-j, temp, h( j, j+1 ), ldh )
520  CALL zscal( j-i1, dconjg( temp ), h( i1, j ), 1 )
521  IF( wantz ) THEN
522  CALL zscal( nz, dconjg( temp ), z( iloz, j ),
523  $ 1 )
524  END IF
525  END IF
526  110 continue
527  END IF
528  120 continue
529 *
530 * Ensure that H(I,I-1) is real.
531 *
532  temp = h( i, i-1 )
533  IF( dimag( temp ).NE.rzero ) THEN
534  rtemp = abs( temp )
535  h( i, i-1 ) = rtemp
536  temp = temp / rtemp
537  IF( i2.GT.i )
538  $ CALL zscal( i2-i, dconjg( temp ), h( i, i+1 ), ldh )
539  CALL zscal( i-i1, temp, h( i1, i ), 1 )
540  IF( wantz ) THEN
541  CALL zscal( nz, temp, z( iloz, i ), 1 )
542  END IF
543  END IF
544 *
545  130 continue
546 *
547 * Failure to converge in remaining number of iterations
548 *
549  info = i
550  return
551 *
552  140 continue
553 *
554 * H(I,I-1) is negligible: one eigenvalue has converged.
555 *
556  w( i ) = h( i, i )
557 *
558 * return to start of the main loop with new value of I.
559 *
560  i = l - 1
561  go to 30
562 *
563  150 continue
564  return
565 *
566 * End of ZLAHQR
567 *
568  END