LAPACK 3.11.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, 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: