LAPACK  3.9.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 *> \date December 2016
175 *
176 *> \ingroup complexOTHERauxiliary
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 CLAHQR 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 clahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
196  $ IHIZ, Z, LDZ, INFO )
197 *
198 * -- LAPACK auxiliary routine (version 3.7.0) --
199 * -- LAPACK is a software package provided by Univ. of Tennessee, --
200 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201 * December 2016
202 *
203 * .. Scalar Arguments ..
204  INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
205  LOGICAL WANTT, WANTZ
206 * ..
207 * .. Array Arguments ..
208  COMPLEX H( LDH, * ), W( * ), Z( LDZ, * )
209 * ..
210 *
211 * =========================================================
212 *
213 * .. Parameters ..
214  COMPLEX ZERO, ONE
215  parameter( zero = ( 0.0e0, 0.0e0 ),
216  $ one = ( 1.0e0, 0.0e0 ) )
217  REAL RZERO, RONE, HALF
218  parameter( rzero = 0.0e0, rone = 1.0e0, half = 0.5e0 )
219  REAL DAT1
220  parameter( dat1 = 3.0e0 / 4.0e0 )
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
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 * The main loop begins here. I is the loop index and decreases from
319 * IHI to ILO in steps of 1. Each iteration of the loop works
320 * with the active submatrix in rows and columns L to I.
321 * Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
322 * H(L,L-1) is negligible so that the matrix splits.
323 *
324  i = ihi
325  30 CONTINUE
326  IF( i.LT.ilo )
327  $ GO TO 150
328 *
329 * Perform QR iterations on rows and columns ILO to I until a
330 * submatrix of order 1 splits off at the bottom because a
331 * subdiagonal element has become negligible.
332 *
333  l = ilo
334  DO 130 its = 0, itmax
335 *
336 * Look for a single small subdiagonal element.
337 *
338  DO 40 k = i, l + 1, -1
339  IF( cabs1( h( k, k-1 ) ).LE.smlnum )
340  $ GO TO 50
341  tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
342  IF( tst.EQ.zero ) THEN
343  IF( k-2.GE.ilo )
344  $ tst = tst + abs( real( h( k-1, k-2 ) ) )
345  IF( k+1.LE.ihi )
346  $ tst = tst + abs( real( h( k+1, k ) ) )
347  END IF
348 * ==== The following is a conservative small subdiagonal
349 * . deflation criterion due to Ahues & Tisseur (LAWN 122,
350 * . 1997). It has better mathematical foundation and
351 * . improves accuracy in some examples. ====
352  IF( abs( real( h( k, k-1 ) ) ).LE.ulp*tst ) THEN
353  ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
354  ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
355  aa = max( cabs1( h( k, k ) ),
356  $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
357  bb = min( cabs1( h( k, k ) ),
358  $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
359  s = aa + ab
360  IF( ba*( ab / s ).LE.max( smlnum,
361  $ ulp*( bb*( aa / s ) ) ) )GO TO 50
362  END IF
363  40 CONTINUE
364  50 CONTINUE
365  l = k
366  IF( l.GT.ilo ) THEN
367 *
368 * H(L,L-1) is negligible
369 *
370  h( l, l-1 ) = zero
371  END IF
372 *
373 * Exit from loop if a submatrix of order 1 has split off.
374 *
375  IF( l.GE.i )
376  $ GO TO 140
377 *
378 * Now the active submatrix is in rows and columns L to I. If
379 * eigenvalues only are being computed, only the active submatrix
380 * need be transformed.
381 *
382  IF( .NOT.wantt ) THEN
383  i1 = l
384  i2 = i
385  END IF
386 *
387  IF( its.EQ.10 ) THEN
388 *
389 * Exceptional shift.
390 *
391  s = dat1*abs( real( h( l+1, l ) ) )
392  t = s + h( l, l )
393  ELSE IF( its.EQ.20 ) THEN
394 *
395 * Exceptional shift.
396 *
397  s = dat1*abs( real( h( i, i-1 ) ) )
398  t = s + h( i, i )
399  ELSE
400 *
401 * Wilkinson's shift.
402 *
403  t = h( i, i )
404  u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
405  s = cabs1( u )
406  IF( s.NE.rzero ) THEN
407  x = half*( h( i-1, i-1 )-t )
408  sx = cabs1( x )
409  s = max( s, cabs1( x ) )
410  y = s*sqrt( ( x / s )**2+( u / s )**2 )
411  IF( sx.GT.rzero ) THEN
412  IF( real( x / sx )*real( y )+aimag( x / sx )*
413  $ aimag( y ).LT.rzero )y = -y
414  END IF
415  t = t - u*cladiv( u, ( x+y ) )
416  END IF
417  END IF
418 *
419 * Look for two consecutive small subdiagonal elements.
420 *
421  DO 60 m = i - 1, l + 1, -1
422 *
423 * Determine the effect of starting the single-shift QR
424 * iteration at row M, and see if this would make H(M,M-1)
425 * negligible.
426 *
427  h11 = h( m, m )
428  h22 = h( m+1, m+1 )
429  h11s = h11 - t
430  h21 = real( h( m+1, m ) )
431  s = cabs1( h11s ) + abs( h21 )
432  h11s = h11s / s
433  h21 = h21 / s
434  v( 1 ) = h11s
435  v( 2 ) = h21
436  h10 = real( h( m, m-1 ) )
437  IF( abs( h10 )*abs( h21 ).LE.ulp*
438  $ ( cabs1( h11s )*( cabs1( h11 )+cabs1( h22 ) ) ) )
439  $ GO TO 70
440  60 CONTINUE
441  h11 = h( l, l )
442  h22 = h( l+1, l+1 )
443  h11s = h11 - t
444  h21 = real( h( l+1, l ) )
445  s = cabs1( h11s ) + abs( h21 )
446  h11s = h11s / s
447  h21 = h21 / s
448  v( 1 ) = h11s
449  v( 2 ) = h21
450  70 CONTINUE
451 *
452 * Single-shift QR step
453 *
454  DO 120 k = m, i - 1
455 *
456 * The first iteration of this loop determines a reflection G
457 * from the vector V and applies it from left and right to H,
458 * thus creating a nonzero bulge below the subdiagonal.
459 *
460 * Each subsequent iteration determines a reflection G to
461 * restore the Hessenberg form in the (K-1)th column, and thus
462 * chases the bulge one step toward the bottom of the active
463 * submatrix.
464 *
465 * V(2) is always real before the call to CLARFG, and hence
466 * after the call T2 ( = T1*V(2) ) is also real.
467 *
468  IF( k.GT.m )
469  $ CALL ccopy( 2, h( k, k-1 ), 1, v, 1 )
470  CALL clarfg( 2, v( 1 ), v( 2 ), 1, t1 )
471  IF( k.GT.m ) THEN
472  h( k, k-1 ) = v( 1 )
473  h( k+1, k-1 ) = zero
474  END IF
475  v2 = v( 2 )
476  t2 = real( t1*v2 )
477 *
478 * Apply G from the left to transform the rows of the matrix
479 * in columns K to I2.
480 *
481  DO 80 j = k, i2
482  sum = conjg( t1 )*h( k, j ) + t2*h( k+1, j )
483  h( k, j ) = h( k, j ) - sum
484  h( k+1, j ) = h( k+1, j ) - sum*v2
485  80 CONTINUE
486 *
487 * Apply G from the right to transform the columns of the
488 * matrix in rows I1 to min(K+2,I).
489 *
490  DO 90 j = i1, min( k+2, i )
491  sum = t1*h( j, k ) + t2*h( j, k+1 )
492  h( j, k ) = h( j, k ) - sum
493  h( j, k+1 ) = h( j, k+1 ) - sum*conjg( v2 )
494  90 CONTINUE
495 *
496  IF( wantz ) THEN
497 *
498 * Accumulate transformations in the matrix Z
499 *
500  DO 100 j = iloz, ihiz
501  sum = t1*z( j, k ) + t2*z( j, k+1 )
502  z( j, k ) = z( j, k ) - sum
503  z( j, k+1 ) = z( j, k+1 ) - sum*conjg( v2 )
504  100 CONTINUE
505  END IF
506 *
507  IF( k.EQ.m .AND. m.GT.l ) THEN
508 *
509 * If the QR step was started at row M > L because two
510 * consecutive small subdiagonals were found, then extra
511 * scaling must be performed to ensure that H(M,M-1) remains
512 * real.
513 *
514  temp = one - t1
515  temp = temp / abs( temp )
516  h( m+1, m ) = h( m+1, m )*conjg( temp )
517  IF( m+2.LE.i )
518  $ h( m+2, m+1 ) = h( m+2, m+1 )*temp
519  DO 110 j = m, i
520  IF( j.NE.m+1 ) THEN
521  IF( i2.GT.j )
522  $ CALL cscal( i2-j, temp, h( j, j+1 ), ldh )
523  CALL cscal( j-i1, conjg( temp ), h( i1, j ), 1 )
524  IF( wantz ) THEN
525  CALL cscal( nz, conjg( temp ), z( iloz, j ), 1 )
526  END IF
527  END IF
528  110 CONTINUE
529  END IF
530  120 CONTINUE
531 *
532 * Ensure that H(I,I-1) is real.
533 *
534  temp = h( i, i-1 )
535  IF( aimag( temp ).NE.rzero ) THEN
536  rtemp = abs( temp )
537  h( i, i-1 ) = rtemp
538  temp = temp / rtemp
539  IF( i2.GT.i )
540  $ CALL cscal( i2-i, conjg( temp ), h( i, i+1 ), ldh )
541  CALL cscal( i-i1, temp, h( i1, i ), 1 )
542  IF( wantz ) THEN
543  CALL cscal( nz, temp, z( iloz, i ), 1 )
544  END IF
545  END IF
546 *
547  130 CONTINUE
548 *
549 * Failure to converge in remaining number of iterations
550 *
551  info = i
552  RETURN
553 *
554  140 CONTINUE
555 *
556 * H(I,I-1) is negligible: one eigenvalue has converged.
557 *
558  w( i ) = h( i, i )
559 *
560 * return to start of the main loop with new value of I.
561 *
562  i = l - 1
563  GO TO 30
564 *
565  150 CONTINUE
566  RETURN
567 *
568 * End of CLAHQR
569 *
570  END
clahqr
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:197
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
clarfg
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:80
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83