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

Go to the source code of this file.

Functions/Subroutines

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, using the double-shift/single-shift QR algorithm. More...
 

Function/Subroutine Documentation

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.

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

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
          .GT. 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 .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 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.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 h( ldh, * ), w( * ), z( ldz, * )
209 * ..
210 *
211 * =========================================================
212 *
213 * .. Parameters ..
214  INTEGER itmax
215  parameter( itmax = 30 )
216  COMPLEX zero, one
217  parameter( zero = ( 0.0e0, 0.0e0 ),
218  $ one = ( 1.0e0, 0.0e0 ) )
219  REAL rzero, rone, half
220  parameter( rzero = 0.0e0, rone = 1.0e0, half = 0.5e0 )
221  REAL dat1
222  parameter( dat1 = 3.0e0 / 4.0e0 )
223 * ..
224 * .. Local Scalars ..
225  COMPLEX cdum, h11, h11s, h22, sc, sum, t, t1, temp, u,
226  $ v2, x, y
227  REAL 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 v( 2 )
233 * ..
234 * .. External Functions ..
235  COMPLEX cladiv
236  REAL slamch
237  EXTERNAL cladiv, slamch
238 * ..
239 * .. External Subroutines ..
240  EXTERNAL ccopy, clarfg, cscal, slabad
241 * ..
242 * .. Statement Functions ..
243  REAL cabs1
244 * ..
245 * .. Intrinsic Functions ..
246  INTRINSIC abs, aimag, conjg, max, min, REAL, sqrt
247 * ..
248 * .. Statement Function definitions ..
249  cabs1( cdum ) = abs( REAL( CDUM ) ) + abs( aimag( 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( aimag( 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 = conjg( sc ) / abs( sc )
286  h( i, i-1 ) = abs( h( i, i-1 ) )
287  CALL cscal( jhi-i+1, sc, h( i, i ), ldh )
288  CALL cscal( min( jhi, i+1 )-jlo+1, conjg( sc ), h( jlo, i ),
289  $ 1 )
290  IF( wantz )
291  $ CALL cscal( ihiz-iloz+1, conjg( 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 = slamch( 'SAFE MINIMUM' )
301  safmax = rone / safmin
302  CALL slabad( safmin, safmax )
303  ulp = slamch( 'PRECISION' )
304  smlnum = safmin*( REAL( 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( REAL( H( K-1, K-2 ) ) )
342  IF( k+1.LE.ihi )
343  $ tst = tst + abs( REAL( 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( REAL( 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( REAL( 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( REAL( 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( REAL( x / sx )*REAL( y )+aimag( x / sx )*
410  $ aimag( y ).LT.rzero )y = -y
411  END IF
412  t = t - u*cladiv( 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 = REAL( 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 = REAL( 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 = REAL( 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 CLARFG, and hence
463 * after the call T2 ( = T1*V(2) ) is also real.
464 *
465  IF( k.GT.m )
466  $ CALL ccopy( 2, h( k, k-1 ), 1, v, 1 )
467  CALL clarfg( 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 = REAL( 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 = conjg( 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*conjg( 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*conjg( 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 )*conjg( 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 cscal( i2-j, temp, h( j, j+1 ), ldh )
520  CALL cscal( j-i1, conjg( temp ), h( i1, j ), 1 )
521  IF( wantz ) THEN
522  CALL cscal( nz, conjg( temp ), z( iloz, j ), 1 )
523  END IF
524  END IF
525  110 CONTINUE
526  END IF
527  120 CONTINUE
528 *
529 * Ensure that H(I,I-1) is real.
530 *
531  temp = h( i, i-1 )
532  IF( aimag( temp ).NE.rzero ) THEN
533  rtemp = abs( temp )
534  h( i, i-1 ) = rtemp
535  temp = temp / rtemp
536  IF( i2.GT.i )
537  $ CALL cscal( i2-i, conjg( temp ), h( i, i+1 ), ldh )
538  CALL cscal( i-i1, temp, h( i1, i ), 1 )
539  IF( wantz ) THEN
540  CALL cscal( nz, temp, z( iloz, i ), 1 )
541  END IF
542  END IF
543 *
544  130 CONTINUE
545 *
546 * Failure to converge in remaining number of iterations
547 *
548  info = i
549  RETURN
550 *
551  140 CONTINUE
552 *
553 * H(I,I-1) is negligible: one eigenvalue has converged.
554 *
555  w( i ) = h( i, i )
556 *
557 * return to start of the main loop with new value of I.
558 *
559  i = l - 1
560  GO TO 30
561 *
562  150 CONTINUE
563  RETURN
564 *
565 * End of CLAHQR
566 *
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
complex function cladiv(X, Y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: cladiv.f:66
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52

Here is the call graph for this function:

Here is the caller graph for this function: