LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
zlahqr.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine zlahqr (WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
 ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. More...
 

Function/Subroutine Documentation

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
September 2012
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.4.2) --
199 * -- LAPACK is a software package provided by Univ. of Tennessee, --
200 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201 * September 2012
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  INTEGER itmax
215  parameter( itmax = 30 )
216  COMPLEX*16 zero, one
217  parameter( zero = ( 0.0d0, 0.0d0 ),
218  $ one = ( 1.0d0, 0.0d0 ) )
219  DOUBLE PRECISION rzero, rone, half
220  parameter( rzero = 0.0d0, rone = 1.0d0, half = 0.5d0 )
221  DOUBLE PRECISION dat1
222  parameter( dat1 = 3.0d0 / 4.0d0 )
223 * ..
224 * .. Local Scalars ..
225  COMPLEX*16 cdum, h11, h11s, h22, sc, sum, t, t1, temp, u,
226  $ v2, x, y
227  DOUBLE PRECISION aa, ab, ba, bb, h10, h21, rtemp, s, safmax,
228  $ safmin, smlnum, sx, t2, tst, ulp
229  INTEGER i, i1, i2, its, j, jhi, jlo, k, l, m, nh, nz
230 * ..
231 * .. Local Arrays ..
232  COMPLEX*16 v( 2 )
233 * ..
234 * .. External Functions ..
235  COMPLEX*16 zladiv
236  DOUBLE PRECISION dlamch
237  EXTERNAL zladiv, dlamch
238 * ..
239 * .. External Subroutines ..
240  EXTERNAL dlabad, zcopy, zlarfg, zscal
241 * ..
242 * .. Statement Functions ..
243  DOUBLE PRECISION cabs1
244 * ..
245 * .. Intrinsic Functions ..
246  INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
247 * ..
248 * .. Statement Function definitions ..
249  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
250 * ..
251 * .. Executable Statements ..
252 *
253  info = 0
254 *
255 * Quick return if possible
256 *
257  IF( n.EQ.0 )
258  $ RETURN
259  IF( ilo.EQ.ihi ) THEN
260  w( ilo ) = h( ilo, ilo )
261  RETURN
262  END IF
263 *
264 * ==== clear out the trash ====
265  DO 10 j = ilo, ihi - 3
266  h( j+2, j ) = zero
267  h( j+3, j ) = zero
268  10 CONTINUE
269  IF( ilo.LE.ihi-2 )
270  $ h( ihi, ihi-2 ) = zero
271 * ==== ensure that subdiagonal entries are real ====
272  IF( wantt ) THEN
273  jlo = 1
274  jhi = n
275  ELSE
276  jlo = ilo
277  jhi = ihi
278  END IF
279  DO 20 i = ilo + 1, ihi
280  IF( dimag( h( i, i-1 ) ).NE.rzero ) THEN
281 * ==== The following redundant normalization
282 * . avoids problems with both gradual and
283 * . sudden underflow in ABS(H(I,I-1)) ====
284  sc = h( i, i-1 ) / cabs1( h( i, i-1 ) )
285  sc = dconjg( sc ) / abs( sc )
286  h( i, i-1 ) = abs( h( i, i-1 ) )
287  CALL zscal( jhi-i+1, sc, h( i, i ), ldh )
288  CALL zscal( min( jhi, i+1 )-jlo+1, dconjg( sc ),
289  $ h( jlo, i ), 1 )
290  IF( wantz )
291  $ CALL zscal( ihiz-iloz+1, dconjg( sc ), z( iloz, i ), 1 )
292  END IF
293  20 CONTINUE
294 *
295  nh = ihi - ilo + 1
296  nz = ihiz - iloz + 1
297 *
298 * Set machine-dependent constants for the stopping criterion.
299 *
300  safmin = dlamch( 'SAFE MINIMUM' )
301  safmax = rone / safmin
302  CALL dlabad( safmin, safmax )
303  ulp = dlamch( 'PRECISION' )
304  smlnum = safmin*( dble( nh ) / ulp )
305 *
306 * I1 and I2 are the indices of the first row and last column of H
307 * to which transformations must be applied. If eigenvalues only are
308 * being computed, I1 and I2 are set inside the main loop.
309 *
310  IF( wantt ) THEN
311  i1 = 1
312  i2 = n
313  END IF
314 *
315 * The main loop begins here. I is the loop index and decreases from
316 * IHI to ILO in steps of 1. Each iteration of the loop works
317 * with the active submatrix in rows and columns L to I.
318 * Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
319 * H(L,L-1) is negligible so that the matrix splits.
320 *
321  i = ihi
322  30 CONTINUE
323  IF( i.LT.ilo )
324  $ GO TO 150
325 *
326 * Perform QR iterations on rows and columns ILO to I until a
327 * submatrix of order 1 splits off at the bottom because a
328 * subdiagonal element has become negligible.
329 *
330  l = ilo
331  DO 130 its = 0, itmax
332 *
333 * Look for a single small subdiagonal element.
334 *
335  DO 40 k = i, l + 1, -1
336  IF( cabs1( h( k, k-1 ) ).LE.smlnum )
337  $ GO TO 50
338  tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
339  IF( tst.EQ.zero ) THEN
340  IF( k-2.GE.ilo )
341  $ tst = tst + abs( dble( h( k-1, k-2 ) ) )
342  IF( k+1.LE.ihi )
343  $ tst = tst + abs( dble( h( k+1, k ) ) )
344  END IF
345 * ==== The following is a conservative small subdiagonal
346 * . deflation criterion due to Ahues & Tisseur (LAWN 122,
347 * . 1997). It has better mathematical foundation and
348 * . improves accuracy in some examples. ====
349  IF( abs( dble( h( k, k-1 ) ) ).LE.ulp*tst ) THEN
350  ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
351  ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
352  aa = max( cabs1( h( k, k ) ),
353  $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
354  bb = min( cabs1( h( k, k ) ),
355  $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
356  s = aa + ab
357  IF( ba*( ab / s ).LE.max( smlnum,
358  $ ulp*( bb*( aa / s ) ) ) )GO TO 50
359  END IF
360  40 CONTINUE
361  50 CONTINUE
362  l = k
363  IF( l.GT.ilo ) THEN
364 *
365 * H(L,L-1) is negligible
366 *
367  h( l, l-1 ) = zero
368  END IF
369 *
370 * Exit from loop if a submatrix of order 1 has split off.
371 *
372  IF( l.GE.i )
373  $ GO TO 140
374 *
375 * Now the active submatrix is in rows and columns L to I. If
376 * eigenvalues only are being computed, only the active submatrix
377 * need be transformed.
378 *
379  IF( .NOT.wantt ) THEN
380  i1 = l
381  i2 = i
382  END IF
383 *
384  IF( its.EQ.10 ) THEN
385 *
386 * Exceptional shift.
387 *
388  s = dat1*abs( dble( h( l+1, l ) ) )
389  t = s + h( l, l )
390  ELSE IF( its.EQ.20 ) THEN
391 *
392 * Exceptional shift.
393 *
394  s = dat1*abs( dble( h( i, i-1 ) ) )
395  t = s + h( i, i )
396  ELSE
397 *
398 * Wilkinson's shift.
399 *
400  t = h( i, i )
401  u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
402  s = cabs1( u )
403  IF( s.NE.rzero ) THEN
404  x = half*( h( i-1, i-1 )-t )
405  sx = cabs1( x )
406  s = max( s, cabs1( x ) )
407  y = s*sqrt( ( x / s )**2+( u / s )**2 )
408  IF( sx.GT.rzero ) THEN
409  IF( dble( x / sx )*dble( y )+dimag( x / sx )*
410  $ dimag( y ).LT.rzero )y = -y
411  END IF
412  t = t - u*zladiv( u, ( x+y ) )
413  END IF
414  END IF
415 *
416 * Look for two consecutive small subdiagonal elements.
417 *
418  DO 60 m = i - 1, l + 1, -1
419 *
420 * Determine the effect of starting the single-shift QR
421 * iteration at row M, and see if this would make H(M,M-1)
422 * negligible.
423 *
424  h11 = h( m, m )
425  h22 = h( m+1, m+1 )
426  h11s = h11 - t
427  h21 = dble( h( m+1, m ) )
428  s = cabs1( h11s ) + abs( h21 )
429  h11s = h11s / s
430  h21 = h21 / s
431  v( 1 ) = h11s
432  v( 2 ) = h21
433  h10 = dble( h( m, m-1 ) )
434  IF( abs( h10 )*abs( h21 ).LE.ulp*
435  $ ( cabs1( h11s )*( cabs1( h11 )+cabs1( h22 ) ) ) )
436  $ GO TO 70
437  60 CONTINUE
438  h11 = h( l, l )
439  h22 = h( l+1, l+1 )
440  h11s = h11 - t
441  h21 = dble( h( l+1, l ) )
442  s = cabs1( h11s ) + abs( h21 )
443  h11s = h11s / s
444  h21 = h21 / s
445  v( 1 ) = h11s
446  v( 2 ) = h21
447  70 CONTINUE
448 *
449 * Single-shift QR step
450 *
451  DO 120 k = m, i - 1
452 *
453 * The first iteration of this loop determines a reflection G
454 * from the vector V and applies it from left and right to H,
455 * thus creating a nonzero bulge below the subdiagonal.
456 *
457 * Each subsequent iteration determines a reflection G to
458 * restore the Hessenberg form in the (K-1)th column, and thus
459 * chases the bulge one step toward the bottom of the active
460 * submatrix.
461 *
462 * V(2) is always real before the call to ZLARFG, and hence
463 * after the call T2 ( = T1*V(2) ) is also real.
464 *
465  IF( k.GT.m )
466  $ CALL zcopy( 2, h( k, k-1 ), 1, v, 1 )
467  CALL zlarfg( 2, v( 1 ), v( 2 ), 1, t1 )
468  IF( k.GT.m ) THEN
469  h( k, k-1 ) = v( 1 )
470  h( k+1, k-1 ) = zero
471  END IF
472  v2 = v( 2 )
473  t2 = dble( t1*v2 )
474 *
475 * Apply G from the left to transform the rows of the matrix
476 * in columns K to I2.
477 *
478  DO 80 j = k, i2
479  sum = dconjg( t1 )*h( k, j ) + t2*h( k+1, j )
480  h( k, j ) = h( k, j ) - sum
481  h( k+1, j ) = h( k+1, j ) - sum*v2
482  80 CONTINUE
483 *
484 * Apply G from the right to transform the columns of the
485 * matrix in rows I1 to min(K+2,I).
486 *
487  DO 90 j = i1, min( k+2, i )
488  sum = t1*h( j, k ) + t2*h( j, k+1 )
489  h( j, k ) = h( j, k ) - sum
490  h( j, k+1 ) = h( j, k+1 ) - sum*dconjg( v2 )
491  90 CONTINUE
492 *
493  IF( wantz ) THEN
494 *
495 * Accumulate transformations in the matrix Z
496 *
497  DO 100 j = iloz, ihiz
498  sum = t1*z( j, k ) + t2*z( j, k+1 )
499  z( j, k ) = z( j, k ) - sum
500  z( j, k+1 ) = z( j, k+1 ) - sum*dconjg( v2 )
501  100 CONTINUE
502  END IF
503 *
504  IF( k.EQ.m .AND. m.GT.l ) THEN
505 *
506 * If the QR step was started at row M > L because two
507 * consecutive small subdiagonals were found, then extra
508 * scaling must be performed to ensure that H(M,M-1) remains
509 * real.
510 *
511  temp = one - t1
512  temp = temp / abs( temp )
513  h( m+1, m ) = h( m+1, m )*dconjg( temp )
514  IF( m+2.LE.i )
515  $ h( m+2, m+1 ) = h( m+2, m+1 )*temp
516  DO 110 j = m, i
517  IF( j.NE.m+1 ) THEN
518  IF( i2.GT.j )
519  $ CALL zscal( i2-j, temp, h( j, j+1 ), ldh )
520  CALL zscal( j-i1, dconjg( temp ), h( i1, j ), 1 )
521  IF( wantz ) THEN
522  CALL zscal( nz, dconjg( temp ), z( iloz, j ),
523  $ 1 )
524  END IF
525  END IF
526  110 CONTINUE
527  END IF
528  120 CONTINUE
529 *
530 * Ensure that H(I,I-1) is real.
531 *
532  temp = h( i, i-1 )
533  IF( dimag( temp ).NE.rzero ) THEN
534  rtemp = abs( temp )
535  h( i, i-1 ) = rtemp
536  temp = temp / rtemp
537  IF( i2.GT.i )
538  $ CALL zscal( i2-i, dconjg( temp ), h( i, i+1 ), ldh )
539  CALL zscal( i-i1, temp, h( i1, i ), 1 )
540  IF( wantz ) THEN
541  CALL zscal( nz, temp, z( iloz, i ), 1 )
542  END IF
543  END IF
544 *
545  130 CONTINUE
546 *
547 * Failure to converge in remaining number of iterations
548 *
549  info = i
550  RETURN
551 *
552  140 CONTINUE
553 *
554 * H(I,I-1) is negligible: one eigenvalue has converged.
555 *
556  w( i ) = h( i, i )
557 *
558 * return to start of the main loop with new value of I.
559 *
560  i = l - 1
561  GO TO 30
562 *
563  150 CONTINUE
564  RETURN
565 *
566 * End of ZLAHQR
567 *
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
complex *16 function zladiv(X, Y)
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: zladiv.f:66
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108
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: