LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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, "An Algorithm for Generalized Matrix
      Eigenvalue Problems", 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 301 of file dhgeqz.f.

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