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

◆ dhgeqz()

subroutine dhgeqz ( character job,
character compq,
character compz,
integer n,
integer ilo,
integer ihi,
double precision, dimension( ldh, * ) h,
integer ldh,
double precision, dimension( ldt, * ) t,
integer ldt,
double precision, dimension( * ) alphar,
double precision, dimension( * ) alphai,
double precision, dimension( * ) beta,
double precision, dimension( ldq, * ) q,
integer ldq,
double precision, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer lwork,
integer info )

DHGEQZ

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

Purpose:
!>
!> DHGEQZ computes the eigenvalues of a real matrix pair (H,T),
!> where H is an upper Hessenberg matrix and T is upper triangular,
!> using the double-shift QZ method.
!> Matrix pairs of this type are produced by the reduction to
!> generalized upper Hessenberg form of a real matrix pair (A,B):
!>
!>    A = Q1*H*Z1**T,  B = Q1*T*Z1**T,
!>
!> as computed by DGGHRD.
!>
!> If JOB='S', then the Hessenberg-triangular pair (H,T) is
!> also reduced to generalized Schur form,
!>
!>    H = Q*S*Z**T,  T = Q*P*Z**T,
!>
!> where Q and Z are orthogonal matrices, P is an upper triangular
!> matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
!> diagonal blocks.
!>
!> The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
!> (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
!> eigenvalues.
!>
!> Additionally, the 2-by-2 upper triangular diagonal blocks of P
!> corresponding to 2-by-2 blocks of S are reduced to positive diagonal
!> form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
!> P(j,j) > 0, and P(j+1,j+1) > 0.
!>
!> Optionally, the orthogonal matrix Q from the generalized Schur
!> factorization may be postmultiplied into an input matrix Q1, and the
!> orthogonal matrix Z may be postmultiplied into an input matrix Z1.
!> If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced
!> the matrix pair (A,B) to generalized upper Hessenberg form, then the
!> output matrices Q1*Q and Z1*Z are the orthogonal factors from the
!> generalized Schur factorization of (A,B):
!>
!>    A = (Q1*Q)*S*(Z1*Z)**T,  B = (Q1*Q)*P*(Z1*Z)**T.
!>
!> To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
!> of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
!> complex and beta real.
!> If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
!> generalized nonsymmetric eigenvalue problem (GNEP)
!>    A*x = lambda*B*x
!> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
!> alternate form of the GNEP
!>    mu*A*y = B*y.
!> Real eigenvalues can be read directly from the generalized Schur
!> form:
!>   alpha = S(i,i), beta = P(i,i).
!>
!> Ref: C.B. Moler & G.W. Stewart, , SIAM J. Numer. Anal., 10(1973),
!>      pp. 241--256.
!>
Parameters
[in]JOB
!>          JOB is CHARACTER*1
!>          = 'E': Compute eigenvalues only;
!>          = 'S': Compute eigenvalues and the Schur form.
!>
[in]COMPQ
!>          COMPQ is CHARACTER*1
!>          = 'N': Left Schur vectors (Q) are not computed;
!>          = 'I': Q is initialized to the unit matrix and the matrix Q
!>                 of left Schur vectors of (H,T) is returned;
!>          = 'V': Q must contain an orthogonal matrix Q1 on entry and
!>                 the product Q1*Q is returned.
!>
[in]COMPZ
!>          COMPZ is CHARACTER*1
!>          = 'N': Right Schur vectors (Z) are not computed;
!>          = 'I': Z is initialized to the unit matrix and the matrix Z
!>                 of right Schur vectors of (H,T) is returned;
!>          = 'V': Z must contain an orthogonal matrix Z1 on entry and
!>                 the product Z1*Z is returned.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices H, T, Q, and Z.  N >= 0.
!>
[in]ILO
!>          ILO is INTEGER
!>
[in]IHI
!>          IHI is INTEGER
!>          ILO and IHI mark the rows and columns of H which are in
!>          Hessenberg form.  It is assumed that A is already upper
!>          triangular in rows and columns 1:ILO-1 and IHI+1:N.
!>          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
!>
[in,out]H
!>          H is DOUBLE PRECISION array, dimension (LDH, N)
!>          On entry, the N-by-N upper Hessenberg matrix H.
!>          On exit, if JOB = 'S', H contains the upper quasi-triangular
!>          matrix S from the generalized Schur factorization.
!>          If JOB = 'E', the diagonal blocks of H match those of S, but
!>          the rest of H is unspecified.
!>
[in]LDH
!>          LDH is INTEGER
!>          The leading dimension of the array H.  LDH >= max( 1, N ).
!>
[in,out]T
!>          T is DOUBLE PRECISION array, dimension (LDT, N)
!>          On entry, the N-by-N upper triangular matrix T.
!>          On exit, if JOB = 'S', T contains the upper triangular
!>          matrix P from the generalized Schur factorization;
!>          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
!>          are reduced to positive diagonal form, i.e., if H(j+1,j) is
!>          non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
!>          T(j+1,j+1) > 0.
!>          If JOB = 'E', the diagonal blocks of T match those of P, but
!>          the rest of T is unspecified.
!>
[in]LDT
!>          LDT is INTEGER
!>          The leading dimension of the array T.  LDT >= max( 1, N ).
!>
[out]ALPHAR
!>          ALPHAR is DOUBLE PRECISION array, dimension (N)
!>          The real parts of each scalar alpha defining an eigenvalue
!>          of GNEP.
!>
[out]ALPHAI
!>          ALPHAI is DOUBLE PRECISION array, dimension (N)
!>          The imaginary parts of each scalar alpha defining an
!>          eigenvalue of GNEP.
!>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
!>          positive, then the j-th and (j+1)-st eigenvalues are a
!>          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>          The scalars beta that define the eigenvalues of GNEP.
!>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
!>          beta = BETA(j) represent the j-th eigenvalue of the matrix
!>          pair (A,B), in one of the forms lambda = alpha/beta or
!>          mu = beta/alpha.  Since either lambda or mu may overflow,
!>          they should not, in general, be computed.
!>
[in,out]Q
!>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
!>          On entry, if COMPQ = 'V', the orthogonal matrix Q1 used in
!>          the reduction of (A,B) to generalized Hessenberg form.
!>          On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
!>          vectors of (H,T), and if COMPQ = 'V', the orthogonal matrix
!>          of left Schur vectors of (A,B).
!>          Not referenced if COMPQ = 'N'.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q.  LDQ >= 1.
!>          If COMPQ='V' or 'I', then LDQ >= N.
!>
[in,out]Z
!>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
!>          On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
!>          the reduction of (A,B) to generalized Hessenberg form.
!>          On exit, if COMPZ = 'I', the orthogonal matrix of
!>          right Schur vectors of (H,T), and if COMPZ = 'V', the
!>          orthogonal matrix of right Schur vectors of (A,B).
!>          Not referenced if COMPZ = 'N'.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1.
!>          If COMPZ='V' or 'I', then LDZ >= N.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,N).
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
!>                     in Schur form, but ALPHAR(i), ALPHAI(i), and
!>                     BETA(i), i=INFO+1,...,N should be correct.
!>          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
!>                     in Schur form, but ALPHAR(i), ALPHAI(i), and
!>                     BETA(i), i=INFO-N+1,...,N should be correct.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Iteration counters:
!>
!>  JITER  -- counts iterations.
!>  IITER  -- counts iterations run since ILAST was last
!>            changed.  This is therefore reset only when a 1-by-1 or
!>            2-by-2 block deflates off the bottom.
!>

