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

◆ dlahqr()

subroutine dlahqr ( logical wantt,
logical wantz,
integer n,
integer ilo,
integer ihi,
double precision, dimension( ldh, * ) h,
integer ldh,
double precision, dimension( * ) wr,
double precision, dimension( * ) wi,
integer iloz,
integer ihiz,
double precision, dimension( ldz, * ) z,
integer ldz,
integer info )

DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.

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

Purpose:
!>
!>    DLAHQR is an auxiliary routine called by DHSEQR to update the
!>    eigenvalues and Schur decomposition already computed by DHSEQR, 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 quasi-triangular in
!>          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
!>          ILO = 1). DLAHQR 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 DOUBLE PRECISION array, dimension (LDH,N)
!>          On entry, the upper Hessenberg matrix H.
!>          On exit, if INFO is zero and if WANTT is .TRUE., H is upper
!>          quasi-triangular in rows and columns ILO:IHI, with any
!>          2-by-2 diagonal blocks in standard form. If INFO is zero
!>          and WANTT is .FALSE., the contents of H are unspecified on
!>          exit.  The output state of H if INFO is nonzero is given
!>          below under the description of INFO.
!>
[in]LDH
!>          LDH is INTEGER
!>          The leading dimension of the array H. LDH >= max(1,N).
!>
[out]WR
!>          WR is DOUBLE PRECISION array, dimension (N)
!>
[out]WI
!>          WI is DOUBLE PRECISION array, dimension (N)
!>          The real and imaginary parts, respectively, of the computed
!>          eigenvalues ILO to IHI are stored in the corresponding
!>          elements of WR and WI. If two eigenvalues are computed as a
!>          complex conjugate pair, they are stored in consecutive
!>          elements of WR and WI, say the i-th and (i+1)th, with
!>          WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
!>          eigenvalues are stored in the same order as on the diagonal
!>          of the Schur form returned in H, with WR(i) = H(i,i), and, if
!>          H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
!>          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(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 DOUBLE PRECISION array, dimension (LDZ,N)
!>          If WANTZ is .TRUE., on entry Z must contain the current
!>          matrix Z of transformations accumulated by DHSEQR, 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, DLAHQR failed to compute all the
!>                  eigenvalues ILO to IHI in a total of 30 iterations
!>                  per eigenvalue; elements i+1:ihi of WR and WI
!>                  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.
Further Details:
!>
!>     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 DLAHQR 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 203 of file dlahqr.f.

205 IMPLICIT NONE
206*
207* -- LAPACK auxiliary routine --
208* -- LAPACK is a software package provided by Univ. of Tennessee, --
209* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210*
211* .. Scalar Arguments ..
212 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
213 LOGICAL WANTT, WANTZ
214* ..
215* .. Array Arguments ..
216 DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
217* ..
218*
219* =========================================================
220*
221* .. Parameters ..
222 DOUBLE PRECISION ZERO, ONE, TWO
223 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
224 DOUBLE PRECISION DAT1, DAT2
225 parameter( dat1 = 3.0d0 / 4.0d0, dat2 = -0.4375d0 )
226 INTEGER KEXSH
227 parameter( kexsh = 10 )
228* ..
229* .. Local Scalars ..
230 DOUBLE PRECISION AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
231 $ H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX,
232 $ SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST,
233 $ ULP, V2, V3
234 INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ,
235 $ KDEFL
236* ..
237* .. Local Arrays ..
238 DOUBLE PRECISION V( 3 )
239* ..
240* .. External Functions ..
241 DOUBLE PRECISION DLAMCH
242 EXTERNAL dlamch
243* ..
244* .. External Subroutines ..
245 EXTERNAL dcopy, dlanv2, dlarfg, drot
246* ..
247* .. Intrinsic Functions ..
248 INTRINSIC abs, dble, max, min, sqrt
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 wr( ilo ) = h( ilo, ilo )
260 wi( ilo ) = zero
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*
272 nh = ihi - ilo + 1
273 nz = ihiz - iloz + 1
274*
275* Set machine-dependent constants for the stopping criterion.
276*
277 safmin = dlamch( 'SAFE MINIMUM' )
278 safmax = one / safmin
279 ulp = dlamch( 'PRECISION' )
280 smlnum = safmin*( dble( nh ) / ulp )
281*
282* I1 and I2 are the indices of the first row and last column of H
283* to which transformations must be applied. If eigenvalues only are
284* being computed, I1 and I2 are set inside the main loop.
285*
286 IF( wantt ) THEN
287 i1 = 1
288 i2 = n
289 END IF
290*
291* ITMAX is the total number of QR iterations allowed.
292*
293 itmax = 30 * max( 10, nh )
294*
295* KDEFL counts the number of iterations since a deflation
296*
297 kdefl = 0
298*
299* The main loop begins here. I is the loop index and decreases from
300* IHI to ILO in steps of 1 or 2. Each iteration of the loop works
301* with the active submatrix in rows and columns L to I.
302* Eigenvalues I+1 to IHI have already converged. Either L = ILO or
303* H(L,L-1) is negligible so that the matrix splits.
304*
305 i = ihi
306 20 CONTINUE
307 l = ilo
308 IF( i.LT.ilo )
309 $ GO TO 160
310*
311* Perform QR iterations on rows and columns ILO to I until a
312* submatrix of order 1 or 2 splits off at the bottom because a
313* subdiagonal element has become negligible.
314*
315 DO 140 its = 0, itmax
316*
317* Look for a single small subdiagonal element.
318*
319 DO 30 k = i, l + 1, -1
320 IF( abs( h( k, k-1 ) ).LE.smlnum )
321 $ GO TO 40
322 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
323 IF( tst.EQ.zero ) THEN
324 IF( k-2.GE.ilo )
325 $ tst = tst + abs( h( k-1, k-2 ) )
326 IF( k+1.LE.ihi )
327 $ tst = tst + abs( h( k+1, k ) )
328 END IF
329* ==== The following is a conservative small subdiagonal
330* . deflation criterion due to Ahues & Tisseur (LAWN 122,
331* . 1997). It has better mathematical foundation and
332* . improves accuracy in some cases. ====
333 IF( abs( h( k, k-1 ) ).LE.ulp*tst ) THEN
334 ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
335 ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
336 aa = max( abs( h( k, k ) ),
337 $ abs( h( k-1, k-1 )-h( k, k ) ) )
338 bb = min( abs( h( k, k ) ),
339 $ abs( h( k-1, k-1 )-h( k, k ) ) )
340 s = aa + ab
341 IF( ba*( ab / s ).LE.max( smlnum,
342 $ ulp*( bb*( aa / s ) ) ) )GO TO 40
343 END IF
344 30 CONTINUE
345 40 CONTINUE
346 l = k
347 IF( l.GT.ilo ) THEN
348*
349* H(L,L-1) is negligible
350*
351 h( l, l-1 ) = zero
352 END IF
353*
354* Exit from loop if a submatrix of order 1 or 2 has split off.
355*
356 IF( l.GE.i-1 )
357 $ GO TO 150
358 kdefl = kdefl + 1
359*
360* Now the active submatrix is in rows and columns L to I. If
361* eigenvalues only are being computed, only the active submatrix
362* need be transformed.
363*
364 IF( .NOT.wantt ) THEN
365 i1 = l
366 i2 = i
367 END IF
368*
369 IF( mod(kdefl,2*kexsh).EQ.0 ) THEN
370*
371* Exceptional shift.
372*
373 s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
374 h11 = dat1*s + h( i, i )
375 h12 = dat2*s
376 h21 = s
377 h22 = h11
378 ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
379*
380* Exceptional shift.
381*
382 s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
383 h11 = dat1*s + h( l, l )
384 h12 = dat2*s
385 h21 = s
386 h22 = h11
387 ELSE
388*
389* Prepare to use Francis' double shift
390* (i.e. 2nd degree generalized Rayleigh quotient)
391*
392 h11 = h( i-1, i-1 )
393 h21 = h( i, i-1 )
394 h12 = h( i-1, i )
395 h22 = h( i, i )
396 END IF
397 s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
398 IF( s.EQ.zero ) THEN
399 rt1r = zero
400 rt1i = zero
401 rt2r = zero
402 rt2i = zero
403 ELSE
404 h11 = h11 / s
405 h21 = h21 / s
406 h12 = h12 / s
407 h22 = h22 / s
408 tr = ( h11+h22 ) / two
409 det = ( h11-tr )*( h22-tr ) - h12*h21
410 rtdisc = sqrt( abs( det ) )
411 IF( det.GE.zero ) THEN
412*
413* ==== complex conjugate shifts ====
414*
415 rt1r = tr*s
416 rt2r = rt1r
417 rt1i = rtdisc*s
418 rt2i = -rt1i
419 ELSE
420*
421* ==== real shifts (use only one of them) ====
422*
423 rt1r = tr + rtdisc
424 rt2r = tr - rtdisc
425 IF( abs( rt1r-h22 ).LE.abs( rt2r-h22 ) ) THEN
426 rt1r = rt1r*s
427 rt2r = rt1r
428 ELSE
429 rt2r = rt2r*s
430 rt1r = rt2r
431 END IF
432 rt1i = zero
433 rt2i = zero
434 END IF
435 END IF
436*
437* Look for two consecutive small subdiagonal elements.
438*
439 DO 50 m = i - 2, l, -1
440* Determine the effect of starting the double-shift QR
441* iteration at row M, and see if this would make H(M,M-1)
442* negligible. (The following uses scaling to avoid
443* overflows and most underflows.)
444*
445 h21s = h( m+1, m )
446 s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
447 h21s = h( m+1, m ) / s
448 v( 1 ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*
449 $ ( ( h( m, m )-rt2r ) / s ) - rt1i*( rt2i / s )
450 v( 2 ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
451 v( 3 ) = h21s*h( m+2, m+1 )
452 s = abs( v( 1 ) ) + abs( v( 2 ) ) + abs( v( 3 ) )
453 v( 1 ) = v( 1 ) / s
454 v( 2 ) = v( 2 ) / s
455 v( 3 ) = v( 3 ) / s
456 IF( m.EQ.l )
457 $ GO TO 60
458 IF( abs( h( m, m-1 ) )*( abs( v( 2 ) )+abs( v( 3 ) ) ).LE.
459 $ ulp*abs( v( 1 ) )*( abs( h( m-1, m-1 ) )+abs( h( m,
460 $ m ) )+abs( h( m+1, m+1 ) ) ) )GO TO 60
461 50 CONTINUE
462 60 CONTINUE
463*
464* Double-shift QR step
465*
466 DO 130 k = m, i - 1
467*
468* The first iteration of this loop determines a reflection G
469* from the vector V and applies it from left and right to H,
470* thus creating a nonzero bulge below the subdiagonal.
471*
472* Each subsequent iteration determines a reflection G to
473* restore the Hessenberg form in the (K-1)th column, and thus
474* chases the bulge one step toward the bottom of the active
475* submatrix. NR is the order of G.
476*
477 nr = min( 3, i-k+1 )
478 IF( k.GT.m )
479 $ CALL dcopy( nr, h( k, k-1 ), 1, v, 1 )
480 CALL dlarfg( nr, v( 1 ), v( 2 ), 1, t1 )
481 IF( k.GT.m ) THEN
482 h( k, k-1 ) = v( 1 )
483 h( k+1, k-1 ) = zero
484 IF( k.LT.i-1 )
485 $ h( k+2, k-1 ) = zero
486 ELSE IF( m.GT.l ) THEN
487* ==== Use the following instead of
488* . H( K, K-1 ) = -H( K, K-1 ) to
489* . avoid a bug when v(2) and v(3)
490* . underflow. ====
491 h( k, k-1 ) = h( k, k-1 )*( one-t1 )
492 END IF
493 v2 = v( 2 )
494 t2 = t1*v2
495 IF( nr.EQ.3 ) THEN
496 v3 = v( 3 )
497 t3 = t1*v3
498*
499* Apply G from the left to transform the rows of the matrix
500* in columns K to I2.
501*
502 DO 70 j = k, i2
503 sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
504 h( k, j ) = h( k, j ) - sum*t1
505 h( k+1, j ) = h( k+1, j ) - sum*t2
506 h( k+2, j ) = h( k+2, j ) - sum*t3
507 70 CONTINUE
508*
509* Apply G from the right to transform the columns of the
510* matrix in rows I1 to min(K+3,I).
511*
512 DO 80 j = i1, min( k+3, i )
513 sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
514 h( j, k ) = h( j, k ) - sum*t1
515 h( j, k+1 ) = h( j, k+1 ) - sum*t2
516 h( j, k+2 ) = h( j, k+2 ) - sum*t3
517 80 CONTINUE
518*
519 IF( wantz ) THEN
520*
521* Accumulate transformations in the matrix Z
522*
523 DO 90 j = iloz, ihiz
524 sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
525 z( j, k ) = z( j, k ) - sum*t1
526 z( j, k+1 ) = z( j, k+1 ) - sum*t2
527 z( j, k+2 ) = z( j, k+2 ) - sum*t3
528 90 CONTINUE
529 END IF
530 ELSE IF( nr.EQ.2 ) THEN
531*
532* Apply G from the left to transform the rows of the matrix
533* in columns K to I2.
534*
535 DO 100 j = k, i2
536 sum = h( k, j ) + v2*h( k+1, j )
537 h( k, j ) = h( k, j ) - sum*t1
538 h( k+1, j ) = h( k+1, j ) - sum*t2
539 100 CONTINUE
540*
541* Apply G from the right to transform the columns of the
542* matrix in rows I1 to min(K+3,I).
543*
544 DO 110 j = i1, i
545 sum = h( j, k ) + v2*h( j, k+1 )
546 h( j, k ) = h( j, k ) - sum*t1
547 h( j, k+1 ) = h( j, k+1 ) - sum*t2
548 110 CONTINUE
549*
550 IF( wantz ) THEN
551*
552* Accumulate transformations in the matrix Z
553*
554 DO 120 j = iloz, ihiz
555 sum = z( j, k ) + v2*z( j, k+1 )
556 z( j, k ) = z( j, k ) - sum*t1
557 z( j, k+1 ) = z( j, k+1 ) - sum*t2
558 120 CONTINUE
559 END IF
560 END IF
561 130 CONTINUE
562*
563 140 CONTINUE
564*
565* Failure to converge in remaining number of iterations
566*
567 info = i
568 RETURN
569*
570 150 CONTINUE
571*
572 IF( l.EQ.i ) THEN
573*
574* H(I,I-1) is negligible: one eigenvalue has converged.
575*
576 wr( i ) = h( i, i )
577 wi( i ) = zero
578 ELSE IF( l.EQ.i-1 ) THEN
579*
580* H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
581*
582* Transform the 2-by-2 submatrix to standard Schur form,
583* and compute and store the eigenvalues.
584*
585 CALL dlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
586 $ h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
587 $ cs, sn )
588*
589 IF( wantt ) THEN
590*
591* Apply the transformation to the rest of H.
592*
593 IF( i2.GT.i )
594 $ CALL drot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
595 $ cs, sn )
596 CALL drot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs,
597 $ sn )
598 END IF
599 IF( wantz ) THEN
600*
601* Apply the transformation to Z.
602*
603 CALL drot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs,
604 $ sn )
605 END IF
606 END IF
607* reset deflation counter
608 kdefl = 0
609*
610* return to start of the main loop with new value of I.
611*
612 i = l - 1
613 GO TO 20
614*
615 160 CONTINUE
616 RETURN
617*
618* End of DLAHQR
619*
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanv2
subroutine dlanv2(a, b, c, d, rt1r, rt1i, rt2r, rt2i, cs, sn)
DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition dlanv2.f:125
dlarfg
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:104
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: