LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine chgeqz ( character  JOB,
character  COMPQ,
character  COMPZ,
integer  N,
integer  ILO,
integer  IHI,
complex, dimension( ldh, * )  H,
integer  LDH,
complex, dimension( ldt, * )  T,
integer  LDT,
complex, dimension( * )  ALPHA,
complex, dimension( * )  BETA,
complex, dimension( ldq, * )  Q,
integer  LDQ,
complex, dimension( ldz, * )  Z,
integer  LDZ,
complex, dimension( * )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CHGEQZ

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

Purpose:
 CHGEQZ computes the eigenvalues of a complex matrix pair (H,T),
 where H is an upper Hessenberg matrix and T is upper triangular,
 using the single-shift QZ method.
 Matrix pairs of this type are produced by the reduction to
 generalized upper Hessenberg form of a complex matrix pair (A,B):
 
    A = Q1*H*Z1**H,  B = Q1*T*Z1**H,
 
 as computed by CGGHRD.
 
 If JOB='S', then the Hessenberg-triangular pair (H,T) is
 also reduced to generalized Schur form,
 
    H = Q*S*Z**H,  T = Q*P*Z**H,
 
 where Q and Z are unitary matrices and S and P are upper triangular.
 
 Optionally, the unitary matrix Q from the generalized Schur
 factorization may be postmultiplied into an input matrix Q1, and the
 unitary matrix Z may be postmultiplied into an input matrix Z1.
 If Q1 and Z1 are the unitary matrices from CGGHRD that reduced
 the matrix pair (A,B) to generalized Hessenberg form, then the output
 matrices Q1*Q and Z1*Z are the unitary factors from the generalized
 Schur factorization of (A,B):
 
    A = (Q1*Q)*S*(Z1*Z)**H,  B = (Q1*Q)*P*(Z1*Z)**H.
 
 To avoid overflow, eigenvalues of the matrix pair (H,T)
 (equivalently, of (A,B)) are computed as a pair of complex values
 (alpha,beta).  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.
 The values of alpha and beta for the i-th eigenvalue 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': Computer 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 a unitary 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': Q is initialized to the unit matrix and the matrix Z
                 of right Schur vectors of (H,T) is returned;
          = 'V': Z must contain a unitary 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 COMPLEX array, dimension (LDH, N)
          On entry, the N-by-N upper Hessenberg matrix H.
          On exit, if JOB = 'S', H contains the upper triangular
          matrix S from the generalized Schur factorization.
          If JOB = 'E', the diagonal of H matches that 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 COMPLEX 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.
          If JOB = 'E', the diagonal of T matches that 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]ALPHA
          ALPHA is COMPLEX array, dimension (N)
          The complex scalars alpha that define the eigenvalues of
          GNEP.  ALPHA(i) = S(i,i) in the generalized Schur
          factorization.
[out]BETA
          BETA is COMPLEX array, dimension (N)
          The real non-negative scalars beta that define the
          eigenvalues of GNEP.  BETA(i) = P(i,i) in the generalized
          Schur factorization.

          Together, the quantities alpha = ALPHA(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 COMPLEX array, dimension (LDQ, N)
          On entry, if COMPQ = 'V', the unitary matrix Q1 used in the
          reduction of (A,B) to generalized Hessenberg form.
          On exit, if COMPQ = 'I', the unitary matrix of left Schur
          vectors of (H,T), and if COMPQ = 'V', the unitary 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 COMPLEX array, dimension (LDZ, N)
          On entry, if COMPZ = 'V', the unitary matrix Z1 used in the
          reduction of (A,B) to generalized Hessenberg form.
          On exit, if COMPZ = 'I', the unitary matrix of right Schur
          vectors of (H,T), and if COMPZ = 'V', the unitary 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 COMPLEX 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]RWORK
          RWORK is REAL array, dimension (N)