Definition at line 299 of file dhgeqz.f.

303*
304* -- LAPACK computational routine --
305* -- LAPACK is a software package provided by Univ. of Tennessee, --
306* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
307*
308* .. Scalar Arguments ..
309 CHARACTER COMPQ, COMPZ, JOB
310 INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
311* ..
312* .. Array Arguments ..
313 DOUBLE PRECISION ALPHAI( * ), ALPHAR( * ), BETA( * ),
314 $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
315 $ WORK( * ), Z( LDZ, * )
316* ..
317*
318* =====================================================================
319*
320* .. Parameters ..
321* $ SAFETY = 1.0E+0 )
322 DOUBLE PRECISION HALF, ZERO, ONE, SAFETY
323 parameter( half = 0.5d+0, zero = 0.0d+0, one = 1.0d+0,
324 $ safety = 1.0d+2 )
325* ..
326* .. Local Scalars ..
327 LOGICAL ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
328 $ LQUERY
329 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
330 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
331 $ JR, MAXIT
332 DOUBLE PRECISION A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
333 $ AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
334 $ AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I,
335 $ B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE,
336 $ BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL,
337 $ CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX,
338 $ SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1,
339 $ T2, T3, TAU, TEMP, TEMP2, TEMPI, TEMPR, U1,
340 $ U12, U12L, U2, ULP, VS, W11, W12, W21, W22,
341 $ WABS, WI, WR, WR2
342* ..
343* .. Local Arrays ..
344 DOUBLE PRECISION V( 3 )
345* ..
346* .. External Functions ..
347 LOGICAL LSAME
348 DOUBLE PRECISION DLAMCH, DLANHS, DLAPY2, DLAPY3
349 EXTERNAL lsame, dlamch, dlanhs, dlapy2,
350 $ dlapy3
351* ..
352* .. External Subroutines ..
353 EXTERNAL dlag2, dlarfg, dlartg, dlaset, dlasv2,
354 $ drot,
355 $ xerbla
356* ..
357* .. Intrinsic Functions ..
358 INTRINSIC abs, dble, max, min, sqrt
359* ..
360* .. Executable Statements ..
361*
362* Decode JOB, COMPQ, COMPZ
363*
364 IF( lsame( job, 'E' ) ) THEN
365 ilschr = .false.
366 ischur = 1
367 ELSE IF( lsame( job, 'S' ) ) THEN
368 ilschr = .true.
369 ischur = 2
370 ELSE
371 ischur = 0
372 END IF
373*
374 IF( lsame( compq, 'N' ) ) THEN
375 ilq = .false.
376 icompq = 1
377 ELSE IF( lsame( compq, 'V' ) ) THEN
378 ilq = .true.
379 icompq = 2
380 ELSE IF( lsame( compq, 'I' ) ) THEN
381 ilq = .true.
382 icompq = 3
383 ELSE
384 icompq = 0
385 END IF
386*
387 IF( lsame( compz, 'N' ) ) THEN
388 ilz = .false.
389 icompz = 1
390 ELSE IF( lsame( compz, 'V' ) ) THEN
391 ilz = .true.
392 icompz = 2
393 ELSE IF( lsame( compz, 'I' ) ) THEN
394 ilz = .true.
395 icompz = 3
396 ELSE
397 icompz = 0
398 END IF
399*
400* Check Argument Values
401*
402 info = 0
403 work( 1 ) = max( 1, n )
404 lquery = ( lwork.EQ.-1 )
405 IF( ischur.EQ.0 ) THEN
406 info = -1
407 ELSE IF( icompq.EQ.0 ) THEN
408 info = -2
409 ELSE IF( icompz.EQ.0 ) THEN
410 info = -3
411 ELSE IF( n.LT.0 ) THEN
412 info = -4
413 ELSE IF( ilo.LT.1 ) THEN
414 info = -5
415 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
416 info = -6
417 ELSE IF( ldh.LT.n ) THEN
418 info = -8
419 ELSE IF( ldt.LT.n ) THEN
420 info = -10
421 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
422 info = -15
423 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
424 info = -17
425 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
426 info = -19
427 END IF
428 IF( info.NE.0 ) THEN
429 CALL xerbla( 'DHGEQZ', -info )
430 RETURN
431 ELSE IF( lquery ) THEN
432 RETURN
433 END IF
434*
435* Quick return if possible
436*
437 IF( n.LE.0 ) THEN
438 work( 1 ) = dble( 1 )
439 RETURN
440 END IF
441*
442* Initialize Q and Z
443*
444 IF( icompq.EQ.3 )
445 $ CALL dlaset( 'Full', n, n, zero, one, q, ldq )
446 IF( icompz.EQ.3 )
447 $ CALL dlaset( 'Full', n, n, zero, one, z, ldz )
448*
449* Machine Constants
450*
451 in = ihi + 1 - ilo
452 safmin = dlamch( 'S' )
453 safmax = one / safmin
454 ulp = dlamch( 'E' )*dlamch( 'B' )
455 anorm = dlanhs( 'F', in, h( ilo, ilo ), ldh, work )
456 bnorm = dlanhs( 'F', in, t( ilo, ilo ), ldt, work )
457 atol = max( safmin, ulp*anorm )
458 btol = max( safmin, ulp*bnorm )
459 ascale = one / max( safmin, anorm )
460 bscale = one / max( safmin, bnorm )
461*
462* Set Eigenvalues IHI+1:N
463*
464 DO 30 j = ihi + 1, n
465 IF( t( j, j ).LT.zero ) THEN
466 IF( ilschr ) THEN
467 DO 10 jr = 1, j
468 h( jr, j ) = -h( jr, j )
469 t( jr, j ) = -t( jr, j )
470 10 CONTINUE
471 ELSE
472 h( j, j ) = -h( j, j )
473 t( j, j ) = -t( j, j )
474 END IF
475 IF( ilz ) THEN
476 DO 20 jr = 1, n
477 z( jr, j ) = -z( jr, j )
478 20 CONTINUE
479 END IF
480 END IF
481 alphar( j ) = h( j, j )
482 alphai( j ) = zero
483 beta( j ) = t( j, j )
484 30 CONTINUE
485*
486* If IHI < ILO, skip QZ steps
487*
488 IF( ihi.LT.ilo )
489 $ GO TO 380
490*
491* MAIN QZ ITERATION LOOP
492*
493* Initialize dynamic indices
494*
495* Eigenvalues ILAST+1:N have been found.
496* Column operations modify rows IFRSTM:whatever.
497* Row operations modify columns whatever:ILASTM.
498*
499* If only eigenvalues are being computed, then
500* IFRSTM is the row of the last splitting row above row ILAST;
501* this is always at least ILO.
502* IITER counts iterations since the last eigenvalue was found,
503* to tell when to use an extraordinary shift.
504* MAXIT is the maximum number of QZ sweeps allowed.
505*
506 ilast = ihi
507 IF( ilschr ) THEN
508 ifrstm = 1
509 ilastm = n
510 ELSE
511 ifrstm = ilo
512 ilastm = ihi
513 END IF
514 iiter = 0
515 eshift = zero
516 maxit = 30*( ihi-ilo+1 )
517*
518 DO 360 jiter = 1, maxit
519*
520* Split the matrix if possible.
521*
522* Two tests:
523* 1: H(j,j-1)=0 or j=ILO
524* 2: T(j,j)=0
525*
526 IF( ilast.EQ.ilo ) THEN
527*
528* Special case: j=ILAST
529*
530 GO TO 80
531 ELSE
532 IF( abs( h( ilast, ilast-1 ) ).LE.max( safmin, ulp*(
533 $ abs( h( ilast, ilast ) ) + abs( h( ilast-1, ilast-1 ) )
534 $ ) ) ) THEN
535 h( ilast, ilast-1 ) = zero
536 GO TO 80
537 END IF
538 END IF
539*
540 IF( abs( t( ilast, ilast ) ).LE.btol ) THEN
541 t( ilast, ilast ) = zero
542 GO TO 70
543 END IF
544*
545* General case: j<ILAST
546*
547 DO 60 j = ilast - 1, ilo, -1
548*
549* Test 1: for H(j,j-1)=0 or j=ILO
550*
551 IF( j.EQ.ilo ) THEN
552 ilazro = .true.
553 ELSE
554 IF( abs( h( j, j-1 ) ).LE.max( safmin, ulp*(
555 $ abs( h( j, j ) ) + abs( h( j-1, j-1 ) )
556 $ ) ) ) THEN
557 h( j, j-1 ) = zero
558 ilazro = .true.
559 ELSE
560 ilazro = .false.
561 END IF
562 END IF
563*
564* Test 2: for T(j,j)=0
565*
566 IF( abs( t( j, j ) ).LT.btol ) THEN
567 t( j, j ) = zero
568*
569* Test 1a: Check for 2 consecutive small subdiagonals in A
570*
571 ilazr2 = .false.
572 IF( .NOT.ilazro ) THEN
573 temp = abs( h( j, j-1 ) )
574 temp2 = abs( h( j, j ) )
575 tempr = max( temp, temp2 )
576 IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
577 temp = temp / tempr
578 temp2 = temp2 / tempr
579 END IF
580 IF( temp*( ascale*abs( h( j+1, j ) ) ).LE.temp2*
581 $ ( ascale*atol ) )ilazr2 = .true.
582 END IF
583*
584* If both tests pass (1 & 2), i.e., the leading diagonal
585* element of B in the block is zero, split a 1x1 block off
586* at the top. (I.e., at the J-th row/column) The leading
587* diagonal element of the remainder can also be zero, so
588* this may have to be done repeatedly.
589*
590 IF( ilazro .OR. ilazr2 ) THEN
591 DO 40 jch = j, ilast - 1
592 temp = h( jch, jch )
593 CALL dlartg( temp, h( jch+1, jch ), c, s,
594 $ h( jch, jch ) )
595 h( jch+1, jch ) = zero
596 CALL drot( ilastm-jch, h( jch, jch+1 ), ldh,
597 $ h( jch+1, jch+1 ), ldh, c, s )
598 CALL drot( ilastm-jch, t( jch, jch+1 ), ldt,
599 $ t( jch+1, jch+1 ), ldt, c, s )
600 IF( ilq )
601 $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ),
602 $ 1,
603 $ c, s )
604 IF( ilazr2 )
605 $ h( jch, jch-1 ) = h( jch, jch-1 )*c
606 ilazr2 = .false.
607 IF( abs( t( jch+1, jch+1 ) ).GE.btol ) THEN
608 IF( jch+1.GE.ilast ) THEN
609 GO TO 80
610 ELSE
611 ifirst = jch + 1
612 GO TO 110
613 END IF
614 END IF
615 t( jch+1, jch+1 ) = zero
616 40 CONTINUE
617 GO TO 70
618 ELSE
619*
620* Only test 2 passed -- chase the zero to T(ILAST,ILAST)
621* Then process as in the case T(ILAST,ILAST)=0
622*
623 DO 50 jch = j, ilast - 1
624 temp = t( jch, jch+1 )
625 CALL dlartg( temp, t( jch+1, jch+1 ), c, s,
626 $ t( jch, jch+1 ) )
627 t( jch+1, jch+1 ) = zero
628 IF( jch.LT.ilastm-1 )
629 $ CALL drot( ilastm-jch-1, t( jch, jch+2 ),
630 $ ldt,
631 $ t( jch+1, jch+2 ), ldt, c, s )
632 CALL drot( ilastm-jch+2, h( jch, jch-1 ), ldh,
633 $ h( jch+1, jch-1 ), ldh, c, s )
634 IF( ilq )
635 $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ),
636 $ 1,
637 $ c, s )
638 temp = h( jch+1, jch )
639 CALL dlartg( temp, h( jch+1, jch-1 ), c, s,
640 $ h( jch+1, jch ) )
641 h( jch+1, jch-1 ) = zero
642 CALL drot( jch+1-ifrstm, h( ifrstm, jch ), 1,
643 $ h( ifrstm, jch-1 ), 1, c, s )
644 CALL drot( jch-ifrstm, t( ifrstm, jch ), 1,
645 $ t( ifrstm, jch-1 ), 1, c, s )
646 IF( ilz )
647 $ CALL drot( n, z( 1, jch ), 1, z( 1, jch-1 ),
648 $ 1,
649 $ c, s )
650 50 CONTINUE
651 GO TO 70
652 END IF
653 ELSE IF( ilazro ) THEN
654*
655* Only test 1 passed -- work on J:ILAST
656*
657 ifirst = j
658 GO TO 110
659 END IF
660*
661* Neither test passed -- try next J
662*
663 60 CONTINUE
664*
665* (Drop-through is "impossible")
666*
667 info = n + 1
668 GO TO 420
669*
670* T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
671* 1x1 block.
672*
673 70 CONTINUE
674 temp = h( ilast, ilast )
675 CALL dlartg( temp, h( ilast, ilast-1 ), c, s,
676 $ h( ilast, ilast ) )
677 h( ilast, ilast-1 ) = zero
678 CALL drot( ilast-ifrstm, h( ifrstm, ilast ), 1,
679 $ h( ifrstm, ilast-1 ), 1, c, s )
680 CALL drot( ilast-ifrstm, t( ifrstm, ilast ), 1,
681 $ t( ifrstm, ilast-1 ), 1, c, s )
682 IF( ilz )
683 $ CALL drot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c,
684 $ s )
685*
686* H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
687* and BETA
688*
689 80 CONTINUE
690 IF( t( ilast, ilast ).LT.zero ) THEN
691 IF( ilschr ) THEN
692 DO 90 j = ifrstm, ilast
693 h( j, ilast ) = -h( j, ilast )
694 t( j, ilast ) = -t( j, ilast )
695 90 CONTINUE
696 ELSE
697 h( ilast, ilast ) = -h( ilast, ilast )
698 t( ilast, ilast ) = -t( ilast, ilast )
699 END IF
700 IF( ilz ) THEN
701 DO 100 j = 1, n
702 z( j, ilast ) = -z( j, ilast )
703 100 CONTINUE
704 END IF
705 END IF
706 alphar( ilast ) = h( ilast, ilast )
707 alphai( ilast ) = zero
708 beta( ilast ) = t( ilast, ilast )
709*
710* Go to next block -- exit if finished.
711*
712 ilast = ilast - 1
713 IF( ilast.LT.ilo )
714 $ GO TO 380
715*
716* Reset counters
717*
718 iiter = 0
719 eshift = zero
720 IF( .NOT.ilschr ) THEN
721 ilastm = ilast
722 IF( ifrstm.GT.ilast )
723 $ ifrstm = ilo
724 END IF
725 GO TO 350
726*
727* QZ step
728*
729* This iteration only involves rows/columns IFIRST:ILAST. We
730* assume IFIRST < ILAST, and that the diagonal of B is non-zero.
731*
732 110 CONTINUE
733 iiter = iiter + 1
734 IF( .NOT.ilschr ) THEN
735 ifrstm = ifirst
736 END IF
737*
738* Compute single shifts.
739*
740* At this point, IFIRST < ILAST, and the diagonal elements of
741* T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
742* magnitude)
743*
744 IF( ( iiter / 10 )*10.EQ.iiter ) THEN
745*
746* Exceptional shift. Chosen for no particularly good reason.
747* (Single shift only.)
748*
749 IF( ( dble( maxit )*safmin )*abs( h( ilast, ilast-1 ) ).LT.
750 $ abs( t( ilast-1, ilast-1 ) ) ) THEN
751 eshift = h( ilast, ilast-1 ) /
752 $ t( ilast-1, ilast-1 )
753 ELSE
754 eshift = eshift + one / ( safmin*dble( maxit ) )
755 END IF
756 s1 = one
757 wr = eshift
758*
759 ELSE
760*
761* Shifts based on the generalized eigenvalues of the
762* bottom-right 2x2 block of A and B. The first eigenvalue
763* returned by DLAG2 is the Wilkinson shift (AEP p.512),
764*
765 CALL dlag2( h( ilast-1, ilast-1 ), ldh,
766 $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
767 $ s2, wr, wr2, wi )
768*
769 IF ( abs( (wr/s1)*t( ilast, ilast ) - h( ilast, ilast ) )
770 $ .GT. abs( (wr2/s2)*t( ilast, ilast )
771 $ - h( ilast, ilast ) ) ) THEN
772 temp = wr
773 wr = wr2
774 wr2 = temp
775 temp = s1
776 s1 = s2
777 s2 = temp
778 END IF
779 temp = max( s1, safmin*max( one, abs( wr ), abs( wi ) ) )
780 IF( wi.NE.zero )
781 $ GO TO 200
782 END IF
783*
784* Fiddle with shift to avoid overflow
785*
786 temp = min( ascale, one )*( half*safmax )
787 IF( s1.GT.temp ) THEN
788 scale = temp / s1
789 ELSE
790 scale = one
791 END IF
792*
793 temp = min( bscale, one )*( half*safmax )
794 IF( abs( wr ).GT.temp )
795 $ scale = min( scale, temp / abs( wr ) )
796 s1 = scale*s1
797 wr = scale*wr
798*
799* Now check for two consecutive small subdiagonals.
800*
801 DO 120 j = ilast - 1, ifirst + 1, -1
802 istart = j
803 temp = abs( s1*h( j, j-1 ) )
804 temp2 = abs( s1*h( j, j )-wr*t( j, j ) )
805 tempr = max( temp, temp2 )
806 IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
807 temp = temp / tempr
808 temp2 = temp2 / tempr
809 END IF
810 IF( abs( ( ascale*h( j+1, j ) )*temp ).LE.( ascale*atol )*
811 $ temp2 )GO TO 130
812 120 CONTINUE
813*
814 istart = ifirst
815 130 CONTINUE
816*
817* Do an implicit single-shift QZ sweep.
818*
819* Initial Q
820*
821 temp = s1*h( istart, istart ) - wr*t( istart, istart )
822 temp2 = s1*h( istart+1, istart )
823 CALL dlartg( temp, temp2, c, s, tempr )
824*
825* Sweep
826*
827 DO 190 j = istart, ilast - 1
828 IF( j.GT.istart ) THEN
829 temp = h( j, j-1 )
830 CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
831 h( j+1, j-1 ) = zero
832 END IF
833*
834 DO 140 jc = j, ilastm
835 temp = c*h( j, jc ) + s*h( j+1, jc )
836 h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
837 h( j, jc ) = temp
838 temp2 = c*t( j, jc ) + s*t( j+1, jc )
839 t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
840 t( j, jc ) = temp2
841 140 CONTINUE
842 IF( ilq ) THEN
843 DO 150 jr = 1, n
844 temp = c*q( jr, j ) + s*q( jr, j+1 )
845 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
846 q( jr, j ) = temp
847 150 CONTINUE
848 END IF
849*
850 temp = t( j+1, j+1 )
851 CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
852 t( j+1, j ) = zero
853*
854 DO 160 jr = ifrstm, min( j+2, ilast )
855 temp = c*h( jr, j+1 ) + s*h( jr, j )
856 h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
857 h( jr, j+1 ) = temp
858 160 CONTINUE
859 DO 170 jr = ifrstm, j
860 temp = c*t( jr, j+1 ) + s*t( jr, j )
861 t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
862 t( jr, j+1 ) = temp
863 170 CONTINUE
864 IF( ilz ) THEN
865 DO 180 jr = 1, n
866 temp = c*z( jr, j+1 ) + s*z( jr, j )
867 z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
868 z( jr, j+1 ) = temp
869 180 CONTINUE
870 END IF
871 190 CONTINUE
872*
873 GO TO 350
874*
875* Use Francis double-shift
876*
877* Note: the Francis double-shift should work with real shifts,
878* but only if the block is at least 3x3.
879* This code may break if this point is reached with
880* a 2x2 block with real eigenvalues.
881*
882 200 CONTINUE
883 IF( ifirst+1.EQ.ilast ) THEN
884*
885* Special case -- 2x2 block with complex eigenvectors
886*
887* Step 1: Standardize, that is, rotate so that
888*
889* ( B11 0 )
890* B = ( ) with B11 non-negative.
891* ( 0 B22 )
892*
893 CALL dlasv2( t( ilast-1, ilast-1 ), t( ilast-1, ilast ),
894 $ t( ilast, ilast ), b22, b11, sr, cr, sl, cl )
895*
896 IF( b11.LT.zero ) THEN
897 cr = -cr
898 sr = -sr
899 b11 = -b11
900 b22 = -b22
901 END IF
902*
903 CALL drot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
904 $ h( ilast, ilast-1 ), ldh, cl, sl )
905 CALL drot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
906 $ h( ifrstm, ilast ), 1, cr, sr )
907*
908 IF( ilast.LT.ilastm )
909 $ CALL drot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
910 $ t( ilast, ilast+1 ), ldt, cl, sl )
911 IF( ifrstm.LT.ilast-1 )
912 $ CALL drot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
913 $ t( ifrstm, ilast ), 1, cr, sr )
914*
915 IF( ilq )
916 $ CALL drot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1,
917 $ cl,
918 $ sl )
919 IF( ilz )
920 $ CALL drot( n, z( 1, ilast-1 ), 1, z( 1, ilast ), 1,
921 $ cr,
922 $ sr )
923*
924 t( ilast-1, ilast-1 ) = b11
925 t( ilast-1, ilast ) = zero
926 t( ilast, ilast-1 ) = zero
927 t( ilast, ilast ) = b22
928*
929* If B22 is negative, negate column ILAST
930*
931 IF( b22.LT.zero ) THEN
932 DO 210 j = ifrstm, ilast
933 h( j, ilast ) = -h( j, ilast )
934 t( j, ilast ) = -t( j, ilast )
935 210 CONTINUE
936*
937 IF( ilz ) THEN
938 DO 220 j = 1, n
939 z( j, ilast ) = -z( j, ilast )
940 220 CONTINUE
941 END IF
942 b22 = -b22
943 END IF
944*
945* Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
946*
947* Recompute shift
948*
949 CALL dlag2( h( ilast-1, ilast-1 ), ldh,
950 $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
951 $ temp, wr, temp2, wi )
952*
953* If standardization has perturbed the shift onto real line,
954* do another (real single-shift) QR step.
955*
956 IF( wi.EQ.zero )
957 $ GO TO 350
958 s1inv = one / s1
959*
960* Do EISPACK (QZVAL) computation of alpha and beta
961*
962 a11 = h( ilast-1, ilast-1 )
963 a21 = h( ilast, ilast-1 )
964 a12 = h( ilast-1, ilast )
965 a22 = h( ilast, ilast )
966*
967* Compute complex Givens rotation on right
968* (Assume some element of C = (sA - wB) > unfl )
969* __
970* (sA - wB) ( CZ -SZ )
971* ( SZ CZ )
972*
973 c11r = s1*a11 - wr*b11
974 c11i = -wi*b11
975 c12 = s1*a12
976 c21 = s1*a21
977 c22r = s1*a22 - wr*b22
978 c22i = -wi*b22
979*
980 IF( abs( c11r )+abs( c11i )+abs( c12 ).GT.abs( c21 )+
981 $ abs( c22r )+abs( c22i ) ) THEN
982 t1 = dlapy3( c12, c11r, c11i )
983 cz = c12 / t1
984 szr = -c11r / t1
985 szi = -c11i / t1
986 ELSE
987 cz = dlapy2( c22r, c22i )
988 IF( cz.LE.safmin ) THEN
989 cz = zero
990 szr = one
991 szi = zero
992 ELSE
993 tempr = c22r / cz
994 tempi = c22i / cz
995 t1 = dlapy2( cz, c21 )
996 cz = cz / t1
997 szr = -c21*tempr / t1
998 szi = c21*tempi / t1
999 END IF
1000 END IF
1001*
1002* Compute Givens rotation on left
1003*
1004* ( CQ SQ )
1005* ( __ ) A or B
1006* ( -SQ CQ )
1007*
1008 an = abs( a11 ) + abs( a12 ) + abs( a21 ) + abs( a22 )
1009 bn = abs( b11 ) + abs( b22 )
1010 wabs = abs( wr ) + abs( wi )
1011 IF( s1*an.GT.wabs*bn ) THEN
1012 cq = cz*b11
1013 sqr = szr*b22
1014 sqi = -szi*b22
1015 ELSE
1016 a1r = cz*a11 + szr*a12
1017 a1i = szi*a12
1018 a2r = cz*a21 + szr*a22
1019 a2i = szi*a22
1020 cq = dlapy2( a1r, a1i )
1021 IF( cq.LE.safmin ) THEN
1022 cq = zero
1023 sqr = one
1024 sqi = zero
1025 ELSE
1026 tempr = a1r / cq
1027 tempi = a1i / cq
1028 sqr = tempr*a2r + tempi*a2i
1029 sqi = tempi*a2r - tempr*a2i
1030 END IF
1031 END IF
1032 t1 = dlapy3( cq, sqr, sqi )
1033 cq = cq / t1
1034 sqr = sqr / t1
1035 sqi = sqi / t1
1036*
1037* Compute diagonal elements of QBZ
1038*
1039 tempr = sqr*szr - sqi*szi
1040 tempi = sqr*szi + sqi*szr
1041 b1r = cq*cz*b11 + tempr*b22
1042 b1i = tempi*b22
1043 b1a = dlapy2( b1r, b1i )
1044 b2r = cq*cz*b22 + tempr*b11
1045 b2i = -tempi*b11
1046 b2a = dlapy2( b2r, b2i )
1047*
1048* Normalize so beta > 0, and Im( alpha1 ) > 0
1049*
1050 beta( ilast-1 ) = b1a
1051 beta( ilast ) = b2a
1052 alphar( ilast-1 ) = ( wr*b1a )*s1inv
1053 alphai( ilast-1 ) = ( wi*b1a )*s1inv
1054 alphar( ilast ) = ( wr*b2a )*s1inv
1055 alphai( ilast ) = -( wi*b2a )*s1inv
1056*
1057* Step 3: Go to next block -- exit if finished.
1058*
1059 ilast = ifirst - 1
1060 IF( ilast.LT.ilo )
1061 $ GO TO 380
1062*
1063* Reset counters
1064*
1065 iiter = 0
1066 eshift = zero
1067 IF( .NOT.ilschr ) THEN
1068 ilastm = ilast
1069 IF( ifrstm.GT.ilast )
1070 $ ifrstm = ilo
1071 END IF
1072 GO TO 350
1073 ELSE
1074*
1075* Usual case: 3x3 or larger block, using Francis implicit
1076* double-shift
1077*
1078* 2
1079* Eigenvalue equation is w - c w + d = 0,
1080*
1081* -1 2 -1
1082* so compute 1st column of (A B ) - c A B + d
1083* using the formula in QZIT (from EISPACK)
1084*
1085* We assume that the block is at least 3x3
1086*
1087 ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
1088 $ ( bscale*t( ilast-1, ilast-1 ) )
1089 ad21 = ( ascale*h( ilast, ilast-1 ) ) /
1090 $ ( bscale*t( ilast-1, ilast-1 ) )
1091 ad12 = ( ascale*h( ilast-1, ilast ) ) /
1092 $ ( bscale*t( ilast, ilast ) )
1093 ad22 = ( ascale*h( ilast, ilast ) ) /
1094 $ ( bscale*t( ilast, ilast ) )
1095 u12 = t( ilast-1, ilast ) / t( ilast, ilast )
1096 ad11l = ( ascale*h( ifirst, ifirst ) ) /
1097 $ ( bscale*t( ifirst, ifirst ) )
1098 ad21l = ( ascale*h( ifirst+1, ifirst ) ) /
1099 $ ( bscale*t( ifirst, ifirst ) )
1100 ad12l = ( ascale*h( ifirst, ifirst+1 ) ) /
1101 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1102 ad22l = ( ascale*h( ifirst+1, ifirst+1 ) ) /
1103 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1104 ad32l = ( ascale*h( ifirst+2, ifirst+1 ) ) /
1105 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1106 u12l = t( ifirst, ifirst+1 ) / t( ifirst+1, ifirst+1 )
1107*
1108 v( 1 ) = ( ad11-ad11l )*( ad22-ad11l ) - ad12*ad21 +
1109 $ ad21*u12*ad11l + ( ad12l-ad11l*u12l )*ad21l
1110 v( 2 ) = ( ( ad22l-ad11l )-ad21l*u12l-( ad11-ad11l )-
1111 $ ( ad22-ad11l )+ad21*u12 )*ad21l
1112 v( 3 ) = ad32l*ad21l
1113*
1114 istart = ifirst
1115*
1116 CALL dlarfg( 3, v( 1 ), v( 2 ), 1, tau )
1117 v( 1 ) = one
1118*
1119* Sweep
1120*
1121 DO 290 j = istart, ilast - 2
1122*
1123* All but last elements: use 3x3 Householder transforms.
1124*
1125* Zero (j-1)st column of A
1126*
1127 IF( j.GT.istart ) THEN
1128 v( 1 ) = h( j, j-1 )
1129 v( 2 ) = h( j+1, j-1 )
1130 v( 3 ) = h( j+2, j-1 )
1131*
1132 CALL dlarfg( 3, h( j, j-1 ), v( 2 ), 1, tau )
1133 v( 1 ) = one
1134 h( j+1, j-1 ) = zero
1135 h( j+2, j-1 ) = zero
1136 END IF
1137*
1138 t2 = tau*v( 2 )
1139 t3 = tau*v( 3 )
1140 DO 230 jc = j, ilastm
1141 temp = h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
1142 $ h( j+2, jc )
1143 h( j, jc ) = h( j, jc ) - temp*tau
1144 h( j+1, jc ) = h( j+1, jc ) - temp*t2
1145 h( j+2, jc ) = h( j+2, jc ) - temp*t3
1146 temp2 = t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
1147 $ t( j+2, jc )
1148 t( j, jc ) = t( j, jc ) - temp2*tau
1149 t( j+1, jc ) = t( j+1, jc ) - temp2*t2
1150 t( j+2, jc ) = t( j+2, jc ) - temp2*t3
1151 230 CONTINUE
1152 IF( ilq ) THEN
1153 DO 240 jr = 1, n
1154 temp = q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
1155 $ q( jr, j+2 )
1156 q( jr, j ) = q( jr, j ) - temp*tau
1157 q( jr, j+1 ) = q( jr, j+1 ) - temp*t2
1158 q( jr, j+2 ) = q( jr, j+2 ) - temp*t3
1159 240 CONTINUE
1160 END IF
1161*
1162* Zero j-th column of B (see DLAGBC for details)
1163*
1164* Swap rows to pivot
1165*
1166 ilpivt = .false.
1167 temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
1168 temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
1169 IF( max( temp, temp2 ).LT.safmin ) THEN
1170 scale = zero
1171 u1 = one
1172 u2 = zero
1173 GO TO 250
1174 ELSE IF( temp.GE.temp2 ) THEN
1175 w11 = t( j+1, j+1 )
1176 w21 = t( j+2, j+1 )
1177 w12 = t( j+1, j+2 )
1178 w22 = t( j+2, j+2 )
1179 u1 = t( j+1, j )
1180 u2 = t( j+2, j )
1181 ELSE
1182 w21 = t( j+1, j+1 )
1183 w11 = t( j+2, j+1 )
1184 w22 = t( j+1, j+2 )
1185 w12 = t( j+2, j+2 )
1186 u2 = t( j+1, j )
1187 u1 = t( j+2, j )
1188 END IF
1189*
1190* Swap columns if nec.
1191*
1192 IF( abs( w12 ).GT.abs( w11 ) ) THEN
1193 ilpivt = .true.
1194 temp = w12
1195 temp2 = w22
1196 w12 = w11
1197 w22 = w21
1198 w11 = temp
1199 w21 = temp2
1200 END IF
1201*
1202* LU-factor
1203*
1204 temp = w21 / w11
1205 u2 = u2 - temp*u1
1206 w22 = w22 - temp*w12
1207 w21 = zero
1208*
1209* Compute SCALE
1210*
1211 scale = one
1212 IF( abs( w22 ).LT.safmin ) THEN
1213 scale = zero
1214 u2 = one
1215 u1 = -w12 / w11
1216 GO TO 250
1217 END IF
1218 IF( abs( w22 ).LT.abs( u2 ) )
1219 $ scale = abs( w22 / u2 )
1220 IF( abs( w11 ).LT.abs( u1 ) )
1221 $ scale = min( scale, abs( w11 / u1 ) )
1222*
1223* Solve
1224*
1225 u2 = ( scale*u2 ) / w22
1226 u1 = ( scale*u1-w12*u2 ) / w11
1227*
1228 250 CONTINUE
1229 IF( ilpivt ) THEN
1230 temp = u2
1231 u2 = u1
1232 u1 = temp
1233 END IF
1234*
1235* Compute Householder Vector
1236*
1237 t1 = sqrt( scale**2+u1**2+u2**2 )
1238 tau = one + scale / t1
1239 vs = -one / ( scale+t1 )
1240 v( 1 ) = one
1241 v( 2 ) = vs*u1
1242 v( 3 ) = vs*u2
1243*
1244* Apply transformations from the right.
1245*
1246 t2 = tau*v(2)
1247 t3 = tau*v(3)
1248 DO 260 jr = ifrstm, min( j+3, ilast )
1249 temp = h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
1250 $ h( jr, j+2 )
1251 h( jr, j ) = h( jr, j ) - temp*tau
1252 h( jr, j+1 ) = h( jr, j+1 ) - temp*t2
1253 h( jr, j+2 ) = h( jr, j+2 ) - temp*t3
1254 260 CONTINUE
1255 DO 270 jr = ifrstm, j + 2
1256 temp = t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
1257 $ t( jr, j+2 )
1258 t( jr, j ) = t( jr, j ) - temp*tau
1259 t( jr, j+1 ) = t( jr, j+1 ) - temp*t2
1260 t( jr, j+2 ) = t( jr, j+2 ) - temp*t3
1261 270 CONTINUE
1262 IF( ilz ) THEN
1263 DO 280 jr = 1, n
1264 temp = z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
1265 $ z( jr, j+2 )
1266 z( jr, j ) = z( jr, j ) - temp*tau
1267 z( jr, j+1 ) = z( jr, j+1 ) - temp*t2
1268 z( jr, j+2 ) = z( jr, j+2 ) - temp*t3
1269 280 CONTINUE
1270 END IF
1271 t( j+1, j ) = zero
1272 t( j+2, j ) = zero
1273 290 CONTINUE
1274*
1275* Last elements: Use Givens rotations
1276*
1277* Rotations from the left
1278*
1279 j = ilast - 1
1280 temp = h( j, j-1 )
1281 CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
1282 h( j+1, j-1 ) = zero
1283*
1284 DO 300 jc = j, ilastm
1285 temp = c*h( j, jc ) + s*h( j+1, jc )
1286 h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
1287 h( j, jc ) = temp
1288 temp2 = c*t( j, jc ) + s*t( j+1, jc )
1289 t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
1290 t( j, jc ) = temp2
1291 300 CONTINUE
1292 IF( ilq ) THEN
1293 DO 310 jr = 1, n
1294 temp = c*q( jr, j ) + s*q( jr, j+1 )
1295 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
1296 q( jr, j ) = temp
1297 310 CONTINUE
1298 END IF
1299*
1300* Rotations from the right.
1301*
1302 temp = t( j+1, j+1 )
1303 CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
1304 t( j+1, j ) = zero
1305*
1306 DO 320 jr = ifrstm, ilast
1307 temp = c*h( jr, j+1 ) + s*h( jr, j )
1308 h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
1309 h( jr, j+1 ) = temp
1310 320 CONTINUE
1311 DO 330 jr = ifrstm, ilast - 1
1312 temp = c*t( jr, j+1 ) + s*t( jr, j )
1313 t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
1314 t( jr, j+1 ) = temp
1315 330 CONTINUE
1316 IF( ilz ) THEN
1317 DO 340 jr = 1, n
1318 temp = c*z( jr, j+1 ) + s*z( jr, j )
1319 z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
1320 z( jr, j+1 ) = temp
1321 340 CONTINUE
1322 END IF
1323*
1324* End of Double-Shift code
1325*
1326 END IF
1327*
1328 GO TO 350
1329*
1330* End of iteration loop
1331*
1332 350 CONTINUE
1333 360 CONTINUE
1334*
1335* Drop-through = non-convergence
1336*
1337 info = ilast
1338 GO TO 420
1339*
1340* Successful completion of all QZ steps
1341*
1342 380 CONTINUE
1343*
1344* Set Eigenvalues 1:ILO-1
1345*
1346 DO 410 j = 1, ilo - 1
1347 IF( t( j, j ).LT.zero ) THEN
1348 IF( ilschr ) THEN
1349 DO 390 jr = 1, j
1350 h( jr, j ) = -h( jr, j )
1351 t( jr, j ) = -t( jr, j )
1352 390 CONTINUE
1353 ELSE
1354 h( j, j ) = -h( j, j )
1355 t( j, j ) = -t( j, j )
1356 END IF
1357 IF( ilz ) THEN
1358 DO 400 jr = 1, n
1359 z( jr, j ) = -z( jr, j )
1360 400 CONTINUE
1361 END IF
1362 END IF
1363 alphar( j ) = h( j, j )
1364 alphai( j ) = zero
1365 beta( j ) = t( j, j )
1366 410 CONTINUE
1367*
1368* Normal Termination
1369*
1370 info = 0
1371*
1372* Exit (other than argument error) -- return optimal workspace size
1373*
1374 420 CONTINUE
1375 work( 1 ) = dble( n )
1376 RETURN
1377*
1378* End of DHGEQZ
1379*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dlag2
subroutine dlag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition dlag2.f:154
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanhs
double precision function dlanhs(norm, n, a, lda, work)
DLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlanhs.f:106
dlapy2
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:61
dlapy3
double precision function dlapy3(x, y, z)
DLAPY3 returns sqrt(x2+y2+z2).
Definition dlapy3.f:66
dlarfg
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:104
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:108
dlasv2
subroutine dlasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition dlasv2.f:134
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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: