LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *> > 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 > 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 through INFO of the final, output
164 *> value of H.
165 *>
166 *> If INFO > 0 and WANTT is .TRUE., then on exit
167 *> (*) (initial value of H)*U = U*(final value of H)
168 *> where U is an orthogonal 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 > 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 *> \ingroup doubleOTHERauxiliary
187 *
188 *> \par Further Details:
189 * =====================
190 *>
191 *> \verbatim
192 *>
193 *> 02-96 Based on modifications by
194 *> David Day, Sandia National Laboratory, USA
195 *>
196 *> 12-04 Further modifications by
197 *> Ralph Byers, University of Kansas, USA
198 *> This is a modified version of DLAHQR from LAPACK version 3.0.
199 *> It is (1) more robust against overflow and underflow and
200 *> (2) adopts the more conservative Ahues & Tisseur stopping
201 *> criterion (LAWN 122, 1997).
202 *> \endverbatim
203 *>
204 * =====================================================================
205  SUBROUTINE dlahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
206  $ ILOZ, IHIZ, Z, LDZ, INFO )
207  IMPLICIT NONE
208 *
209 * -- LAPACK auxiliary routine --
210 * -- LAPACK is a software package provided by Univ. of Tennessee, --
211 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
212 *
213 * .. Scalar Arguments ..
214  INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
215  LOGICAL WANTT, WANTZ
216 * ..
217 * .. Array Arguments ..
218  DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
219 * ..
220 *
221 * =========================================================
222 *
223 * .. Parameters ..
224  DOUBLE PRECISION ZERO, ONE, TWO
225  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
226  DOUBLE PRECISION DAT1, DAT2
227  parameter( dat1 = 3.0d0 / 4.0d0, dat2 = -0.4375d0 )
228  INTEGER KEXSH
229  parameter( kexsh = 10 )
230 * ..
231 * .. Local Scalars ..
232  DOUBLE PRECISION AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
233  $ h22, rt1i, rt1r, rt2i, rt2r, rtdisc, s, safmax,
234  $ safmin, smlnum, sn, sum, t1, t2, t3, tr, tst,
235  $ ulp, v2, v3
236  INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ,
237  $ kdefl
238 * ..
239 * .. Local Arrays ..
240  DOUBLE PRECISION V( 3 )
241 * ..
242 * .. External Functions ..
243  DOUBLE PRECISION DLAMCH
244  EXTERNAL dlamch
245 * ..
246 * .. External Subroutines ..
247  EXTERNAL dcopy, dlabad, dlanv2, dlarfg, drot
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, dble, max, min, sqrt
251 * ..
252 * .. Executable Statements ..
253 *
254  info = 0
255 *
256 * Quick return if possible
257 *
258  IF( n.EQ.0 )
259  $ RETURN
260  IF( ilo.EQ.ihi ) THEN
261  wr( ilo ) = h( ilo, ilo )
262  wi( ilo ) = zero
263  RETURN
264  END IF
265 *
266 * ==== clear out the trash ====
267  DO 10 j = ilo, ihi - 3
268  h( j+2, j ) = zero
269  h( j+3, j ) = zero
270  10 CONTINUE
271  IF( ilo.LE.ihi-2 )
272  $ h( ihi, ihi-2 ) = zero
273 *
274  nh = ihi - ilo + 1
275  nz = ihiz - iloz + 1
276 *
277 * Set machine-dependent constants for the stopping criterion.
278 *
279  safmin = dlamch( 'SAFE MINIMUM' )
280  safmax = one / safmin
281  CALL dlabad( safmin, safmax )
282  ulp = dlamch( 'PRECISION' )
283  smlnum = safmin*( dble( nh ) / ulp )
284 *
285 * I1 and I2 are the indices of the first row and last column of H
286 * to which transformations must be applied. If eigenvalues only are
287 * being computed, I1 and I2 are set inside the main loop.
288 *
289  IF( wantt ) THEN
290  i1 = 1
291  i2 = n
292  END IF
293 *
294 * ITMAX is the total number of QR iterations allowed.
295 *
296  itmax = 30 * max( 10, nh )
297 *
298 * KDEFL counts the number of iterations since a deflation
299 *
300  kdefl = 0
301 *
302 * The main loop begins here. I is the loop index and decreases from
303 * IHI to ILO in steps of 1 or 2. Each iteration of the loop works
304 * with the active submatrix in rows and columns L to I.
305 * Eigenvalues I+1 to IHI have already converged. Either L = ILO or
306 * H(L,L-1) is negligible so that the matrix splits.
307 *
308  i = ihi
309  20 CONTINUE
310  l = ilo
311  IF( i.LT.ilo )
312  $ GO TO 160
313 *
314 * Perform QR iterations on rows and columns ILO to I until a
315 * submatrix of order 1 or 2 splits off at the bottom because a
316 * subdiagonal element has become negligible.
317 *
318  DO 140 its = 0, itmax
319 *
320 * Look for a single small subdiagonal element.
321 *
322  DO 30 k = i, l + 1, -1
323  IF( abs( h( k, k-1 ) ).LE.smlnum )
324  $ GO TO 40
325  tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
326  IF( tst.EQ.zero ) THEN
327  IF( k-2.GE.ilo )
328  $ tst = tst + abs( h( k-1, k-2 ) )
329  IF( k+1.LE.ihi )
330  $ tst = tst + abs( h( k+1, k ) )
331  END IF
332 * ==== The following is a conservative small subdiagonal
333 * . deflation criterion due to Ahues & Tisseur (LAWN 122,
334 * . 1997). It has better mathematical foundation and
335 * . improves accuracy in some cases. ====
336  IF( abs( h( k, k-1 ) ).LE.ulp*tst ) THEN
337  ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
338  ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
339  aa = max( abs( h( k, k ) ),
340  $ abs( h( k-1, k-1 )-h( k, k ) ) )
341  bb = min( abs( h( k, k ) ),
342  $ abs( h( k-1, k-1 )-h( k, k ) ) )
343  s = aa + ab
344  IF( ba*( ab / s ).LE.max( smlnum,
345  $ ulp*( bb*( aa / s ) ) ) )GO TO 40
346  END IF
347  30 CONTINUE
348  40 CONTINUE
349  l = k
350  IF( l.GT.ilo ) THEN
351 *
352 * H(L,L-1) is negligible
353 *
354  h( l, l-1 ) = zero
355  END IF
356 *
357 * Exit from loop if a submatrix of order 1 or 2 has split off.
358 *
359  IF( l.GE.i-1 )
360  $ GO TO 150
361  kdefl = kdefl + 1
362 *
363 * Now the active submatrix is in rows and columns L to I. If
364 * eigenvalues only are being computed, only the active submatrix
365 * need be transformed.
366 *
367  IF( .NOT.wantt ) THEN
368  i1 = l
369  i2 = i
370  END IF
371 *
372  IF( mod(kdefl,2*kexsh).EQ.0 ) THEN
373 *
374 * Exceptional shift.
375 *
376  s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
377  h11 = dat1*s + h( i, i )
378  h12 = dat2*s
379  h21 = s
380  h22 = h11
381  ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
382 *
383 * Exceptional shift.
384 *
385  s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
386  h11 = dat1*s + h( l, l )
387  h12 = dat2*s
388  h21 = s
389  h22 = h11
390  ELSE
391 *
392 * Prepare to use Francis' double shift
393 * (i.e. 2nd degree generalized Rayleigh quotient)
394 *
395  h11 = h( i-1, i-1 )
396  h21 = h( i, i-1 )
397  h12 = h( i-1, i )
398  h22 = h( i, i )
399  END IF
400  s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
401  IF( s.EQ.zero ) THEN
402  rt1r = zero
403  rt1i = zero
404  rt2r = zero
405  rt2i = zero
406  ELSE
407  h11 = h11 / s
408  h21 = h21 / s
409  h12 = h12 / s
410  h22 = h22 / s
411  tr = ( h11+h22 ) / two
412  det = ( h11-tr )*( h22-tr ) - h12*h21
413  rtdisc = sqrt( abs( det ) )
414  IF( det.GE.zero ) THEN
415 *
416 * ==== complex conjugate shifts ====
417 *
418  rt1r = tr*s
419  rt2r = rt1r
420  rt1i = rtdisc*s
421  rt2i = -rt1i
422  ELSE
423 *
424 * ==== real shifts (use only one of them) ====
425 *
426  rt1r = tr + rtdisc
427  rt2r = tr - rtdisc
428  IF( abs( rt1r-h22 ).LE.abs( rt2r-h22 ) ) THEN
429  rt1r = rt1r*s
430  rt2r = rt1r
431  ELSE
432  rt2r = rt2r*s
433  rt1r = rt2r
434  END IF
435  rt1i = zero
436  rt2i = zero
437  END IF
438  END IF
439 *
440 * Look for two consecutive small subdiagonal elements.
441 *
442  DO 50 m = i - 2, l, -1
443 * Determine the effect of starting the double-shift QR
444 * iteration at row M, and see if this would make H(M,M-1)
445 * negligible. (The following uses scaling to avoid
446 * overflows and most underflows.)
447 *
448  h21s = h( m+1, m )
449  s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
450  h21s = h( m+1, m ) / s
451  v( 1 ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*
452  $ ( ( h( m, m )-rt2r ) / s ) - rt1i*( rt2i / s )
453  v( 2 ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
454  v( 3 ) = h21s*h( m+2, m+1 )
455  s = abs( v( 1 ) ) + abs( v( 2 ) ) + abs( v( 3 ) )
456  v( 1 ) = v( 1 ) / s
457  v( 2 ) = v( 2 ) / s
458  v( 3 ) = v( 3 ) / s
459  IF( m.EQ.l )
460  $ GO TO 60
461  IF( abs( h( m, m-1 ) )*( abs( v( 2 ) )+abs( v( 3 ) ) ).LE.
462  $ ulp*abs( v( 1 ) )*( abs( h( m-1, m-1 ) )+abs( h( m,
463  $ m ) )+abs( h( m+1, m+1 ) ) ) )GO TO 60
464  50 CONTINUE
465  60 CONTINUE
466 *
467 * Double-shift QR step
468 *
469  DO 130 k = m, i - 1
470 *
471 * The first iteration of this loop determines a reflection G
472 * from the vector V and applies it from left and right to H,
473 * thus creating a nonzero bulge below the subdiagonal.
474 *
475 * Each subsequent iteration determines a reflection G to
476 * restore the Hessenberg form in the (K-1)th column, and thus
477 * chases the bulge one step toward the bottom of the active
478 * submatrix. NR is the order of G.
479 *
480  nr = min( 3, i-k+1 )
481  IF( k.GT.m )
482  $ CALL dcopy( nr, h( k, k-1 ), 1, v, 1 )
483  CALL dlarfg( nr, v( 1 ), v( 2 ), 1, t1 )
484  IF( k.GT.m ) THEN
485  h( k, k-1 ) = v( 1 )
486  h( k+1, k-1 ) = zero
487  IF( k.LT.i-1 )
488  $ h( k+2, k-1 ) = zero
489  ELSE IF( m.GT.l ) THEN
490 * ==== Use the following instead of
491 * . H( K, K-1 ) = -H( K, K-1 ) to
492 * . avoid a bug when v(2) and v(3)
493 * . underflow. ====
494  h( k, k-1 ) = h( k, k-1 )*( one-t1 )
495  END IF
496  v2 = v( 2 )
497  t2 = t1*v2
498  IF( nr.EQ.3 ) THEN
499  v3 = v( 3 )
500  t3 = t1*v3
501 *
502 * Apply G from the left to transform the rows of the matrix
503 * in columns K to I2.
504 *
505  DO 70 j = k, i2
506  sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
507  h( k, j ) = h( k, j ) - sum*t1
508  h( k+1, j ) = h( k+1, j ) - sum*t2
509  h( k+2, j ) = h( k+2, j ) - sum*t3
510  70 CONTINUE
511 *
512 * Apply G from the right to transform the columns of the
513 * matrix in rows I1 to min(K+3,I).
514 *
515  DO 80 j = i1, min( k+3, i )
516  sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
517  h( j, k ) = h( j, k ) - sum*t1
518  h( j, k+1 ) = h( j, k+1 ) - sum*t2
519  h( j, k+2 ) = h( j, k+2 ) - sum*t3
520  80 CONTINUE
521 *
522  IF( wantz ) THEN
523 *
524 * Accumulate transformations in the matrix Z
525 *
526  DO 90 j = iloz, ihiz
527  sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
528  z( j, k ) = z( j, k ) - sum*t1
529  z( j, k+1 ) = z( j, k+1 ) - sum*t2
530  z( j, k+2 ) = z( j, k+2 ) - sum*t3
531  90 CONTINUE
532  END IF
533  ELSE IF( nr.EQ.2 ) THEN
534 *
535 * Apply G from the left to transform the rows of the matrix
536 * in columns K to I2.
537 *
538  DO 100 j = k, i2
539  sum = h( k, j ) + v2*h( k+1, j )
540  h( k, j ) = h( k, j ) - sum*t1
541  h( k+1, j ) = h( k+1, j ) - sum*t2
542  100 CONTINUE
543 *
544 * Apply G from the right to transform the columns of the
545 * matrix in rows I1 to min(K+3,I).
546 *
547  DO 110 j = i1, i
548  sum = h( j, k ) + v2*h( j, k+1 )
549  h( j, k ) = h( j, k ) - sum*t1
550  h( j, k+1 ) = h( j, k+1 ) - sum*t2
551  110 CONTINUE
552 *
553  IF( wantz ) THEN
554 *
555 * Accumulate transformations in the matrix Z
556 *
557  DO 120 j = iloz, ihiz
558  sum = z( j, k ) + v2*z( j, k+1 )
559  z( j, k ) = z( j, k ) - sum*t1
560  z( j, k+1 ) = z( j, k+1 ) - sum*t2
561  120 CONTINUE
562  END IF
563  END IF
564  130 CONTINUE
565 *
566  140 CONTINUE
567 *
568 * Failure to converge in remaining number of iterations
569 *
570  info = i
571  RETURN
572 *
573  150 CONTINUE
574 *
575  IF( l.EQ.i ) THEN
576 *
577 * H(I,I-1) is negligible: one eigenvalue has converged.
578 *
579  wr( i ) = h( i, i )
580  wi( i ) = zero
581  ELSE IF( l.EQ.i-1 ) THEN
582 *
583 * H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
584 *
585 * Transform the 2-by-2 submatrix to standard Schur form,
586 * and compute and store the eigenvalues.
587 *
588  CALL dlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
589  $ h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
590  $ cs, sn )
591 *
592  IF( wantt ) THEN
593 *
594 * Apply the transformation to the rest of H.
595 *
596  IF( i2.GT.i )
597  $ CALL drot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
598  $ cs, sn )
599  CALL drot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs, sn )
600  END IF
601  IF( wantz ) THEN
602 *
603 * Apply the transformation to Z.
604 *
605  CALL drot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs, sn )
606  END IF
607  END IF
608 * reset deflation counter
609  kdefl = 0
610 *
611 * return to start of the main loop with new value of I.
612 *
613  i = l - 1
614  GO TO 20
615 *
616  160 CONTINUE
617  RETURN
618 *
619 * End of DLAHQR
620 *
621  END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dlanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition: dlanv2.f:127
subroutine dlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition: dlahqr.f:207
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:106