LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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
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 ulp = slamch( 'PRECISION' )
302 smlnum = safmin*( real( nh ) / ulp )
303*
304* I1 and I2 are the indices of the first row and last column of H
305* to which transformations must be applied. If eigenvalues only are
306* being computed, I1 and I2 are set inside the main loop.
307*
308 IF( wantt ) THEN
309 i1 = 1
310 i2 = n
311 END IF
312*
313* ITMAX is the total number of QR iterations allowed.
314*
315 itmax = 30 * max( 10, nh )
316*
317* KDEFL counts the number of iterations since a deflation
318*
319 kdefl = 0
320*
321* The main loop begins here. I is the loop index and decreases from
322* IHI to ILO in steps of 1. Each iteration of the loop works
323* with the active submatrix in rows and columns L to I.
324* Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
325* H(L,L-1) is negligible so that the matrix splits.
326*
327 i = ihi
328 30 CONTINUE
329 IF( i.LT.ilo )
330 $ GO TO 150
331*
332* Perform QR iterations on rows and columns ILO to I until a
333* submatrix of order 1 splits off at the bottom because a
334* subdiagonal element has become negligible.
335*
336 l = ilo
337 DO 130 its = 0, itmax
338*
339* Look for a single small subdiagonal element.
340*
341 DO 40 k = i, l + 1, -1
342 IF( cabs1( h( k, k-1 ) ).LE.smlnum )
343 $ GO TO 50
344 tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
345 IF( tst.EQ.zero ) THEN
346 IF( k-2.GE.ilo )
347 $ tst = tst + abs( real( h( k-1, k-2 ) ) )
348 IF( k+1.LE.ihi )
349 $ tst = tst + abs( real( h( k+1, k ) ) )
350 END IF
351* ==== The following is a conservative small subdiagonal
352* . deflation criterion due to Ahues & Tisseur (LAWN 122,
353* . 1997). It has better mathematical foundation and
354* . improves accuracy in some examples. ====
355 IF( abs( real( h( k, k-1 ) ) ).LE.ulp*tst ) THEN
356 ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
357 ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
358 aa = max( cabs1( h( k, k ) ),
359 $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
360 bb = min( cabs1( h( k, k ) ),
361 $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
362 s = aa + ab
363 IF( ba*( ab / s ).LE.max( smlnum,
364 $ ulp*( bb*( aa / s ) ) ) )GO TO 50
365 END IF
366 40 CONTINUE
367 50 CONTINUE
368 l = k
369 IF( l.GT.ilo ) THEN
370*
371* H(L,L-1) is negligible
372*
373 h( l, l-1 ) = zero
374 END IF
375*
376* Exit from loop if a submatrix of order 1 has split off.
377*
378 IF( l.GE.i )
379 $ GO TO 140
380 kdefl = kdefl + 1
381*
382* Now the active submatrix is in rows and columns L to I. If
383* eigenvalues only are being computed, only the active submatrix
384* need be transformed.
385*
386 IF( .NOT.wantt ) THEN
387 i1 = l
388 i2 = i
389 END IF
390*
391 IF( mod(kdefl,2*kexsh).EQ.0 ) THEN
392*
393* Exceptional shift.
394*
395 s = dat1*abs( real( h( i, i-1 ) ) )
396 t = s + h( i, i )
397 ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
398*
399* Exceptional shift.
400*
401 s = dat1*abs( real( h( l+1, l ) ) )
402 t = s + h( l, l )
403 ELSE
404*
405* Wilkinson's shift.
406*
407 t = h( i, i )
408 u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
409 s = cabs1( u )
410 IF( s.NE.rzero ) THEN
411 x = half*( h( i-1, i-1 )-t )
412 sx = cabs1( x )
413 s = max( s, cabs1( x ) )
414 y = s*sqrt( ( x / s )**2+( u / s )**2 )
415 IF( sx.GT.rzero ) THEN
416 IF( real( x / sx )*real( y )+aimag( x / sx )*
417 $ aimag( y ).LT.rzero )y = -y
418 END IF
419 t = t - u*cladiv( u, ( x+y ) )
420 END IF
421 END IF
422*
423* Look for two consecutive small subdiagonal elements.
424*
425 DO 60 m = i - 1, l + 1, -1
426*
427* Determine the effect of starting the single-shift QR
428* iteration at row M, and see if this would make H(M,M-1)
429* negligible.
430*
431 h11 = h( m, m )
432 h22 = h( m+1, m+1 )
433 h11s = h11 - t
434 h21 = real( h( m+1, m ) )
435 s = cabs1( h11s ) + abs( h21 )
436 h11s = h11s / s
437 h21 = h21 / s
438 v( 1 ) = h11s
439 v( 2 ) = h21
440 h10 = real( h( m, m-1 ) )
441 IF( abs( h10 )*abs( h21 ).LE.ulp*
442 $ ( cabs1( h11s )*( cabs1( h11 )+cabs1( h22 ) ) ) )
443 $ GO TO 70
444 60 CONTINUE
445 h11 = h( l, l )
446 h22 = h( l+1, l+1 )
447 h11s = h11 - t
448 h21 = real( h( l+1, l ) )
449 s = cabs1( h11s ) + abs( h21 )
450 h11s = h11s / s
451 h21 = h21 / s
452 v( 1 ) = h11s
453 v( 2 ) = h21
454 70 CONTINUE
455*
456* Single-shift QR step
457*
458 DO 120 k = m, i - 1
459*
460* The first iteration of this loop determines a reflection G
461* from the vector V and applies it from left and right to H,
462* thus creating a nonzero bulge below the subdiagonal.
463*
464* Each subsequent iteration determines a reflection G to
465* restore the Hessenberg form in the (K-1)th column, and thus
466* chases the bulge one step toward the bottom of the active
467* submatrix.
468*
469* V(2) is always real before the call to CLARFG, and hence
470* after the call T2 ( = T1*V(2) ) is also real.
471*
472 IF( k.GT.m )
473 $ CALL ccopy( 2, h( k, k-1 ), 1, v, 1 )
474 CALL clarfg( 2, v( 1 ), v( 2 ), 1, t1 )
475 IF( k.GT.m ) THEN
476 h( k, k-1 ) = v( 1 )
477 h( k+1, k-1 ) = zero
478 END IF
479 v2 = v( 2 )
480 t2 = real( t1*v2 )
481*
482* Apply G from the left to transform the rows of the matrix
483* in columns K to I2.
484*
485 DO 80 j = k, i2
486 sum = conjg( t1 )*h( k, j ) + t2*h( k+1, j )
487 h( k, j ) = h( k, j ) - sum
488 h( k+1, j ) = h( k+1, j ) - sum*v2
489 80 CONTINUE
490*
491* Apply G from the right to transform the columns of the
492* matrix in rows I1 to min(K+2,I).
493*
494 DO 90 j = i1, min( k+2, i )
495 sum = t1*h( j, k ) + t2*h( j, k+1 )
496 h( j, k ) = h( j, k ) - sum
497 h( j, k+1 ) = h( j, k+1 ) - sum*conjg( v2 )
498 90 CONTINUE
499*
500 IF( wantz ) THEN
501*
502* Accumulate transformations in the matrix Z
503*
504 DO 100 j = iloz, ihiz
505 sum = t1*z( j, k ) + t2*z( j, k+1 )
506 z( j, k ) = z( j, k ) - sum
507 z( j, k+1 ) = z( j, k+1 ) - sum*conjg( v2 )
508 100 CONTINUE
509 END IF
510*
511 IF( k.EQ.m .AND. m.GT.l ) THEN
512*
513* If the QR step was started at row M > L because two
514* consecutive small subdiagonals were found, then extra
515* scaling must be performed to ensure that H(M,M-1) remains
516* real.
517*
518 temp = one - t1
519 temp = temp / abs( temp )
520 h( m+1, m ) = h( m+1, m )*conjg( temp )
521 IF( m+2.LE.i )
522 $ h( m+2, m+1 ) = h( m+2, m+1 )*temp
523 DO 110 j = m, i
524 IF( j.NE.m+1 ) THEN
525 IF( i2.GT.j )
526 $ CALL cscal( i2-j, temp, h( j, j+1 ), ldh )
527 CALL cscal( j-i1, conjg( temp ), h( i1, j ), 1 )
528 IF( wantz ) THEN
529 CALL cscal( nz, conjg( temp ), z( iloz, j ), 1 )
530 END IF
531 END IF
532 110 CONTINUE
533 END IF
534 120 CONTINUE
535*
536* Ensure that H(I,I-1) is real.
537*
538 temp = h( i, i-1 )
539 IF( aimag( temp ).NE.rzero ) THEN
540 rtemp = abs( temp )
541 h( i, i-1 ) = rtemp
542 temp = temp / rtemp
543 IF( i2.GT.i )
544 $ CALL cscal( i2-i, conjg( temp ), h( i, i+1 ), ldh )
545 CALL cscal( i-i1, temp, h( i1, i ), 1 )
546 IF( wantz ) THEN
547 CALL cscal( nz, temp, z( iloz, i ), 1 )
548 END IF
549 END IF
550*
551 130 CONTINUE
552*
553* Failure to converge in remaining number of iterations
554*
555 info = i
556 RETURN
557*
558 140 CONTINUE
559*
560* H(I,I-1) is negligible: one eigenvalue has converged.
561*
562 w( i ) = h( i, i )
563* reset deflation counter
564 kdefl = 0
565*
566* return to start of the main loop with new value of I.
567*
568 i = l - 1
569 GO TO 30
570*
571 150 CONTINUE
572 RETURN
573*
574* End of CLAHQR
575*
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
complex function cladiv(x, y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition cladiv.f:64
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
subroutine clarfg(n, alpha, x, incx, tau)
CLARFG generates an elementary reflector (Householder matrix).
Definition clarfg.f:106
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
Here is the call graph for this function:
Here is the caller graph for this function: