 LAPACK  3.9.0 LAPACK: Linear Algebra PACKage

## ◆ clahqr()

 subroutine clahqr ( logical WANTT, logical WANTZ, integer N, integer ILO, integer IHI, complex, dimension( ldh, * ) H, integer LDH, complex, dimension( * ) W, integer ILOZ, integer IHIZ, complex, dimension( ldz, * ) Z, integer LDZ, integer INFO )

CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.

Purpose:
```    CLAHQR is an auxiliary routine called by CHSEQR to update the
eigenvalues and Schur decomposition already computed by CHSEQR, by
dealing with the Hessenberg submatrix in rows and columns ILO to
IHI.```
Parameters
 [in] WANTT ``` WANTT is LOGICAL = .TRUE. : the full Schur form T is required; = .FALSE.: only eigenvalues are required.``` [in] WANTZ ``` WANTZ is LOGICAL = .TRUE. : the matrix of Schur vectors Z is required; = .FALSE.: Schur vectors are not required.``` [in] N ``` N is INTEGER The order of the matrix H. N >= 0.``` [in] ILO ` ILO is INTEGER` [in] IHI ``` IHI is INTEGER It is assumed that H is already upper triangular in rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). CLAHQR works primarily with the Hessenberg submatrix in rows and columns ILO to IHI, but applies transformations to all of H if WANTT is .TRUE.. 1 <= ILO <= max(1,IHI); IHI <= N.``` [in,out] H ``` H is COMPLEX array, dimension (LDH,N) On entry, the upper Hessenberg matrix H. On exit, if INFO is zero and if WANTT is .TRUE., then H is upper triangular in rows and columns ILO:IHI. If INFO is zero and if WANTT is .FALSE., then the contents of H are unspecified on exit. The output state of H in case INF is positive is below under the description of INFO.``` [in] LDH ``` LDH is INTEGER The leading dimension of the array H. LDH >= max(1,N).``` [out] W ``` W is COMPLEX array, dimension (N) The computed eigenvalues ILO to IHI are stored in the corresponding elements of W. If WANTT is .TRUE., the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with W(i) = H(i,i).``` [in] ILOZ ` ILOZ is INTEGER` [in] IHIZ ``` IHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.``` [in,out] Z ``` Z is COMPLEX array, dimension (LDZ,N) If WANTZ is .TRUE., on entry Z must contain the current matrix Z of transformations accumulated by CHSEQR, and on exit Z has been updated; transformations are applied only to the submatrix Z(ILOZ:IHIZ,ILO:IHI). If WANTZ is .FALSE., Z is not referenced.``` [in] LDZ ``` LDZ is INTEGER The leading dimension of the array Z. LDZ >= max(1,N).``` [out] INFO ``` INFO is INTEGER = 0: successful exit > 0: if INFO = i, CLAHQR failed to compute all the eigenvalues ILO to IHI in a total of 30 iterations per eigenvalue; elements i+1:ihi of W contain those eigenvalues which have been successfully computed. If INFO > 0 and WANTT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H. If INFO > 0 and WANTT is .TRUE., then on exit (*) (initial value of H)*U = U*(final value of H) where U is an orthogonal matrix. The final value of H is upper Hessenberg and triangular in rows and columns INFO+1 through IHI. If INFO > 0 and WANTZ is .TRUE., then on exit (final value of Z) = (initial value of Z)*U where U is the orthogonal matrix in (*) (regardless of the value of WANTT.)```
Date
December 2016
Contributors:
```     02-96 Based on modifications by
David Day, Sandia National Laboratory, USA

12-04 Further modifications by
Ralph Byers, University of Kansas, USA
This is a modified version of CLAHQR from LAPACK version 3.0.
It is (1) more robust against overflow and underflow and
(2) adopts the more conservative Ahues & Tisseur stopping
criterion (LAWN 122, 1997).```

Definition at line 197 of file clahqr.f.

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 ..
235  REAL 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 *
Here is the call graph for this function:
Here is the caller graph for this function:
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
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:70
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83