LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zlahqr ( logical  WANTT,
logical  WANTZ,
integer  N,
integer  ILO,
integer  IHI,
complex*16, dimension( ldh, * )  H,
integer  LDH,
complex*16, dimension( * )  W,
integer  ILOZ,
integer  IHIZ,
complex*16, dimension( ldz, * )  Z,
integer  LDZ,
integer  INFO 
)

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

Download ZLAHQR + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
    ZLAHQR 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).
          ZLAHQR 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*16 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*16 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*16 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
          .GT. 0: if INFO = i, ZLAHQR 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 .GT. 0 and WANTT is .FALSE., then on exit,
                  the remaining unconverged eigenvalues are the
                  eigenvalues of the upper Hessenberg matrix
                  rows and columns ILO thorugh INFO of the final,
                  output value of H.

                  If INFO .GT. 0 and WANTT is .TRUE., then on exit
          (*)       (initial value of H)*U  = U*(final value of H)
                  where U is an orthognal matrix.    The final
                  value of H is upper Hessenberg and triangular in
                  rows and columns INFO+1 through IHI.

                  If INFO .GT. 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.)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
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 ZLAHQR 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 zlahqr.f.

197 *
198 * -- LAPACK auxiliary routine (version 3.6.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 * November 2015
202 *
203 * .. Scalar Arguments ..
204  INTEGER ihi, ihiz, ilo, iloz, info, ldh, ldz, n
205  LOGICAL wantt, wantz
206 * ..
207 * .. Array Arguments ..
208  COMPLEX*16 h( ldh, * ), w( * ), z( ldz, * )
209 * ..
210 *
211 * =========================================================
212 *
213 * .. Parameters ..
214  COMPLEX*16 zero, one
215  parameter ( zero = ( 0.0d0, 0.0d0 ),
216  $ one = ( 1.0d0, 0.0d0 ) )
217  DOUBLE PRECISION rzero, rone, half
218  parameter ( rzero = 0.0d0, rone = 1.0d0, half = 0.5d0 )
219  DOUBLE PRECISION dat1
220  parameter ( dat1 = 3.0d0 / 4.0d0 )
221 * ..
222 * .. Local Scalars ..
223  COMPLEX*16 cdum, h11, h11s, h22, sc, sum, t, t1, temp, u,
224  $ v2, x, y
225  DOUBLE PRECISION 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*16 v( 2 )
232 * ..
233 * .. External Functions ..
234  COMPLEX*16 zladiv
235  DOUBLE PRECISION dlamch
236  EXTERNAL zladiv, dlamch
237 * ..
238 * .. External Subroutines ..
239  EXTERNAL dlabad, zcopy, zlarfg, zscal
240 * ..
241 * .. Statement Functions ..
242  DOUBLE PRECISION cabs1
243 * ..
244 * .. Intrinsic Functions ..
245  INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
246 * ..
247 * .. Statement Function definitions ..
248  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( 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( dimag( 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 = dconjg( sc ) / abs( sc )
285  h( i, i-1 ) = abs( h( i, i-1 ) )
286  CALL zscal( jhi-i+1, sc, h( i, i ), ldh )
287  CALL zscal( min( jhi, i+1 )-jlo+1, dconjg( sc ),
288  $ h( jlo, i ), 1 )
289  IF( wantz )
290  $ CALL zscal( ihiz-iloz+1, dconjg( 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 = dlamch( 'SAFE MINIMUM' )
300  safmax = rone / safmin
301  CALL dlabad( safmin, safmax )
302  ulp = dlamch( 'PRECISION' )
303  smlnum = safmin*( dble( 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( dble( h( k-1, k-2 ) ) )
345  IF( k+1.LE.ihi )
346  $ tst = tst + abs( dble( 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( dble( 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( dble( 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( dble( 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( dble( x / sx )*dble( y )+dimag( x / sx )*
413  $ dimag( y ).LT.rzero )y = -y
414  END IF
415  t = t - u*zladiv( 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 = dble( 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 = dble( 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 = dble( 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 ZLARFG, and hence
466 * after the call T2 ( = T1*V(2) ) is also real.
467 *
468  IF( k.GT.m )
469  $ CALL zcopy( 2, h( k, k-1 ), 1, v, 1 )
470  CALL zlarfg( 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 = dble( 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 = dconjg( 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*dconjg( 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*dconjg( 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 )*dconjg( 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 zscal( i2-j, temp, h( j, j+1 ), ldh )
523  CALL zscal( j-i1, dconjg( temp ), h( i1, j ), 1 )
524  IF( wantz ) THEN
525  CALL zscal( nz, dconjg( temp ), z( iloz, j ),
526  $ 1 )
527  END IF
528  END IF
529  110 CONTINUE
530  END IF
531  120 CONTINUE
532 *
533 * Ensure that H(I,I-1) is real.
534 *
535  temp = h( i, i-1 )
536  IF( dimag( temp ).NE.rzero ) THEN
537  rtemp = abs( temp )
538  h( i, i-1 ) = rtemp
539  temp = temp / rtemp
540  IF( i2.GT.i )
541  $ CALL zscal( i2-i, dconjg( temp ), h( i, i+1 ), ldh )
542  CALL zscal( i-i1, temp, h( i1, i ), 1 )
543  IF( wantz ) THEN
544  CALL zscal( nz, temp, z( iloz, i ), 1 )
545  END IF
546  END IF
547 *
548  130 CONTINUE
549 *
550 * Failure to converge in remaining number of iterations
551 *
552  info = i
553  RETURN
554 *
555  140 CONTINUE
556 *
557 * H(I,I-1) is negligible: one eigenvalue has converged.
558 *
559  w( i ) = h( i, i )
560 *
561 * return to start of the main loop with new value of I.
562 *
563  i = l - 1
564  GO TO 30
565 *
566  150 CONTINUE
567  RETURN
568 *
569 * End of ZLAHQR
570 *
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
zlarfg
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
zladiv
complex *16 function zladiv(X, Y)
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: zladiv.f:66
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: