LAPACK  3.10.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.

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
           > 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.)
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 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 193 of file clahqr.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 H( LDH, * ), W( * ), Z( LDZ, * )
207 * ..
208 *
209 * =========================================================
210 *
211 * .. Parameters ..
212  COMPLEX ZERO, ONE
213  parameter( zero = ( 0.0e0, 0.0e0 ),
214  $ one = ( 1.0e0, 0.0e0 ) )
215  REAL RZERO, RONE, HALF
216  parameter( rzero = 0.0e0, rone = 1.0e0, half = 0.5e0 )
217  REAL DAT1
218  parameter( dat1 = 3.0e0 / 4.0e0 )
219  INTEGER KEXSH
220  parameter( kexsh = 10 )
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, KDEFL
229 * ..
230 * .. Local Arrays ..
231  COMPLEX V( 2 )
232 * ..
233 * .. External Functions ..
234  COMPLEX CLADIV
235  REAL SLAMCH
236  EXTERNAL cladiv, 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 * 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( real( h( k-1, k-2 ) ) )
349  IF( k+1.LE.ihi )
350  $ tst = tst + abs( real( 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( real( 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( real( 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( real( 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( real( x / sx )*real( y )+aimag( x / sx )*
418  $ aimag( y ).LT.rzero )y = -y
419  END IF
420  t = t - u*cladiv( 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 = real( 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 = real( 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 = real( 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 CLARFG, and hence
471 * after the call T2 ( = T1*V(2) ) is also real.
472 *
473  IF( k.GT.m )
474  $ CALL ccopy( 2, h( k, k-1 ), 1, v, 1 )
475  CALL clarfg( 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 = real( 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 = conjg( 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*conjg( 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*conjg( 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 )*conjg( 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 cscal( i2-j, temp, h( j, j+1 ), ldh )
528  CALL cscal( j-i1, conjg( temp ), h( i1, j ), 1 )
529  IF( wantz ) THEN
530  CALL cscal( nz, conjg( temp ), z( iloz, j ), 1 )
531  END IF
532  END IF
533  110 CONTINUE
534  END IF
535  120 CONTINUE
536 *
537 * Ensure that H(I,I-1) is real.
538 *
539  temp = h( i, i-1 )
540  IF( aimag( temp ).NE.rzero ) THEN
541  rtemp = abs( temp )
542  h( i, i-1 ) = rtemp
543  temp = temp / rtemp
544  IF( i2.GT.i )
545  $ CALL cscal( i2-i, conjg( temp ), h( i, i+1 ), ldh )
546  CALL cscal( i-i1, temp, h( i1, i ), 1 )
547  IF( wantz ) THEN
548  CALL cscal( nz, temp, z( iloz, i ), 1 )
549  END IF
550  END IF
551 *
552  130 CONTINUE
553 *
554 * Failure to converge in remaining number of iterations
555 *
556  info = i
557  RETURN
558 *
559  140 CONTINUE
560 *
561 * H(I,I-1) is negligible: one eigenvalue has converged.
562 *
563  w( i ) = h( i, i )
564 * reset deflation counter
565  kdefl = 0
566 *
567 * return to start of the main loop with new value of I.
568 *
569  i = l - 1
570  GO TO 30
571 *
572  150 CONTINUE
573  RETURN
574 *
575 * End of CLAHQR
576 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
complex function cladiv(X, Y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: cladiv.f:64
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:106
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: