LAPACK  3.9.1
LAPACK: Linear Algebra PACKage

◆ zlahqr()

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
           > 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 > 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.)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 193 of file zlahqr.f.

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*16 H( LDH, * ), W( * ), Z( LDZ, * )
207 * ..
208 *
209 * =========================================================
210 *
211 * .. Parameters ..
212  COMPLEX*16 ZERO, ONE
213  parameter( zero = ( 0.0d0, 0.0d0 ),
214  $ one = ( 1.0d0, 0.0d0 ) )
215  DOUBLE PRECISION RZERO, RONE, HALF
216  parameter( rzero = 0.0d0, rone = 1.0d0, half = 0.5d0 )
217  DOUBLE PRECISION DAT1
218  parameter( dat1 = 3.0d0 / 4.0d0 )
219  INTEGER KEXSH
220  parameter( kexsh = 10 )
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, KDEFL
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 * 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( dble( h( k-1, k-2 ) ) )
349  IF( k+1.LE.ihi )
350  $ tst = tst + abs( dble( 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( dble( 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( dble( 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( dble( 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( dble( x / sx )*dble( y )+dimag( x / sx )*
418  $ dimag( y ).LT.rzero )y = -y
419  END IF
420  t = t - u*zladiv( 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 = dble( 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 = dble( 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 = dble( 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 ZLARFG, and hence
471 * after the call T2 ( = T1*V(2) ) is also real.
472 *
473  IF( k.GT.m )
474  $ CALL zcopy( 2, h( k, k-1 ), 1, v, 1 )
475  CALL zlarfg( 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 = dble( 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 = dconjg( 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*dconjg( 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*dconjg( 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 )*dconjg( 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 zscal( i2-j, temp, h( j, j+1 ), ldh )
528  CALL zscal( j-i1, dconjg( temp ), h( i1, j ), 1 )
529  IF( wantz ) THEN
530  CALL zscal( nz, dconjg( temp ), z( iloz, j ),
531  $ 1 )
532  END IF
533  END IF
534  110 CONTINUE
535  END IF
536  120 CONTINUE
537 *
538 * Ensure that H(I,I-1) is real.
539 *
540  temp = h( i, i-1 )
541  IF( dimag( temp ).NE.rzero ) THEN
542  rtemp = abs( temp )
543  h( i, i-1 ) = rtemp
544  temp = temp / rtemp
545  IF( i2.GT.i )
546  $ CALL zscal( i2-i, dconjg( temp ), h( i, i+1 ), ldh )
547  CALL zscal( i-i1, temp, h( i1, i ), 1 )
548  IF( wantz ) THEN
549  CALL zscal( nz, temp, z( iloz, i ), 1 )
550  END IF
551  END IF
552 *
553  130 CONTINUE
554 *
555 * Failure to converge in remaining number of iterations
556 *
557  info = i
558  RETURN
559 *
560  140 CONTINUE
561 *
562 * H(I,I-1) is negligible: one eigenvalue has converged.
563 *
564  w( i ) = h( i, i )
565 * reset deflation counter
566  kdefl = 0
567 *
568 * return to start of the main loop with new value of I.
569 *
570  i = l - 1
571  GO TO 30
572 *
573  150 CONTINUE
574  RETURN
575 *
576 * End of ZLAHQR
577 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
complex *16 function zladiv(X, Y)
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: zladiv.f:64
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:106
Here is the call graph for this function:
Here is the caller graph for this function: