LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
clahqr.f
Go to the documentation of this file.
1 *> \brief \b CLAHQR 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 CLAHQR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clahqr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clahqr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clahqr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLAHQR( 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 H( LDH, * ), W( * ), Z( LDZ, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> CLAHQR 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 *> CLAHQR 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 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 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 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 *> > 0: if INFO = i, CLAHQR 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 > 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 through INFO of the final,
152 *> output value of H.
153 *>
154 *> If INFO > 0 and WANTT is .TRUE., then on exit
155 *> (*) (initial value of H)*U = U*(final value of H)
156 *> where U is an orthogonal 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 > 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 *> \ingroup complexOTHERauxiliary
175 *
176 *> \par Contributors:
177 * ==================
178 *>
179 *> \verbatim
180 *>
181 *> 02-96 Based on modifications by
182 *> David Day, Sandia National Laboratory, USA
183 *>
184 *> 12-04 Further modifications by
185 *> Ralph Byers, University of Kansas, USA
186 *> This is a modified version of CLAHQR from LAPACK version 3.0.
187 *> It is (1) more robust against overflow and underflow and
188 *> (2) adopts the more conservative Ahues & Tisseur stopping
189 *> criterion (LAWN 122, 1997).
190 *> \endverbatim
191 *>
192 * =====================================================================
193  SUBROUTINE clahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
194  $ IHIZ, Z, LDZ, INFO )
195  IMPLICIT NONE
196 *
197 * -- LAPACK auxiliary routine --
198 * -- LAPACK is a software package provided by Univ. of Tennessee, --
199 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
200 *
201 * .. Scalar Arguments ..
202  INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
203  LOGICAL WANTT, WANTZ
204 * ..
205 * .. Array Arguments ..
206  COMPLEX H( LDH, * ), W( * ), Z( LDZ, * )
207 * ..
208 *
209 * =========================================================
210 *
211 * .. Parameters ..
212  COMPLEX ZERO, ONE
213  parameter( zero = ( 0.0e0, 0.0e0 ),
214  $ one = ( 1.0e0, 0.0e0 ) )
215  REAL RZERO, RONE, HALF
216  parameter( rzero = 0.0e0, rone = 1.0e0, half = 0.5e0 )
217  REAL DAT1
218  parameter( dat1 = 3.0e0 / 4.0e0 )
219  INTEGER KEXSH
220  parameter( kexsh = 10 )
221 * ..
222 * .. Local Scalars ..
223  COMPLEX CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
224  $ v2, x, y
225  REAL AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
226  $ safmin, smlnum, sx, t2, tst, ulp
227  INTEGER I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
228  $ nh, nz, kdefl
229 * ..
230 * .. Local Arrays ..
231  COMPLEX V( 2 )
232 * ..
233 * .. External Functions ..
234  COMPLEX CLADIV
235  REAL SLAMCH
236  EXTERNAL cladiv, slamch
237 * ..
238 * .. External Subroutines ..
239  EXTERNAL ccopy, clarfg, cscal, slabad
240 * ..
241 * .. Statement Functions ..
242  REAL CABS1
243 * ..
244 * .. Intrinsic Functions ..
245  INTRINSIC abs, aimag, conjg, max, min, real, sqrt
246 * ..
247 * .. Statement Function definitions ..
248  cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
249 * ..
250 * .. Executable Statements ..
251 *
252  info = 0
253 *
254 * Quick return if possible
255 *
256  IF( n.EQ.0 )
257  $ RETURN
258  IF( ilo.EQ.ihi ) THEN
259  w( ilo ) = h( ilo, ilo )
260  RETURN
261  END IF
262 *
263 * ==== clear out the trash ====
264  DO 10 j = ilo, ihi - 3
265  h( j+2, j ) = zero
266  h( j+3, j ) = zero
267  10 CONTINUE
268  IF( ilo.LE.ihi-2 )
269  $ h( ihi, ihi-2 ) = zero
270 * ==== ensure that subdiagonal entries are real ====
271  IF( wantt ) THEN
272  jlo = 1
273  jhi = n
274  ELSE
275  jlo = ilo
276  jhi = ihi
277  END IF
278  DO 20 i = ilo + 1, ihi
279  IF( aimag( h( i, i-1 ) ).NE.rzero ) THEN
280 * ==== The following redundant normalization
281 * . avoids problems with both gradual and
282 * . sudden underflow in ABS(H(I,I-1)) ====
283  sc = h( i, i-1 ) / cabs1( h( i, i-1 ) )
284  sc = conjg( sc ) / abs( sc )
285  h( i, i-1 ) = abs( h( i, i-1 ) )
286  CALL cscal( jhi-i+1, sc, h( i, i ), ldh )
287  CALL cscal( min( jhi, i+1 )-jlo+1, conjg( sc ), h( jlo, i ),
288  $ 1 )
289  IF( wantz )
290  $ CALL cscal( ihiz-iloz+1, conjg( sc ), z( iloz, i ), 1 )
291  END IF
292  20 CONTINUE
293 *
294  nh = ihi - ilo + 1
295  nz = ihiz - iloz + 1
296 *
297 * Set machine-dependent constants for the stopping criterion.
298 *
299  safmin = slamch( 'SAFE MINIMUM' )
300  safmax = rone / safmin
301  CALL slabad( safmin, safmax )
302  ulp = slamch( 'PRECISION' )
303  smlnum = safmin*( real( nh ) / ulp )
304 *
305 * I1 and I2 are the indices of the first row and last column of H
306 * to which transformations must be applied. If eigenvalues only are
307 * being computed, I1 and I2 are set inside the main loop.
308 *
309  IF( wantt ) THEN
310  i1 = 1
311  i2 = n
312  END IF
313 *
314 * ITMAX is the total number of QR iterations allowed.
315 *
316  itmax = 30 * max( 10, nh )
317 *
318 * KDEFL counts the number of iterations since a deflation
319 *
320  kdefl = 0
321 *
322 * The main loop begins here. I is the loop index and decreases from
323 * IHI to ILO in steps of 1. Each iteration of the loop works
324 * with the active submatrix in rows and columns L to I.
325 * Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
326 * H(L,L-1) is negligible so that the matrix splits.
327 *
328  i = ihi
329  30 CONTINUE
330  IF( i.LT.ilo )
331  $ GO TO 150
332 *
333 * Perform QR iterations on rows and columns ILO to I until a
334 * submatrix of order 1 splits off at the bottom because a
335 * subdiagonal element has become negligible.
336 *
337  l = ilo
338  DO 130 its = 0, itmax
339 *
340 * Look for a single small subdiagonal element.
341 *
342  DO 40 k = i, l + 1, -1
343  IF( cabs1( h( k, k-1 ) ).LE.smlnum )
344  $ GO TO 50
345  tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
346  IF( tst.EQ.zero ) THEN
347  IF( k-2.GE.ilo )
348  $ tst = tst + abs( real( h( k-1, k-2 ) ) )
349  IF( k+1.LE.ihi )
350  $ tst = tst + abs( real( h( k+1, k ) ) )
351  END IF
352 * ==== The following is a conservative small subdiagonal
353 * . deflation criterion due to Ahues & Tisseur (LAWN 122,
354 * . 1997). It has better mathematical foundation and
355 * . improves accuracy in some examples. ====
356  IF( abs( real( h( k, k-1 ) ) ).LE.ulp*tst ) THEN
357  ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
358  ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
359  aa = max( cabs1( h( k, k ) ),
360  $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
361  bb = min( cabs1( h( k, k ) ),
362  $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
363  s = aa + ab
364  IF( ba*( ab / s ).LE.max( smlnum,
365  $ ulp*( bb*( aa / s ) ) ) )GO TO 50
366  END IF
367  40 CONTINUE
368  50 CONTINUE
369  l = k
370  IF( l.GT.ilo ) THEN
371 *
372 * H(L,L-1) is negligible
373 *
374  h( l, l-1 ) = zero
375  END IF
376 *
377 * Exit from loop if a submatrix of order 1 has split off.
378 *
379  IF( l.GE.i )
380  $ GO TO 140
381  kdefl = kdefl + 1
382 *
383 * Now the active submatrix is in rows and columns L to I. If
384 * eigenvalues only are being computed, only the active submatrix
385 * need be transformed.
386 *
387  IF( .NOT.wantt ) THEN
388  i1 = l
389  i2 = i
390  END IF
391 *
392  IF( mod(kdefl,2*kexsh).EQ.0 ) THEN
393 *
394 * Exceptional shift.
395 *
396  s = dat1*abs( real( h( i, i-1 ) ) )
397  t = s + h( i, i )
398  ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
399 *
400 * Exceptional shift.
401 *
402  s = dat1*abs( real( h( l+1, l ) ) )
403  t = s + h( l, l )
404  ELSE
405 *
406 * Wilkinson's shift.
407 *
408  t = h( i, i )
409  u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
410  s = cabs1( u )
411  IF( s.NE.rzero ) THEN
412  x = half*( h( i-1, i-1 )-t )
413  sx = cabs1( x )
414  s = max( s, cabs1( x ) )
415  y = s*sqrt( ( x / s )**2+( u / s )**2 )
416  IF( sx.GT.rzero ) THEN
417  IF( real( x / sx )*real( y )+aimag( x / sx )*
418  $ aimag( y ).LT.rzero )y = -y
419  END IF
420  t = t - u*cladiv( u, ( x+y ) )
421  END IF
422  END IF
423 *
424 * Look for two consecutive small subdiagonal elements.
425 *
426  DO 60 m = i - 1, l + 1, -1
427 *
428 * Determine the effect of starting the single-shift QR
429 * iteration at row M, and see if this would make H(M,M-1)
430 * negligible.
431 *
432  h11 = h( m, m )
433  h22 = h( m+1, m+1 )
434  h11s = h11 - t
435  h21 = real( h( m+1, m ) )
436  s = cabs1( h11s ) + abs( h21 )
437  h11s = h11s / s
438  h21 = h21 / s
439  v( 1 ) = h11s
440  v( 2 ) = h21
441  h10 = real( h( m, m-1 ) )
442  IF( abs( h10 )*abs( h21 ).LE.ulp*
443  $ ( cabs1( h11s )*( cabs1( h11 )+cabs1( h22 ) ) ) )
444  $ GO TO 70
445  60 CONTINUE
446  h11 = h( l, l )
447  h22 = h( l+1, l+1 )
448  h11s = h11 - t
449  h21 = real( h( l+1, l ) )
450  s = cabs1( h11s ) + abs( h21 )
451  h11s = h11s / s
452  h21 = h21 / s
453  v( 1 ) = h11s
454  v( 2 ) = h21
455  70 CONTINUE
456 *
457 * Single-shift QR step
458 *
459  DO 120 k = m, i - 1
460 *
461 * The first iteration of this loop determines a reflection G
462 * from the vector V and applies it from left and right to H,
463 * thus creating a nonzero bulge below the subdiagonal.
464 *
465 * Each subsequent iteration determines a reflection G to
466 * restore the Hessenberg form in the (K-1)th column, and thus
467 * chases the bulge one step toward the bottom of the active
468 * submatrix.
469 *
470 * V(2) is always real before the call to CLARFG, and hence
471 * after the call T2 ( = T1*V(2) ) is also real.
472 *
473  IF( k.GT.m )
474  $ CALL ccopy( 2, h( k, k-1 ), 1, v, 1 )
475  CALL clarfg( 2, 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  END IF
480  v2 = v( 2 )
481  t2 = real( t1*v2 )
482 *
483 * Apply G from the left to transform the rows of the matrix
484 * in columns K to I2.
485 *
486  DO 80 j = k, i2
487  sum = conjg( t1 )*h( k, j ) + t2*h( k+1, j )
488  h( k, j ) = h( k, j ) - sum
489  h( k+1, j ) = h( k+1, j ) - sum*v2
490  80 CONTINUE
491 *
492 * Apply G from the right to transform the columns of the
493 * matrix in rows I1 to min(K+2,I).
494 *
495  DO 90 j = i1, min( k+2, i )
496  sum = t1*h( j, k ) + t2*h( j, k+1 )
497  h( j, k ) = h( j, k ) - sum
498  h( j, k+1 ) = h( j, k+1 ) - sum*conjg( v2 )
499  90 CONTINUE
500 *
501  IF( wantz ) THEN
502 *
503 * Accumulate transformations in the matrix Z
504 *
505  DO 100 j = iloz, ihiz
506  sum = t1*z( j, k ) + t2*z( j, k+1 )
507  z( j, k ) = z( j, k ) - sum
508  z( j, k+1 ) = z( j, k+1 ) - sum*conjg( v2 )
509  100 CONTINUE
510  END IF
511 *
512  IF( k.EQ.m .AND. m.GT.l ) THEN
513 *
514 * If the QR step was started at row M > L because two
515 * consecutive small subdiagonals were found, then extra
516 * scaling must be performed to ensure that H(M,M-1) remains
517 * real.
518 *
519  temp = one - t1
520  temp = temp / abs( temp )
521  h( m+1, m ) = h( m+1, m )*conjg( temp )
522  IF( m+2.LE.i )
523  $ h( m+2, m+1 ) = h( m+2, m+1 )*temp
524  DO 110 j = m, i
525  IF( j.NE.m+1 ) THEN
526  IF( i2.GT.j )
527  $ CALL cscal( i2-j, temp, h( j, j+1 ), ldh )
528  CALL cscal( j-i1, conjg( temp ), h( i1, j ), 1 )
529  IF( wantz ) THEN
530  CALL cscal( nz, conjg( temp ), z( iloz, j ), 1 )
531  END IF
532  END IF
533  110 CONTINUE
534  END IF
535  120 CONTINUE
536 *
537 * Ensure that H(I,I-1) is real.
538 *
539  temp = h( i, i-1 )
540  IF( aimag( temp ).NE.rzero ) THEN
541  rtemp = abs( temp )
542  h( i, i-1 ) = rtemp
543  temp = temp / rtemp
544  IF( i2.GT.i )
545  $ CALL cscal( i2-i, conjg( temp ), h( i, i+1 ), ldh )
546  CALL cscal( i-i1, temp, h( i1, i ), 1 )
547  IF( wantz ) THEN
548  CALL cscal( nz, temp, z( iloz, i ), 1 )
549  END IF
550  END IF
551 *
552  130 CONTINUE
553 *
554 * Failure to converge in remaining number of iterations
555 *
556  info = i
557  RETURN
558 *
559  140 CONTINUE
560 *
561 * H(I,I-1) is negligible: one eigenvalue has converged.
562 *
563  w( i ) = h( i, i )
564 * reset deflation counter
565  kdefl = 0
566 *
567 * return to start of the main loop with new value of I.
568 *
569  i = l - 1
570  GO TO 30
571 *
572  150 CONTINUE
573  RETURN
574 *
575 * End of CLAHQR
576 *
577  END
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:106
subroutine clahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition: clahqr.f:195