LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlahqr.f
Go to the documentation of this file.
1 *> \brief \b DLAHQR 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 DLAHQR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahqr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahqr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahqr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
22 * ILOZ, 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 * DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> DLAHQR is an auxiliary routine called by DHSEQR to update the
39 *> eigenvalues and Schur decomposition already computed by DHSEQR, 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 quasi-triangular in
76 *> rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
77 *> ILO = 1). DLAHQR works primarily with the Hessenberg
78 *> submatrix in rows and columns ILO to IHI, but applies
79 *> transformations to all of 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 DOUBLE PRECISION array, dimension (LDH,N)
86 *> On entry, the upper Hessenberg matrix H.
87 *> On exit, if INFO is zero and if WANTT is .TRUE., H is upper
88 *> quasi-triangular in rows and columns ILO:IHI, with any
89 *> 2-by-2 diagonal blocks in standard form. If INFO is zero
90 *> and WANTT is .FALSE., the contents of H are unspecified on
91 *> exit. The output state of H if INFO is nonzero is given
92 *> below under the description of INFO.
93 *> \endverbatim
94 *>
95 *> \param[in] LDH
96 *> \verbatim
97 *> LDH is INTEGER
98 *> The leading dimension of the array H. LDH >= max(1,N).
99 *> \endverbatim
100 *>
101 *> \param[out] WR
102 *> \verbatim
103 *> WR is DOUBLE PRECISION array, dimension (N)
104 *> \endverbatim
105 *>
106 *> \param[out] WI
107 *> \verbatim
108 *> WI is DOUBLE PRECISION array, dimension (N)
109 *> The real and imaginary parts, respectively, of the computed
110 *> eigenvalues ILO to IHI are stored in the corresponding
111 *> elements of WR and WI. If two eigenvalues are computed as a
112 *> complex conjugate pair, they are stored in consecutive
113 *> elements of WR and WI, say the i-th and (i+1)th, with
114 *> WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
115 *> eigenvalues are stored in the same order as on the diagonal
116 *> of the Schur form returned in H, with WR(i) = H(i,i), and, if
117 *> H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
118 *> WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
119 *> \endverbatim
120 *>
121 *> \param[in] ILOZ
122 *> \verbatim
123 *> ILOZ is INTEGER
124 *> \endverbatim
125 *>
126 *> \param[in] IHIZ
127 *> \verbatim
128 *> IHIZ is INTEGER
129 *> Specify the rows of Z to which transformations must be
130 *> applied if WANTZ is .TRUE..
131 *> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
132 *> \endverbatim
133 *>
134 *> \param[in,out] Z
135 *> \verbatim
136 *> Z is DOUBLE PRECISION array, dimension (LDZ,N)
137 *> If WANTZ is .TRUE., on entry Z must contain the current
138 *> matrix Z of transformations accumulated by DHSEQR, and on
139 *> exit Z has been updated; transformations are applied only to
140 *> the submatrix Z(ILOZ:IHIZ,ILO:IHI).
141 *> If WANTZ is .FALSE., Z is not referenced.
142 *> \endverbatim
143 *>
144 *> \param[in] LDZ
145 *> \verbatim
146 *> LDZ is INTEGER
147 *> The leading dimension of the array Z. LDZ >= max(1,N).
148 *> \endverbatim
149 *>
150 *> \param[out] INFO
151 *> \verbatim
152 *> INFO is INTEGER
153 *> = 0: successful exit
154 *> .GT. 0: If INFO = i, DLAHQR failed to compute all the
155 *> eigenvalues ILO to IHI in a total of 30 iterations
156 *> per eigenvalue; elements i+1:ihi of WR and WI
157 *> contain those eigenvalues which have been
158 *> successfully computed.
159 *>
160 *> If INFO .GT. 0 and WANTT is .FALSE., then on exit,
161 *> the remaining unconverged eigenvalues are the
162 *> eigenvalues of the upper Hessenberg matrix rows
163 *> and columns ILO thorugh INFO of the final, output
164 *> value of H.
165 *>
166 *> If INFO .GT. 0 and WANTT is .TRUE., then on exit
167 *> (*) (initial value of H)*U = U*(final value of H)
168 *> where U is an orthognal matrix. The final
169 *> value of H is upper Hessenberg and triangular in
170 *> rows and columns INFO+1 through IHI.
171 *>
172 *> If INFO .GT. 0 and WANTZ is .TRUE., then on exit
173 *> (final value of Z) = (initial value of Z)*U
174 *> where U is the orthogonal matrix in (*)
175 *> (regardless of the value of WANTT.)
176 *> \endverbatim
177 *
178 * Authors:
179 * ========
180 *
181 *> \author Univ. of Tennessee
182 *> \author Univ. of California Berkeley
183 *> \author Univ. of Colorado Denver
184 *> \author NAG Ltd.
185 *
186 *> \date September 2012
187 *
188 *> \ingroup doubleOTHERauxiliary
189 *
190 *> \par Further Details:
191 * =====================
192 *>
193 *> \verbatim
194 *>
195 *> 02-96 Based on modifications by
196 *> David Day, Sandia National Laboratory, USA
197 *>
198 *> 12-04 Further modifications by
199 *> Ralph Byers, University of Kansas, USA
200 *> This is a modified version of DLAHQR from LAPACK version 3.0.
201 *> It is (1) more robust against overflow and underflow and
202 *> (2) adopts the more conservative Ahues & Tisseur stopping
203 *> criterion (LAWN 122, 1997).
204 *> \endverbatim
205 *>
206 * =====================================================================
207  SUBROUTINE dlahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
208  $ iloz, ihiz, z, ldz, info )
209 *
210 * -- LAPACK auxiliary routine (version 3.4.2) --
211 * -- LAPACK is a software package provided by Univ. of Tennessee, --
212 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213 * September 2012
214 *
215 * .. Scalar Arguments ..
216  INTEGER ihi, ihiz, ilo, iloz, info, ldh, ldz, n
217  LOGICAL wantt, wantz
218 * ..
219 * .. Array Arguments ..
220  DOUBLE PRECISION h( ldh, * ), wi( * ), wr( * ), z( ldz, * )
221 * ..
222 *
223 * =========================================================
224 *
225 * .. Parameters ..
226  INTEGER itmax
227  parameter( itmax = 30 )
228  DOUBLE PRECISION zero, one, two
229  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
230  DOUBLE PRECISION dat1, dat2
231  parameter( dat1 = 3.0d0 / 4.0d0, dat2 = -0.4375d0 )
232 * ..
233 * .. Local Scalars ..
234  DOUBLE PRECISION aa, ab, ba, bb, cs, det, h11, h12, h21, h21s,
235  $ h22, rt1i, rt1r, rt2i, rt2r, rtdisc, s, safmax,
236  $ safmin, smlnum, sn, sum, t1, t2, t3, tr, tst,
237  $ ulp, v2, v3
238  INTEGER i, i1, i2, its, j, k, l, m, nh, nr, nz
239 * ..
240 * .. Local Arrays ..
241  DOUBLE PRECISION v( 3 )
242 * ..
243 * .. External Functions ..
244  DOUBLE PRECISION dlamch
245  EXTERNAL dlamch
246 * ..
247 * .. External Subroutines ..
248  EXTERNAL dcopy, dlabad, dlanv2, dlarfg, drot
249 * ..
250 * .. Intrinsic Functions ..
251  INTRINSIC abs, dble, max, min, sqrt
252 * ..
253 * .. Executable Statements ..
254 *
255  info = 0
256 *
257 * Quick return if possible
258 *
259  IF( n.EQ.0 )
260  $ return
261  IF( ilo.EQ.ihi ) THEN
262  wr( ilo ) = h( ilo, ilo )
263  wi( ilo ) = zero
264  return
265  END IF
266 *
267 * ==== clear out the trash ====
268  DO 10 j = ilo, ihi - 3
269  h( j+2, j ) = zero
270  h( j+3, j ) = zero
271  10 continue
272  IF( ilo.LE.ihi-2 )
273  $ h( ihi, ihi-2 ) = zero
274 *
275  nh = ihi - ilo + 1
276  nz = ihiz - iloz + 1
277 *
278 * Set machine-dependent constants for the stopping criterion.
279 *
280  safmin = dlamch( 'SAFE MINIMUM' )
281  safmax = one / safmin
282  CALL dlabad( safmin, safmax )
283  ulp = dlamch( 'PRECISION' )
284  smlnum = safmin*( dble( nh ) / ulp )
285 *
286 * I1 and I2 are the indices of the first row and last column of H
287 * to which transformations must be applied. If eigenvalues only are
288 * being computed, I1 and I2 are set inside the main loop.
289 *
290  IF( wantt ) THEN
291  i1 = 1
292  i2 = n
293  END IF
294 *
295 * The main loop begins here. I is the loop index and decreases from
296 * IHI to ILO in steps of 1 or 2. Each iteration of the loop works
297 * with the active submatrix in rows and columns L to I.
298 * Eigenvalues I+1 to IHI have already converged. Either L = ILO or
299 * H(L,L-1) is negligible so that the matrix splits.
300 *
301  i = ihi
302  20 continue
303  l = ilo
304  IF( i.LT.ilo )
305  $ go to 160
306 *
307 * Perform QR iterations on rows and columns ILO to I until a
308 * submatrix of order 1 or 2 splits off at the bottom because a
309 * subdiagonal element has become negligible.
310 *
311  DO 140 its = 0, itmax
312 *
313 * Look for a single small subdiagonal element.
314 *
315  DO 30 k = i, l + 1, -1
316  IF( abs( h( k, k-1 ) ).LE.smlnum )
317  $ go to 40
318  tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
319  IF( tst.EQ.zero ) THEN
320  IF( k-2.GE.ilo )
321  $ tst = tst + abs( h( k-1, k-2 ) )
322  IF( k+1.LE.ihi )
323  $ tst = tst + abs( h( k+1, k ) )
324  END IF
325 * ==== The following is a conservative small subdiagonal
326 * . deflation criterion due to Ahues & Tisseur (LAWN 122,
327 * . 1997). It has better mathematical foundation and
328 * . improves accuracy in some cases. ====
329  IF( abs( h( k, k-1 ) ).LE.ulp*tst ) THEN
330  ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
331  ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
332  aa = max( abs( h( k, k ) ),
333  $ abs( h( k-1, k-1 )-h( k, k ) ) )
334  bb = min( abs( h( k, k ) ),
335  $ abs( h( k-1, k-1 )-h( k, k ) ) )
336  s = aa + ab
337  IF( ba*( ab / s ).LE.max( smlnum,
338  $ ulp*( bb*( aa / s ) ) ) )go to 40
339  END IF
340  30 continue
341  40 continue
342  l = k
343  IF( l.GT.ilo ) THEN
344 *
345 * H(L,L-1) is negligible
346 *
347  h( l, l-1 ) = zero
348  END IF
349 *
350 * Exit from loop if a submatrix of order 1 or 2 has split off.
351 *
352  IF( l.GE.i-1 )
353  $ go to 150
354 *
355 * Now the active submatrix is in rows and columns L to I. If
356 * eigenvalues only are being computed, only the active submatrix
357 * need be transformed.
358 *
359  IF( .NOT.wantt ) THEN
360  i1 = l
361  i2 = i
362  END IF
363 *
364  IF( its.EQ.10 ) THEN
365 *
366 * Exceptional shift.
367 *
368  s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
369  h11 = dat1*s + h( l, l )
370  h12 = dat2*s
371  h21 = s
372  h22 = h11
373  ELSE IF( its.EQ.20 ) THEN
374 *
375 * Exceptional shift.
376 *
377  s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
378  h11 = dat1*s + h( i, i )
379  h12 = dat2*s
380  h21 = s
381  h22 = h11
382  ELSE
383 *
384 * Prepare to use Francis' double shift
385 * (i.e. 2nd degree generalized Rayleigh quotient)
386 *
387  h11 = h( i-1, i-1 )
388  h21 = h( i, i-1 )
389  h12 = h( i-1, i )
390  h22 = h( i, i )
391  END IF
392  s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
393  IF( s.EQ.zero ) THEN
394  rt1r = zero
395  rt1i = zero
396  rt2r = zero
397  rt2i = zero
398  ELSE
399  h11 = h11 / s
400  h21 = h21 / s
401  h12 = h12 / s
402  h22 = h22 / s
403  tr = ( h11+h22 ) / two
404  det = ( h11-tr )*( h22-tr ) - h12*h21
405  rtdisc = sqrt( abs( det ) )
406  IF( det.GE.zero ) THEN
407 *
408 * ==== complex conjugate shifts ====
409 *
410  rt1r = tr*s
411  rt2r = rt1r
412  rt1i = rtdisc*s
413  rt2i = -rt1i
414  ELSE
415 *
416 * ==== real shifts (use only one of them) ====
417 *
418  rt1r = tr + rtdisc
419  rt2r = tr - rtdisc
420  IF( abs( rt1r-h22 ).LE.abs( rt2r-h22 ) ) THEN
421  rt1r = rt1r*s
422  rt2r = rt1r
423  ELSE
424  rt2r = rt2r*s
425  rt1r = rt2r
426  END IF
427  rt1i = zero
428  rt2i = zero
429  END IF
430  END IF
431 *
432 * Look for two consecutive small subdiagonal elements.
433 *
434  DO 50 m = i - 2, l, -1
435 * Determine the effect of starting the double-shift QR
436 * iteration at row M, and see if this would make H(M,M-1)
437 * negligible. (The following uses scaling to avoid
438 * overflows and most underflows.)
439 *
440  h21s = h( m+1, m )
441  s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
442  h21s = h( m+1, m ) / s
443  v( 1 ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*
444  $ ( ( h( m, m )-rt2r ) / s ) - rt1i*( rt2i / s )
445  v( 2 ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
446  v( 3 ) = h21s*h( m+2, m+1 )
447  s = abs( v( 1 ) ) + abs( v( 2 ) ) + abs( v( 3 ) )
448  v( 1 ) = v( 1 ) / s
449  v( 2 ) = v( 2 ) / s
450  v( 3 ) = v( 3 ) / s
451  IF( m.EQ.l )
452  $ go to 60
453  IF( abs( h( m, m-1 ) )*( abs( v( 2 ) )+abs( v( 3 ) ) ).LE.
454  $ ulp*abs( v( 1 ) )*( abs( h( m-1, m-1 ) )+abs( h( m,
455  $ m ) )+abs( h( m+1, m+1 ) ) ) )go to 60
456  50 continue
457  60 continue
458 *
459 * Double-shift QR step
460 *
461  DO 130 k = m, i - 1
462 *
463 * The first iteration of this loop determines a reflection G
464 * from the vector V and applies it from left and right to H,
465 * thus creating a nonzero bulge below the subdiagonal.
466 *
467 * Each subsequent iteration determines a reflection G to
468 * restore the Hessenberg form in the (K-1)th column, and thus
469 * chases the bulge one step toward the bottom of the active
470 * submatrix. NR is the order of G.
471 *
472  nr = min( 3, i-k+1 )
473  IF( k.GT.m )
474  $ CALL dcopy( nr, h( k, k-1 ), 1, v, 1 )
475  CALL dlarfg( nr, v( 1 ), v( 2 ), 1, t1 )
476  IF( k.GT.m ) THEN
477  h( k, k-1 ) = v( 1 )
478  h( k+1, k-1 ) = zero
479  IF( k.LT.i-1 )
480  $ h( k+2, k-1 ) = zero
481  ELSE IF( m.GT.l ) THEN
482 * ==== Use the following instead of
483 * . H( K, K-1 ) = -H( K, K-1 ) to
484 * . avoid a bug when v(2) and v(3)
485 * . underflow. ====
486  h( k, k-1 ) = h( k, k-1 )*( one-t1 )
487  END IF
488  v2 = v( 2 )
489  t2 = t1*v2
490  IF( nr.EQ.3 ) THEN
491  v3 = v( 3 )
492  t3 = t1*v3
493 *
494 * Apply G from the left to transform the rows of the matrix
495 * in columns K to I2.
496 *
497  DO 70 j = k, i2
498  sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
499  h( k, j ) = h( k, j ) - sum*t1
500  h( k+1, j ) = h( k+1, j ) - sum*t2
501  h( k+2, j ) = h( k+2, j ) - sum*t3
502  70 continue
503 *
504 * Apply G from the right to transform the columns of the
505 * matrix in rows I1 to min(K+3,I).
506 *
507  DO 80 j = i1, min( k+3, i )
508  sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
509  h( j, k ) = h( j, k ) - sum*t1
510  h( j, k+1 ) = h( j, k+1 ) - sum*t2
511  h( j, k+2 ) = h( j, k+2 ) - sum*t3
512  80 continue
513 *
514  IF( wantz ) THEN
515 *
516 * Accumulate transformations in the matrix Z
517 *
518  DO 90 j = iloz, ihiz
519  sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
520  z( j, k ) = z( j, k ) - sum*t1
521  z( j, k+1 ) = z( j, k+1 ) - sum*t2
522  z( j, k+2 ) = z( j, k+2 ) - sum*t3
523  90 continue
524  END IF
525  ELSE IF( nr.EQ.2 ) THEN
526 *
527 * Apply G from the left to transform the rows of the matrix
528 * in columns K to I2.
529 *
530  DO 100 j = k, i2
531  sum = h( k, j ) + v2*h( k+1, j )
532  h( k, j ) = h( k, j ) - sum*t1
533  h( k+1, j ) = h( k+1, j ) - sum*t2
534  100 continue
535 *
536 * Apply G from the right to transform the columns of the
537 * matrix in rows I1 to min(K+3,I).
538 *
539  DO 110 j = i1, i
540  sum = h( j, k ) + v2*h( j, k+1 )
541  h( j, k ) = h( j, k ) - sum*t1
542  h( j, k+1 ) = h( j, k+1 ) - sum*t2
543  110 continue
544 *
545  IF( wantz ) THEN
546 *
547 * Accumulate transformations in the matrix Z
548 *
549  DO 120 j = iloz, ihiz
550  sum = z( j, k ) + v2*z( j, k+1 )
551  z( j, k ) = z( j, k ) - sum*t1
552  z( j, k+1 ) = z( j, k+1 ) - sum*t2
553  120 continue
554  END IF
555  END IF
556  130 continue
557 *
558  140 continue
559 *
560 * Failure to converge in remaining number of iterations
561 *
562  info = i
563  return
564 *
565  150 continue
566 *
567  IF( l.EQ.i ) THEN
568 *
569 * H(I,I-1) is negligible: one eigenvalue has converged.
570 *
571  wr( i ) = h( i, i )
572  wi( i ) = zero
573  ELSE IF( l.EQ.i-1 ) THEN
574 *
575 * H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
576 *
577 * Transform the 2-by-2 submatrix to standard Schur form,
578 * and compute and store the eigenvalues.
579 *
580  CALL dlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
581  $ h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
582  $ cs, sn )
583 *
584  IF( wantt ) THEN
585 *
586 * Apply the transformation to the rest of H.
587 *
588  IF( i2.GT.i )
589  $ CALL drot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
590  $ cs, sn )
591  CALL drot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs, sn )
592  END IF
593  IF( wantz ) THEN
594 *
595 * Apply the transformation to Z.
596 *
597  CALL drot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs, sn )
598  END IF
599  END IF
600 *
601 * return to start of the main loop with new value of I.
602 *
603  i = l - 1
604  go to 20
605 *
606  160 continue
607  return
608 *
609 * End of DLAHQR
610 *
611  END