[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 ALPHA(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 ALPHA(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.
Date
April 2012
Further Details:
  We assume that complex ABS works as long as its value is less than
  overflow.

Definition at line 286 of file chgeqz.f.

286 *
287 * -- LAPACK computational routine (version 3.6.1) --
288 * -- LAPACK is a software package provided by Univ. of Tennessee, --
289 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
290 * April 2012
291 *
292 * .. Scalar Arguments ..
293  CHARACTER compq, compz, job
294  INTEGER ihi, ilo, info, ldh, ldq, ldt, ldz, lwork, n
295 * ..
296 * .. Array Arguments ..
297  REAL rwork( * )
298  COMPLEX alpha( * ), beta( * ), h( ldh, * ),
299  $ q( ldq, * ), t( ldt, * ), work( * ),
300  $ z( ldz, * )
301 * ..
302 *
303 * =====================================================================
304 *
305 * .. Parameters ..
306  COMPLEX czero, cone
307  parameter ( czero = ( 0.0e+0, 0.0e+0 ),
308  $ cone = ( 1.0e+0, 0.0e+0 ) )
309  REAL zero, one
310  parameter ( zero = 0.0e+0, one = 1.0e+0 )
311  REAL half
312  parameter ( half = 0.5e+0 )
313 * ..
314 * .. Local Scalars ..
315  LOGICAL ilazr2, ilazro, ilq, ilschr, ilz, lquery
316  INTEGER icompq, icompz, ifirst, ifrstm, iiter, ilast,
317  $ ilastm, in, ischur, istart, j, jc, jch, jiter,
318  $ jr, maxit
319  REAL absb, anorm, ascale, atol, bnorm, bscale, btol,
320  $ c, safmin, temp, temp2, tempr, ulp
321  COMPLEX abi22, ad11, ad12, ad21, ad22, ctemp, ctemp2,
322  $ ctemp3, eshift, rtdisc, s, shift, signbc, t1,
323  $ u12, x
324 * ..
325 * .. External Functions ..
326  LOGICAL lsame
327  REAL clanhs, slamch
328  EXTERNAL lsame, clanhs, slamch
329 * ..
330 * .. External Subroutines ..
331  EXTERNAL clartg, claset, crot, cscal, xerbla
332 * ..
333 * .. Intrinsic Functions ..
334  INTRINSIC abs, aimag, cmplx, conjg, max, min, REAL, sqrt
335 * ..
336 * .. Statement Functions ..
337  REAL abs1
338 * ..
339 * .. Statement Function definitions ..
340  abs1( x ) = abs( REAL( X ) ) + abs( aimag( x ) )
341 * ..
342 * .. Executable Statements ..
343 *
344 * Decode JOB, COMPQ, COMPZ
345 *
346  IF( lsame( job, 'E' ) ) THEN
347  ilschr = .false.
348  ischur = 1
349  ELSE IF( lsame( job, 'S' ) ) THEN
350  ilschr = .true.
351  ischur = 2
352  ELSE
353  ischur = 0
354  END IF
355 *
356  IF( lsame( compq, 'N' ) ) THEN
357  ilq = .false.
358  icompq = 1
359  ELSE IF( lsame( compq, 'V' ) ) THEN
360  ilq = .true.
361  icompq = 2
362  ELSE IF( lsame( compq, 'I' ) ) THEN
363  ilq = .true.
364  icompq = 3
365  ELSE
366  icompq = 0
367  END IF
368 *
369  IF( lsame( compz, 'N' ) ) THEN
370  ilz = .false.
371  icompz = 1
372  ELSE IF( lsame( compz, 'V' ) ) THEN
373  ilz = .true.
374  icompz = 2
375  ELSE IF( lsame( compz, 'I' ) ) THEN
376  ilz = .true.
377  icompz = 3
378  ELSE
379  icompz = 0
380  END IF
381 *
382 * Check Argument Values
383 *
384  info = 0
385  work( 1 ) = max( 1, n )
386  lquery = ( lwork.EQ.-1 )
387  IF( ischur.EQ.0 ) THEN
388  info = -1
389  ELSE IF( icompq.EQ.0 ) THEN
390  info = -2
391  ELSE IF( icompz.EQ.0 ) THEN
392  info = -3
393  ELSE IF( n.LT.0 ) THEN
394  info = -4
395  ELSE IF( ilo.LT.1 ) THEN
396  info = -5
397  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
398  info = -6
399  ELSE IF( ldh.LT.n ) THEN
400  info = -8
401  ELSE IF( ldt.LT.n ) THEN
402  info = -10
403  ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
404  info = -14
405  ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
406  info = -16
407  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
408  info = -18
409  END IF
410  IF( info.NE.0 ) THEN
411  CALL xerbla( 'CHGEQZ', -info )
412  RETURN
413  ELSE IF( lquery ) THEN
414  RETURN
415  END IF
416 *
417 * Quick return if possible
418 *
419 * WORK( 1 ) = CMPLX( 1 )
420  IF( n.LE.0 ) THEN
421  work( 1 ) = cmplx( 1 )
422  RETURN
423  END IF
424 *
425 * Initialize Q and Z
426 *
427  IF( icompq.EQ.3 )
428  $ CALL claset( 'Full', n, n, czero, cone, q, ldq )
429  IF( icompz.EQ.3 )
430  $ CALL claset( 'Full', n, n, czero, cone, z, ldz )
431 *
432 * Machine Constants
433 *
434  in = ihi + 1 - ilo
435  safmin = slamch( 'S' )
436  ulp = slamch( 'E' )*slamch( 'B' )
437  anorm = clanhs( 'F', in, h( ilo, ilo ), ldh, rwork )
438  bnorm = clanhs( 'F', in, t( ilo, ilo ), ldt, rwork )
439  atol = max( safmin, ulp*anorm )
440  btol = max( safmin, ulp*bnorm )
441  ascale = one / max( safmin, anorm )
442  bscale = one / max( safmin, bnorm )
443 *
444 *
445 * Set Eigenvalues IHI+1:N
446 *
447  DO 10 j = ihi + 1, n
448  absb = abs( t( j, j ) )
449  IF( absb.GT.safmin ) THEN
450  signbc = conjg( t( j, j ) / absb )
451  t( j, j ) = absb
452  IF( ilschr ) THEN
453  CALL cscal( j-1, signbc, t( 1, j ), 1 )
454  CALL cscal( j, signbc, h( 1, j ), 1 )
455  ELSE
456  CALL cscal( 1, signbc, h( j, j ), 1 )
457  END IF
458  IF( ilz )
459  $ CALL cscal( n, signbc, z( 1, j ), 1 )
460  ELSE
461  t( j, j ) = czero
462  END IF
463  alpha( j ) = h( j, j )
464  beta( j ) = t( j, j )
465  10 CONTINUE
466 *
467 * If IHI < ILO, skip QZ steps
468 *
469  IF( ihi.LT.ilo )
470  $ GO TO 190
471 *
472 * MAIN QZ ITERATION LOOP
473 *
474 * Initialize dynamic indices
475 *
476 * Eigenvalues ILAST+1:N have been found.
477 * Column operations modify rows IFRSTM:whatever
478 * Row operations modify columns whatever:ILASTM
479 *
480 * If only eigenvalues are being computed, then
481 * IFRSTM is the row of the last splitting row above row ILAST;
482 * this is always at least ILO.
483 * IITER counts iterations since the last eigenvalue was found,
484 * to tell when to use an extraordinary shift.
485 * MAXIT is the maximum number of QZ sweeps allowed.
486 *
487  ilast = ihi
488  IF( ilschr ) THEN
489  ifrstm = 1
490  ilastm = n
491  ELSE
492  ifrstm = ilo
493  ilastm = ihi
494  END IF
495  iiter = 0
496  eshift = czero
497  maxit = 30*( ihi-ilo+1 )
498 *
499  DO 170 jiter = 1, maxit
500 *
501 * Check for too many iterations.
502 *
503  IF( jiter.GT.maxit )
504  $ GO TO 180
505 *
506 * Split the matrix if possible.
507 *
508 * Two tests:
509 * 1: H(j,j-1)=0 or j=ILO
510 * 2: T(j,j)=0
511 *
512 * Special case: j=ILAST
513 *
514  IF( ilast.EQ.ilo ) THEN
515  GO TO 60
516  ELSE
517  IF( abs1( h( ilast, ilast-1 ) ).LE.atol ) THEN
518  h( ilast, ilast-1 ) = czero
519  GO TO 60
520  END IF
521  END IF
522 *
523  IF( abs( t( ilast, ilast ) ).LE.btol ) THEN
524  t( ilast, ilast ) = czero
525  GO TO 50
526  END IF
527 *
528 * General case: j<ILAST
529 *
530  DO 40 j = ilast - 1, ilo, -1
531 *
532 * Test 1: for H(j,j-1)=0 or j=ILO
533 *
534  IF( j.EQ.ilo ) THEN
535  ilazro = .true.
536  ELSE
537  IF( abs1( h( j, j-1 ) ).LE.atol ) THEN
538  h( j, j-1 ) = czero
539  ilazro = .true.
540  ELSE
541  ilazro = .false.
542  END IF
543  END IF
544 *
545 * Test 2: for T(j,j)=0
546 *
547  IF( abs( t( j, j ) ).LT.btol ) THEN
548  t( j, j ) = czero
549 *
550 * Test 1a: Check for 2 consecutive small subdiagonals in A
551 *
552  ilazr2 = .false.
553  IF( .NOT.ilazro ) THEN
554  IF( abs1( h( j, j-1 ) )*( ascale*abs1( h( j+1,
555  $ j ) ) ).LE.abs1( h( j, j ) )*( ascale*atol ) )
556  $ ilazr2 = .true.
557  END IF
558 *
559 * If both tests pass (1 & 2), i.e., the leading diagonal
560 * element of B in the block is zero, split a 1x1 block off
561 * at the top. (I.e., at the J-th row/column) The leading
562 * diagonal element of the remainder can also be zero, so
563 * this may have to be done repeatedly.
564 *
565  IF( ilazro .OR. ilazr2 ) THEN
566  DO 20 jch = j, ilast - 1
567  ctemp = h( jch, jch )
568  CALL clartg( ctemp, h( jch+1, jch ), c, s,
569  $ h( jch, jch ) )
570  h( jch+1, jch ) = czero
571  CALL crot( ilastm-jch, h( jch, jch+1 ), ldh,
572  $ h( jch+1, jch+1 ), ldh, c, s )
573  CALL crot( ilastm-jch, t( jch, jch+1 ), ldt,
574  $ t( jch+1, jch+1 ), ldt, c, s )
575  IF( ilq )
576  $ CALL crot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
577  $ c, conjg( s ) )
578  IF( ilazr2 )
579  $ h( jch, jch-1 ) = h( jch, jch-1 )*c
580  ilazr2 = .false.
581  IF( abs1( t( jch+1, jch+1 ) ).GE.btol ) THEN
582  IF( jch+1.GE.ilast ) THEN
583  GO TO 60
584  ELSE
585  ifirst = jch + 1
586  GO TO 70
587  END IF
588  END IF
589  t( jch+1, jch+1 ) = czero
590  20 CONTINUE
591  GO TO 50
592  ELSE
593 *
594 * Only test 2 passed -- chase the zero to T(ILAST,ILAST)
595 * Then process as in the case T(ILAST,ILAST)=0
596 *
597  DO 30 jch = j, ilast - 1
598  ctemp = t( jch, jch+1 )
599  CALL clartg( ctemp, t( jch+1, jch+1 ), c, s,
600  $ t( jch, jch+1 ) )
601  t( jch+1, jch+1 ) = czero
602  IF( jch.LT.ilastm-1 )
603  $ CALL crot( ilastm-jch-1, t( jch, jch+2 ), ldt,
604  $ t( jch+1, jch+2 ), ldt, c, s )
605  CALL crot( ilastm-jch+2, h( jch, jch-1 ), ldh,
606  $ h( jch+1, jch-1 ), ldh, c, s )
607  IF( ilq )
608  $ CALL crot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
609  $ c, conjg( s ) )
610  ctemp = h( jch+1, jch )
611  CALL clartg( ctemp, h( jch+1, jch-1 ), c, s,
612  $ h( jch+1, jch ) )
613  h( jch+1, jch-1 ) = czero
614  CALL crot( jch+1-ifrstm, h( ifrstm, jch ), 1,
615  $ h( ifrstm, jch-1 ), 1, c, s )
616  CALL crot( jch-ifrstm, t( ifrstm, jch ), 1,
617  $ t( ifrstm, jch-1 ), 1, c, s )
618  IF( ilz )
619  $ CALL crot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
620  $ c, s )
621  30 CONTINUE
622  GO TO 50
623  END IF
624  ELSE IF( ilazro ) THEN
625 *
626 * Only test 1 passed -- work on J:ILAST
627 *
628  ifirst = j
629  GO TO 70
630  END IF
631 *
632 * Neither test passed -- try next J
633 *
634  40 CONTINUE
635 *
636 * (Drop-through is "impossible")
637 *
638  info = 2*n + 1
639  GO TO 210
640 *
641 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
642 * 1x1 block.
643 *
644  50 CONTINUE
645  ctemp = h( ilast, ilast )
646  CALL clartg( ctemp, h( ilast, ilast-1 ), c, s,
647  $ h( ilast, ilast ) )
648  h( ilast, ilast-1 ) = czero
649  CALL crot( ilast-ifrstm, h( ifrstm, ilast ), 1,
650  $ h( ifrstm, ilast-1 ), 1, c, s )
651  CALL crot( ilast-ifrstm, t( ifrstm, ilast ), 1,
652  $ t( ifrstm, ilast-1 ), 1, c, s )
653  IF( ilz )
654  $ CALL crot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
655 *
656 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
657 *
658  60 CONTINUE
659  absb = abs( t( ilast, ilast ) )
660  IF( absb.GT.safmin ) THEN
661  signbc = conjg( t( ilast, ilast ) / absb )
662  t( ilast, ilast ) = absb
663  IF( ilschr ) THEN
664  CALL cscal( ilast-ifrstm, signbc, t( ifrstm, ilast ), 1 )
665  CALL cscal( ilast+1-ifrstm, signbc, h( ifrstm, ilast ),
666  $ 1 )
667  ELSE
668  CALL cscal( 1, signbc, h( ilast, ilast ), 1 )
669  END IF
670  IF( ilz )
671  $ CALL cscal( n, signbc, z( 1, ilast ), 1 )
672  ELSE
673  t( ilast, ilast ) = czero
674  END IF
675  alpha( ilast ) = h( ilast, ilast )
676  beta( ilast ) = t( ilast, ilast )
677 *
678 * Go to next block -- exit if finished.
679 *
680  ilast = ilast - 1
681  IF( ilast.LT.ilo )
682  $ GO TO 190
683 *
684 * Reset counters
685 *
686  iiter = 0
687  eshift = czero
688  IF( .NOT.ilschr ) THEN
689  ilastm = ilast
690  IF( ifrstm.GT.ilast )
691  $ ifrstm = ilo
692  END IF
693  GO TO 160
694 *
695 * QZ step
696 *
697 * This iteration only involves rows/columns IFIRST:ILAST. We
698 * assume IFIRST < ILAST, and that the diagonal of B is non-zero.
699 *
700  70 CONTINUE
701  iiter = iiter + 1
702  IF( .NOT.ilschr ) THEN
703  ifrstm = ifirst
704  END IF
705 *
706 * Compute the Shift.
707 *
708 * At this point, IFIRST < ILAST, and the diagonal elements of
709 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
710 * magnitude)
711 *
712  IF( ( iiter / 10 )*10.NE.iiter ) THEN
713 *
714 * The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
715 * the bottom-right 2x2 block of A inv(B) which is nearest to
716 * the bottom-right element.
717 *
718 * We factor B as U*D, where U has unit diagonals, and
719 * compute (A*inv(D))*inv(U).
720 *
721  u12 = ( bscale*t( ilast-1, ilast ) ) /
722  $ ( bscale*t( ilast, ilast ) )
723  ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
724  $ ( bscale*t( ilast-1, ilast-1 ) )
725  ad21 = ( ascale*h( ilast, ilast-1 ) ) /
726  $ ( bscale*t( ilast-1, ilast-1 ) )
727  ad12 = ( ascale*h( ilast-1, ilast ) ) /
728  $ ( bscale*t( ilast, ilast ) )
729  ad22 = ( ascale*h( ilast, ilast ) ) /
730  $ ( bscale*t( ilast, ilast ) )
731  abi22 = ad22 - u12*ad21
732 *
733  t1 = half*( ad11+abi22 )
734  rtdisc = sqrt( t1**2+ad12*ad21-ad11*ad22 )
735  temp = REAL( t1-abi22 )*REAL( RTDISC ) +
736  $ aimag( t1-abi22 )*aimag( rtdisc )
737  IF( temp.LE.zero ) THEN
738  shift = t1 + rtdisc
739  ELSE
740  shift = t1 - rtdisc
741  END IF
742  ELSE
743 *
744 * Exceptional shift. Chosen for no particularly good reason.
745 *
746  eshift = eshift + (ascale*h(ilast,ilast-1))/
747  $ (bscale*t(ilast-1,ilast-1))
748  shift = eshift
749  END IF
750 *
751 * Now check for two consecutive small subdiagonals.
752 *
753  DO 80 j = ilast - 1, ifirst + 1, -1
754  istart = j
755  ctemp = ascale*h( j, j ) - shift*( bscale*t( j, j ) )
756  temp = abs1( ctemp )
757  temp2 = ascale*abs1( h( j+1, j ) )
758  tempr = max( temp, temp2 )
759  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
760  temp = temp / tempr
761  temp2 = temp2 / tempr
762  END IF
763  IF( abs1( h( j, j-1 ) )*temp2.LE.temp*atol )
764  $ GO TO 90
765  80 CONTINUE
766 *
767  istart = ifirst
768  ctemp = ascale*h( ifirst, ifirst ) -
769  $ shift*( bscale*t( ifirst, ifirst ) )
770  90 CONTINUE
771 *
772 * Do an implicit-shift QZ sweep.
773 *
774 * Initial Q
775 *
776  ctemp2 = ascale*h( istart+1, istart )
777  CALL clartg( ctemp, ctemp2, c, s, ctemp3 )
778 *
779 * Sweep
780 *
781  DO 150 j = istart, ilast - 1
782  IF( j.GT.istart ) THEN
783  ctemp = h( j, j-1 )
784  CALL clartg( ctemp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
785  h( j+1, j-1 ) = czero
786  END IF
787 *
788  DO 100 jc = j, ilastm
789  ctemp = c*h( j, jc ) + s*h( j+1, jc )
790  h( j+1, jc ) = -conjg( s )*h( j, jc ) + c*h( j+1, jc )
791  h( j, jc ) = ctemp
792  ctemp2 = c*t( j, jc ) + s*t( j+1, jc )
793  t( j+1, jc ) = -conjg( s )*t( j, jc ) + c*t( j+1, jc )
794  t( j, jc ) = ctemp2
795  100 CONTINUE
796  IF( ilq ) THEN
797  DO 110 jr = 1, n
798  ctemp = c*q( jr, j ) + conjg( s )*q( jr, j+1 )
799  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
800  q( jr, j ) = ctemp
801  110 CONTINUE
802  END IF
803 *
804  ctemp = t( j+1, j+1 )
805  CALL clartg( ctemp, t( j+1, j ), c, s, t( j+1, j+1 ) )
806  t( j+1, j ) = czero
807 *
808  DO 120 jr = ifrstm, min( j+2, ilast )
809  ctemp = c*h( jr, j+1 ) + s*h( jr, j )
810  h( jr, j ) = -conjg( s )*h( jr, j+1 ) + c*h( jr, j )
811  h( jr, j+1 ) = ctemp
812  120 CONTINUE
813  DO 130 jr = ifrstm, j
814  ctemp = c*t( jr, j+1 ) + s*t( jr, j )
815  t( jr, j ) = -conjg( s )*t( jr, j+1 ) + c*t( jr, j )
816  t( jr, j+1 ) = ctemp
817  130 CONTINUE
818  IF( ilz ) THEN
819  DO 140 jr = 1, n
820  ctemp = c*z( jr, j+1 ) + s*z( jr, j )
821  z( jr, j ) = -conjg( s )*z( jr, j+1 ) + c*z( jr, j )
822  z( jr, j+1 ) = ctemp
823  140 CONTINUE
824  END IF
825  150 CONTINUE
826 *
827  160 CONTINUE
828 *
829  170 CONTINUE
830 *
831 * Drop-through = non-convergence
832 *
833  180 CONTINUE
834  info = ilast
835  GO TO 210
836 *
837 * Successful completion of all QZ steps
838 *
839  190 CONTINUE
840 *
841 * Set Eigenvalues 1:ILO-1
842 *
843  DO 200 j = 1, ilo - 1
844  absb = abs( t( j, j ) )
845  IF( absb.GT.safmin ) THEN
846  signbc = conjg( t( j, j ) / absb )
847  t( j, j ) = absb
848  IF( ilschr ) THEN
849  CALL cscal( j-1, signbc, t( 1, j ), 1 )
850  CALL cscal( j, signbc, h( 1, j ), 1 )
851  ELSE
852  CALL cscal( 1, signbc, h( j, j ), 1 )
853  END IF
854  IF( ilz )
855  $ CALL cscal( n, signbc, z( 1, j ), 1 )
856  ELSE
857  t( j, j ) = czero
858  END IF
859  alpha( j ) = h( j, j )
860  beta( j ) = t( j, j )
861  200 CONTINUE
862 *
863 * Normal Termination
864 *
865  info = 0
866 *
867 * Exit (other than argument error) -- return optimal workspace size
868 *
869  210 CONTINUE
870  work( 1 ) = cmplx( n )
871  RETURN
872 *
873 * End of CHGEQZ
874 *
subroutine clartg(F, G, CS, SN, R)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f:105
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
real function clanhs(NORM, N, A, LDA, WORK)
CLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clanhs.f:111
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
Definition: crot.f:105

Here is the call graph for this function:

Here is the caller graph for this function: