/* zhgeqz.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "blaswrap.h"
/* Table of constant values */
static doublecomplex c_b1 = {0.,0.};
static doublecomplex c_b2 = {1.,0.};
static integer c__1 = 1;
static integer c__2 = 2;
/* Subroutine */ int zhgeqz_(char *job, char *compq, char *compz, integer *n,
integer *ilo, integer *ihi, doublecomplex *h__, integer *ldh,
doublecomplex *t, integer *ldt, doublecomplex *alpha, doublecomplex *
beta, doublecomplex *q, integer *ldq, doublecomplex *z__, integer *
ldz, doublecomplex *work, integer *lwork, doublereal *rwork, integer *
info)
{
/* System generated locals */
integer h_dim1, h_offset, q_dim1, q_offset, t_dim1, t_offset, z_dim1,
z_offset, i__1, i__2, i__3, i__4, i__5, i__6;
doublereal d__1, d__2, d__3, d__4, d__5, d__6;
doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;
/* Builtin functions */
double z_abs(doublecomplex *);
void d_cnjg(doublecomplex *, doublecomplex *);
double d_imag(doublecomplex *);
void z_div(doublecomplex *, doublecomplex *, doublecomplex *), pow_zi(
doublecomplex *, doublecomplex *, integer *), z_sqrt(
doublecomplex *, doublecomplex *);
/* Local variables */
doublereal c__;
integer j;
doublecomplex s, t1;
integer jc, in;
doublecomplex u12;
integer jr;
doublecomplex ad11, ad12, ad21, ad22;
integer jch;
logical ilq, ilz;
doublereal ulp;
doublecomplex abi22;
doublereal absb, atol, btol, temp;
extern /* Subroutine */ int zrot_(integer *, doublecomplex *, integer *,
doublecomplex *, integer *, doublereal *, doublecomplex *);
doublereal temp2;
extern logical lsame_(char *, char *);
doublecomplex ctemp;
integer iiter, ilast, jiter;
doublereal anorm, bnorm;
integer maxit;
doublecomplex shift;
extern /* Subroutine */ int zscal_(integer *, doublecomplex *,
doublecomplex *, integer *);
doublereal tempr;
doublecomplex ctemp2, ctemp3;
logical ilazr2;
doublereal ascale, bscale;
extern doublereal dlamch_(char *);
doublecomplex signbc;
doublereal safmin;
extern /* Subroutine */ int xerbla_(char *, integer *);
doublecomplex eshift;
logical ilschr;
integer icompq, ilastm;
doublecomplex rtdisc;
integer ischur;
extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *,
doublereal *);
logical ilazro;
integer icompz, ifirst;
extern /* Subroutine */ int zlartg_(doublecomplex *, doublecomplex *,
doublereal *, doublecomplex *, doublecomplex *);
integer ifrstm;
extern /* Subroutine */ int zlaset_(char *, integer *, integer *,
doublecomplex *, doublecomplex *, doublecomplex *, integer *);
integer istart;
logical lquery;
/* -- LAPACK routine (version 3.2) -- */
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZHGEQZ 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 ZGGHRD. */
/* 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 ZGGHRD 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. */
/* Arguments */
/* ========= */
/* JOB (input) CHARACTER*1 */
/* = 'E': Compute eigenvalues only; */
/* = 'S': Computer eigenvalues and the Schur form. */
/* COMPQ (input) 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. */
/* COMPZ (input) 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. */
/* N (input) INTEGER */
/* The order of the matrices H, T, Q, and Z. N >= 0. */
/* ILO (input) INTEGER */
/* IHI (input) 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. */
/* H (input/output) COMPLEX*16 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. */
/* LDH (input) INTEGER */
/* The leading dimension of the array H. LDH >= max( 1, N ). */
/* T (input/output) COMPLEX*16 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. */
/* LDT (input) INTEGER */
/* The leading dimension of the array T. LDT >= max( 1, N ). */
/* ALPHA (output) COMPLEX*16 array, dimension (N) */
/* The complex scalars alpha that define the eigenvalues of */
/* GNEP. ALPHA(i) = S(i,i) in the generalized Schur */
/* factorization. */
/* BETA (output) COMPLEX*16 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. */
/* Q (input/output) COMPLEX*16 array, dimension (LDQ, N) */
/* On entry, if COMPZ = 'V', the unitary matrix Q1 used in the */
/* reduction of (A,B) to generalized Hessenberg form. */
/* On exit, if COMPZ = 'I', the unitary matrix of left Schur */
/* vectors of (H,T), and if COMPZ = 'V', the unitary matrix of */
/* left Schur vectors of (A,B). */
/* Not referenced if COMPZ = 'N'. */
/* LDQ (input) INTEGER */
/* The leading dimension of the array Q. LDQ >= 1. */
/* If COMPQ='V' or 'I', then LDQ >= N. */
/* Z (input/output) COMPLEX*16 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'. */
/* LDZ (input) INTEGER */
/* The leading dimension of the array Z. LDZ >= 1. */
/* If COMPZ='V' or 'I', then LDZ >= N. */
/* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) */
/* On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. */
/* LWORK (input) 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. */
/* RWORK (workspace) DOUBLE PRECISION array, dimension (N) */
/* INFO (output) 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. */
/* Further Details */
/* =============== */
/* We assume that complex ABS works as long as its value is less than */
/* overflow. */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Statement Functions .. */
/* .. */
/* .. Statement Function definitions .. */
/* .. */
/* .. Executable Statements .. */
/* Decode JOB, COMPQ, COMPZ */
/* Parameter adjustments */
h_dim1 = *ldh;
h_offset = 1 + h_dim1;
h__ -= h_offset;
t_dim1 = *ldt;
t_offset = 1 + t_dim1;
t -= t_offset;
--alpha;
--beta;
q_dim1 = *ldq;
q_offset = 1 + q_dim1;
q -= q_offset;
z_dim1 = *ldz;
z_offset = 1 + z_dim1;
z__ -= z_offset;
--work;
--rwork;
/* Function Body */
if (lsame_(job, "E")) {
ilschr = FALSE_;
ischur = 1;
} else if (lsame_(job, "S")) {
ilschr = TRUE_;
ischur = 2;
} else {
ischur = 0;
}
if (lsame_(compq, "N")) {
ilq = FALSE_;
icompq = 1;
} else if (lsame_(compq, "V")) {
ilq = TRUE_;
icompq = 2;
} else if (lsame_(compq, "I")) {
ilq = TRUE_;
icompq = 3;
} else {
icompq = 0;
}
if (lsame_(compz, "N")) {
ilz = FALSE_;
icompz = 1;
} else if (lsame_(compz, "V")) {
ilz = TRUE_;
icompz = 2;
} else if (lsame_(compz, "I")) {
ilz = TRUE_;
icompz = 3;
} else {
icompz = 0;
}
/* Check Argument Values */
*info = 0;
i__1 = max(1,*n);
work[1].r = (doublereal) i__1, work[1].i = 0.;
lquery = *lwork == -1;
if (ischur == 0) {
*info = -1;
} else if (icompq == 0) {
*info = -2;
} else if (icompz == 0) {
*info = -3;
} else if (*n < 0) {
*info = -4;
} else if (*ilo < 1) {
*info = -5;
} else if (*ihi > *n || *ihi < *ilo - 1) {
*info = -6;
} else if (*ldh < *n) {
*info = -8;
} else if (*ldt < *n) {
*info = -10;
} else if (*ldq < 1 || ilq && *ldq < *n) {
*info = -14;
} else if (*ldz < 1 || ilz && *ldz < *n) {
*info = -16;
} else if (*lwork < max(1,*n) && ! lquery) {
*info = -18;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("ZHGEQZ", &i__1);
return 0;
} else if (lquery) {
return 0;
}
/* Quick return if possible */
/* WORK( 1 ) = CMPLX( 1 ) */
if (*n <= 0) {
work[1].r = 1., work[1].i = 0.;
return 0;
}
/* Initialize Q and Z */
if (icompq == 3) {
zlaset_("Full", n, n, &c_b1, &c_b2, &q[q_offset], ldq);
}
if (icompz == 3) {
zlaset_("Full", n, n, &c_b1, &c_b2, &z__[z_offset], ldz);
}
/* Machine Constants */
in = *ihi + 1 - *ilo;
safmin = dlamch_("S");
ulp = dlamch_("E") * dlamch_("B");
anorm = zlanhs_("F", &in, &h__[*ilo + *ilo * h_dim1], ldh, &rwork[1]);
bnorm = zlanhs_("F", &in, &t[*ilo + *ilo * t_dim1], ldt, &rwork[1]);
/* Computing MAX */
d__1 = safmin, d__2 = ulp * anorm;
atol = max(d__1,d__2);
/* Computing MAX */
d__1 = safmin, d__2 = ulp * bnorm;
btol = max(d__1,d__2);
ascale = 1. / max(safmin,anorm);
bscale = 1. / max(safmin,bnorm);
/* Set Eigenvalues IHI+1:N */
i__1 = *n;
for (j = *ihi + 1; j <= i__1; ++j) {
absb = z_abs(&t[j + j * t_dim1]);
if (absb > safmin) {
i__2 = j + j * t_dim1;
z__2.r = t[i__2].r / absb, z__2.i = t[i__2].i / absb;
d_cnjg(&z__1, &z__2);
signbc.r = z__1.r, signbc.i = z__1.i;
i__2 = j + j * t_dim1;
t[i__2].r = absb, t[i__2].i = 0.;
if (ilschr) {
i__2 = j - 1;
zscal_(&i__2, &signbc, &t[j * t_dim1 + 1], &c__1);
zscal_(&j, &signbc, &h__[j * h_dim1 + 1], &c__1);
} else {
i__2 = j + j * h_dim1;
i__3 = j + j * h_dim1;
z__1.r = h__[i__3].r * signbc.r - h__[i__3].i * signbc.i,
z__1.i = h__[i__3].r * signbc.i + h__[i__3].i *
signbc.r;
h__[i__2].r = z__1.r, h__[i__2].i = z__1.i;
}
if (ilz) {
zscal_(n, &signbc, &z__[j * z_dim1 + 1], &c__1);
}
} else {
i__2 = j + j * t_dim1;
t[i__2].r = 0., t[i__2].i = 0.;
}
i__2 = j;
i__3 = j + j * h_dim1;
alpha[i__2].r = h__[i__3].r, alpha[i__2].i = h__[i__3].i;
i__2 = j;
i__3 = j + j * t_dim1;
beta[i__2].r = t[i__3].r, beta[i__2].i = t[i__3].i;
/* L10: */
}
/* If IHI < ILO, skip QZ steps */
if (*ihi < *ilo) {
goto L190;
}
/* MAIN QZ ITERATION LOOP */
/* Initialize dynamic indices */
/* Eigenvalues ILAST+1:N have been found. */
/* Column operations modify rows IFRSTM:whatever */
/* Row operations modify columns whatever:ILASTM */
/* If only eigenvalues are being computed, then */
/* IFRSTM is the row of the last splitting row above row ILAST; */
/* this is always at least ILO. */
/* IITER counts iterations since the last eigenvalue was found, */
/* to tell when to use an extraordinary shift. */
/* MAXIT is the maximum number of QZ sweeps allowed. */
ilast = *ihi;
if (ilschr) {
ifrstm = 1;
ilastm = *n;
} else {
ifrstm = *ilo;
ilastm = *ihi;
}
iiter = 0;
eshift.r = 0., eshift.i = 0.;
maxit = (*ihi - *ilo + 1) * 30;
i__1 = maxit;
for (jiter = 1; jiter <= i__1; ++jiter) {
/* Check for too many iterations. */
if (jiter > maxit) {
goto L180;
}
/* Split the matrix if possible. */
/* Two tests: */
/* 1: H(j,j-1)=0 or j=ILO */
/* 2: T(j,j)=0 */
/* Special case: j=ILAST */
if (ilast == *ilo) {
goto L60;
} else {
i__2 = ilast + (ilast - 1) * h_dim1;
if ((d__1 = h__[i__2].r, abs(d__1)) + (d__2 = d_imag(&h__[ilast +
(ilast - 1) * h_dim1]), abs(d__2)) <= atol) {
i__2 = ilast + (ilast - 1) * h_dim1;
h__[i__2].r = 0., h__[i__2].i = 0.;
goto L60;
}
}
if (z_abs(&t[ilast + ilast * t_dim1]) <= btol) {
i__2 = ilast + ilast * t_dim1;
t[i__2].r = 0., t[i__2].i = 0.;
goto L50;
}
/* General case: j<ILAST */
i__2 = *ilo;
for (j = ilast - 1; j >= i__2; --j) {
/* Test 1: for H(j,j-1)=0 or j=ILO */
if (j == *ilo) {
ilazro = TRUE_;
} else {
i__3 = j + (j - 1) * h_dim1;
if ((d__1 = h__[i__3].r, abs(d__1)) + (d__2 = d_imag(&h__[j +
(j - 1) * h_dim1]), abs(d__2)) <= atol) {
i__3 = j + (j - 1) * h_dim1;
h__[i__3].r = 0., h__[i__3].i = 0.;
ilazro = TRUE_;
} else {
ilazro = FALSE_;
}
}
/* Test 2: for T(j,j)=0 */
if (z_abs(&t[j + j * t_dim1]) < btol) {
i__3 = j + j * t_dim1;
t[i__3].r = 0., t[i__3].i = 0.;
/* Test 1a: Check for 2 consecutive small subdiagonals in A */
ilazr2 = FALSE_;
if (! ilazro) {
i__3 = j + (j - 1) * h_dim1;
i__4 = j + 1 + j * h_dim1;
i__5 = j + j * h_dim1;
if (((d__1 = h__[i__3].r, abs(d__1)) + (d__2 = d_imag(&
h__[j + (j - 1) * h_dim1]), abs(d__2))) * (ascale
* ((d__3 = h__[i__4].r, abs(d__3)) + (d__4 =
d_imag(&h__[j + 1 + j * h_dim1]), abs(d__4)))) <=
((d__5 = h__[i__5].r, abs(d__5)) + (d__6 = d_imag(
&h__[j + j * h_dim1]), abs(d__6))) * (ascale *
atol)) {
ilazr2 = TRUE_;
}
}
/* If both tests pass (1 & 2), i.e., the leading diagonal */
/* element of B in the block is zero, split a 1x1 block off */
/* at the top. (I.e., at the J-th row/column) The leading */
/* diagonal element of the remainder can also be zero, so */
/* this may have to be done repeatedly. */
if (ilazro || ilazr2) {
i__3 = ilast - 1;
for (jch = j; jch <= i__3; ++jch) {
i__4 = jch + jch * h_dim1;
ctemp.r = h__[i__4].r, ctemp.i = h__[i__4].i;
zlartg_(&ctemp, &h__[jch + 1 + jch * h_dim1], &c__, &
s, &h__[jch + jch * h_dim1]);
i__4 = jch + 1 + jch * h_dim1;
h__[i__4].r = 0., h__[i__4].i = 0.;
i__4 = ilastm - jch;
zrot_(&i__4, &h__[jch + (jch + 1) * h_dim1], ldh, &
h__[jch + 1 + (jch + 1) * h_dim1], ldh, &c__,
&s);
i__4 = ilastm - jch;
zrot_(&i__4, &t[jch + (jch + 1) * t_dim1], ldt, &t[
jch + 1 + (jch + 1) * t_dim1], ldt, &c__, &s);
if (ilq) {
d_cnjg(&z__1, &s);
zrot_(n, &q[jch * q_dim1 + 1], &c__1, &q[(jch + 1)
* q_dim1 + 1], &c__1, &c__, &z__1);
}
if (ilazr2) {
i__4 = jch + (jch - 1) * h_dim1;
i__5 = jch + (jch - 1) * h_dim1;
z__1.r = c__ * h__[i__5].r, z__1.i = c__ * h__[
i__5].i;
h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
}
ilazr2 = FALSE_;
i__4 = jch + 1 + (jch + 1) * t_dim1;
if ((d__1 = t[i__4].r, abs(d__1)) + (d__2 = d_imag(&t[
jch + 1 + (jch + 1) * t_dim1]), abs(d__2)) >=
btol) {
if (jch + 1 >= ilast) {
goto L60;
} else {
ifirst = jch + 1;
goto L70;
}
}
i__4 = jch + 1 + (jch + 1) * t_dim1;
t[i__4].r = 0., t[i__4].i = 0.;
/* L20: */
}
goto L50;
} else {
/* Only test 2 passed -- chase the zero to T(ILAST,ILAST) */
/* Then process as in the case T(ILAST,ILAST)=0 */
i__3 = ilast - 1;
for (jch = j; jch <= i__3; ++jch) {
i__4 = jch + (jch + 1) * t_dim1;
ctemp.r = t[i__4].r, ctemp.i = t[i__4].i;
zlartg_(&ctemp, &t[jch + 1 + (jch + 1) * t_dim1], &
c__, &s, &t[jch + (jch + 1) * t_dim1]);
i__4 = jch + 1 + (jch + 1) * t_dim1;
t[i__4].r = 0., t[i__4].i = 0.;
if (jch < ilastm - 1) {
i__4 = ilastm - jch - 1;
zrot_(&i__4, &t[jch + (jch + 2) * t_dim1], ldt, &
t[jch + 1 + (jch + 2) * t_dim1], ldt, &
c__, &s);
}
i__4 = ilastm - jch + 2;
zrot_(&i__4, &h__[jch + (jch - 1) * h_dim1], ldh, &
h__[jch + 1 + (jch - 1) * h_dim1], ldh, &c__,
&s);
if (ilq) {
d_cnjg(&z__1, &s);
zrot_(n, &q[jch * q_dim1 + 1], &c__1, &q[(jch + 1)
* q_dim1 + 1], &c__1, &c__, &z__1);
}
i__4 = jch + 1 + jch * h_dim1;
ctemp.r = h__[i__4].r, ctemp.i = h__[i__4].i;
zlartg_(&ctemp, &h__[jch + 1 + (jch - 1) * h_dim1], &
c__, &s, &h__[jch + 1 + jch * h_dim1]);
i__4 = jch + 1 + (jch - 1) * h_dim1;
h__[i__4].r = 0., h__[i__4].i = 0.;
i__4 = jch + 1 - ifrstm;
zrot_(&i__4, &h__[ifrstm + jch * h_dim1], &c__1, &h__[
ifrstm + (jch - 1) * h_dim1], &c__1, &c__, &s)
;
i__4 = jch - ifrstm;
zrot_(&i__4, &t[ifrstm + jch * t_dim1], &c__1, &t[
ifrstm + (jch - 1) * t_dim1], &c__1, &c__, &s)
;
if (ilz) {
zrot_(n, &z__[jch * z_dim1 + 1], &c__1, &z__[(jch
- 1) * z_dim1 + 1], &c__1, &c__, &s);
}
/* L30: */
}
goto L50;
}
} else if (ilazro) {
/* Only test 1 passed -- work on J:ILAST */
ifirst = j;
goto L70;
}
/* Neither test passed -- try next J */
/* L40: */
}
/* (Drop-through is "impossible") */
*info = (*n << 1) + 1;
goto L210;
/* T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a */
/* 1x1 block. */
L50:
i__2 = ilast + ilast * h_dim1;
ctemp.r = h__[i__2].r, ctemp.i = h__[i__2].i;
zlartg_(&ctemp, &h__[ilast + (ilast - 1) * h_dim1], &c__, &s, &h__[
ilast + ilast * h_dim1]);
i__2 = ilast + (ilast - 1) * h_dim1;
h__[i__2].r = 0., h__[i__2].i = 0.;
i__2 = ilast - ifrstm;
zrot_(&i__2, &h__[ifrstm + ilast * h_dim1], &c__1, &h__[ifrstm + (
ilast - 1) * h_dim1], &c__1, &c__, &s);
i__2 = ilast - ifrstm;
zrot_(&i__2, &t[ifrstm + ilast * t_dim1], &c__1, &t[ifrstm + (ilast -
1) * t_dim1], &c__1, &c__, &s);
if (ilz) {
zrot_(n, &z__[ilast * z_dim1 + 1], &c__1, &z__[(ilast - 1) *
z_dim1 + 1], &c__1, &c__, &s);
}
/* H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA */
L60:
absb = z_abs(&t[ilast + ilast * t_dim1]);
if (absb > safmin) {
i__2 = ilast + ilast * t_dim1;
z__2.r = t[i__2].r / absb, z__2.i = t[i__2].i / absb;
d_cnjg(&z__1, &z__2);
signbc.r = z__1.r, signbc.i = z__1.i;
i__2 = ilast + ilast * t_dim1;
t[i__2].r = absb, t[i__2].i = 0.;
if (ilschr) {
i__2 = ilast - ifrstm;
zscal_(&i__2, &signbc, &t[ifrstm + ilast * t_dim1], &c__1);
i__2 = ilast + 1 - ifrstm;
zscal_(&i__2, &signbc, &h__[ifrstm + ilast * h_dim1], &c__1);
} else {
i__2 = ilast + ilast * h_dim1;
i__3 = ilast + ilast * h_dim1;
z__1.r = h__[i__3].r * signbc.r - h__[i__3].i * signbc.i,
z__1.i = h__[i__3].r * signbc.i + h__[i__3].i *
signbc.r;
h__[i__2].r = z__1.r, h__[i__2].i = z__1.i;
}
if (ilz) {
zscal_(n, &signbc, &z__[ilast * z_dim1 + 1], &c__1);
}
} else {
i__2 = ilast + ilast * t_dim1;
t[i__2].r = 0., t[i__2].i = 0.;
}
i__2 = ilast;
i__3 = ilast + ilast * h_dim1;
alpha[i__2].r = h__[i__3].r, alpha[i__2].i = h__[i__3].i;
i__2 = ilast;
i__3 = ilast + ilast * t_dim1;
beta[i__2].r = t[i__3].r, beta[i__2].i = t[i__3].i;
/* Go to next block -- exit if finished. */
--ilast;
if (ilast < *ilo) {
goto L190;
}
/* Reset counters */
iiter = 0;
eshift.r = 0., eshift.i = 0.;
if (! ilschr) {
ilastm = ilast;
if (ifrstm > ilast) {
ifrstm = *ilo;
}
}
goto L160;
/* QZ step */
/* This iteration only involves rows/columns IFIRST:ILAST. We */
/* assume IFIRST < ILAST, and that the diagonal of B is non-zero. */
L70:
++iiter;
if (! ilschr) {
ifrstm = ifirst;
}
/* Compute the Shift. */
/* At this point, IFIRST < ILAST, and the diagonal elements of */
/* T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in */
/* magnitude) */
if (iiter / 10 * 10 != iiter) {
/* The Wilkinson shift (AEP p.512), i.e., the eigenvalue of */
/* the bottom-right 2x2 block of A inv(B) which is nearest to */
/* the bottom-right element. */
/* We factor B as U*D, where U has unit diagonals, and */
/* compute (A*inv(D))*inv(U). */
i__2 = ilast - 1 + ilast * t_dim1;
z__2.r = bscale * t[i__2].r, z__2.i = bscale * t[i__2].i;
i__3 = ilast + ilast * t_dim1;
z__3.r = bscale * t[i__3].r, z__3.i = bscale * t[i__3].i;
z_div(&z__1, &z__2, &z__3);
u12.r = z__1.r, u12.i = z__1.i;
i__2 = ilast - 1 + (ilast - 1) * h_dim1;
z__2.r = ascale * h__[i__2].r, z__2.i = ascale * h__[i__2].i;
i__3 = ilast - 1 + (ilast - 1) * t_dim1;
z__3.r = bscale * t[i__3].r, z__3.i = bscale * t[i__3].i;
z_div(&z__1, &z__2, &z__3);
ad11.r = z__1.r, ad11.i = z__1.i;
i__2 = ilast + (ilast - 1) * h_dim1;
z__2.r = ascale * h__[i__2].r, z__2.i = ascale * h__[i__2].i;
i__3 = ilast - 1 + (ilast - 1) * t_dim1;
z__3.r = bscale * t[i__3].r, z__3.i = bscale * t[i__3].i;
z_div(&z__1, &z__2, &z__3);
ad21.r = z__1.r, ad21.i = z__1.i;
i__2 = ilast - 1 + ilast * h_dim1;
z__2.r = ascale * h__[i__2].r, z__2.i = ascale * h__[i__2].i;
i__3 = ilast + ilast * t_dim1;
z__3.r = bscale * t[i__3].r, z__3.i = bscale * t[i__3].i;
z_div(&z__1, &z__2, &z__3);
ad12.r = z__1.r, ad12.i = z__1.i;
i__2 = ilast + ilast * h_dim1;
z__2.r = ascale * h__[i__2].r, z__2.i = ascale * h__[i__2].i;
i__3 = ilast + ilast * t_dim1;
z__3.r = bscale * t[i__3].r, z__3.i = bscale * t[i__3].i;
z_div(&z__1, &z__2, &z__3);
ad22.r = z__1.r, ad22.i = z__1.i;
z__2.r = u12.r * ad21.r - u12.i * ad21.i, z__2.i = u12.r * ad21.i
+ u12.i * ad21.r;
z__1.r = ad22.r - z__2.r, z__1.i = ad22.i - z__2.i;
abi22.r = z__1.r, abi22.i = z__1.i;
z__2.r = ad11.r + abi22.r, z__2.i = ad11.i + abi22.i;
z__1.r = z__2.r * .5, z__1.i = z__2.i * .5;
t1.r = z__1.r, t1.i = z__1.i;
pow_zi(&z__4, &t1, &c__2);
z__5.r = ad12.r * ad21.r - ad12.i * ad21.i, z__5.i = ad12.r *
ad21.i + ad12.i * ad21.r;
z__3.r = z__4.r + z__5.r, z__3.i = z__4.i + z__5.i;
z__6.r = ad11.r * ad22.r - ad11.i * ad22.i, z__6.i = ad11.r *
ad22.i + ad11.i * ad22.r;
z__2.r = z__3.r - z__6.r, z__2.i = z__3.i - z__6.i;
z_sqrt(&z__1, &z__2);
rtdisc.r = z__1.r, rtdisc.i = z__1.i;
z__1.r = t1.r - abi22.r, z__1.i = t1.i - abi22.i;
z__2.r = t1.r - abi22.r, z__2.i = t1.i - abi22.i;
temp = z__1.r * rtdisc.r + d_imag(&z__2) * d_imag(&rtdisc);
if (temp <= 0.) {
z__1.r = t1.r + rtdisc.r, z__1.i = t1.i + rtdisc.i;
shift.r = z__1.r, shift.i = z__1.i;
} else {
z__1.r = t1.r - rtdisc.r, z__1.i = t1.i - rtdisc.i;
shift.r = z__1.r, shift.i = z__1.i;
}
} else {
/* Exceptional shift. Chosen for no particularly good reason. */
i__2 = ilast - 1 + ilast * h_dim1;
z__4.r = ascale * h__[i__2].r, z__4.i = ascale * h__[i__2].i;
i__3 = ilast - 1 + (ilast - 1) * t_dim1;
z__5.r = bscale * t[i__3].r, z__5.i = bscale * t[i__3].i;
z_div(&z__3, &z__4, &z__5);
d_cnjg(&z__2, &z__3);
z__1.r = eshift.r + z__2.r, z__1.i = eshift.i + z__2.i;
eshift.r = z__1.r, eshift.i = z__1.i;
shift.r = eshift.r, shift.i = eshift.i;
}
/* Now check for two consecutive small subdiagonals. */
i__2 = ifirst + 1;
for (j = ilast - 1; j >= i__2; --j) {
istart = j;
i__3 = j + j * h_dim1;
z__2.r = ascale * h__[i__3].r, z__2.i = ascale * h__[i__3].i;
i__4 = j + j * t_dim1;
z__4.r = bscale * t[i__4].r, z__4.i = bscale * t[i__4].i;
z__3.r = shift.r * z__4.r - shift.i * z__4.i, z__3.i = shift.r *
z__4.i + shift.i * z__4.r;
z__1.r = z__2.r - z__3.r, z__1.i = z__2.i - z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
temp = (d__1 = ctemp.r, abs(d__1)) + (d__2 = d_imag(&ctemp), abs(
d__2));
i__3 = j + 1 + j * h_dim1;
temp2 = ascale * ((d__1 = h__[i__3].r, abs(d__1)) + (d__2 =
d_imag(&h__[j + 1 + j * h_dim1]), abs(d__2)));
tempr = max(temp,temp2);
if (tempr < 1. && tempr != 0.) {
temp /= tempr;
temp2 /= tempr;
}
i__3 = j + (j - 1) * h_dim1;
if (((d__1 = h__[i__3].r, abs(d__1)) + (d__2 = d_imag(&h__[j + (j
- 1) * h_dim1]), abs(d__2))) * temp2 <= temp * atol) {
goto L90;
}
/* L80: */
}
istart = ifirst;
i__2 = ifirst + ifirst * h_dim1;
z__2.r = ascale * h__[i__2].r, z__2.i = ascale * h__[i__2].i;
i__3 = ifirst + ifirst * t_dim1;
z__4.r = bscale * t[i__3].r, z__4.i = bscale * t[i__3].i;
z__3.r = shift.r * z__4.r - shift.i * z__4.i, z__3.i = shift.r *
z__4.i + shift.i * z__4.r;
z__1.r = z__2.r - z__3.r, z__1.i = z__2.i - z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
L90:
/* Do an implicit-shift QZ sweep. */
/* Initial Q */
i__2 = istart + 1 + istart * h_dim1;
z__1.r = ascale * h__[i__2].r, z__1.i = ascale * h__[i__2].i;
ctemp2.r = z__1.r, ctemp2.i = z__1.i;
zlartg_(&ctemp, &ctemp2, &c__, &s, &ctemp3);
/* Sweep */
i__2 = ilast - 1;
for (j = istart; j <= i__2; ++j) {
if (j > istart) {
i__3 = j + (j - 1) * h_dim1;
ctemp.r = h__[i__3].r, ctemp.i = h__[i__3].i;
zlartg_(&ctemp, &h__[j + 1 + (j - 1) * h_dim1], &c__, &s, &
h__[j + (j - 1) * h_dim1]);
i__3 = j + 1 + (j - 1) * h_dim1;
h__[i__3].r = 0., h__[i__3].i = 0.;
}
i__3 = ilastm;
for (jc = j; jc <= i__3; ++jc) {
i__4 = j + jc * h_dim1;
z__2.r = c__ * h__[i__4].r, z__2.i = c__ * h__[i__4].i;
i__5 = j + 1 + jc * h_dim1;
z__3.r = s.r * h__[i__5].r - s.i * h__[i__5].i, z__3.i = s.r *
h__[i__5].i + s.i * h__[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__4 = j + 1 + jc * h_dim1;
d_cnjg(&z__4, &s);
z__3.r = -z__4.r, z__3.i = -z__4.i;
i__5 = j + jc * h_dim1;
z__2.r = z__3.r * h__[i__5].r - z__3.i * h__[i__5].i, z__2.i =
z__3.r * h__[i__5].i + z__3.i * h__[i__5].r;
i__6 = j + 1 + jc * h_dim1;
z__5.r = c__ * h__[i__6].r, z__5.i = c__ * h__[i__6].i;
z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
i__4 = j + jc * h_dim1;
h__[i__4].r = ctemp.r, h__[i__4].i = ctemp.i;
i__4 = j + jc * t_dim1;
z__2.r = c__ * t[i__4].r, z__2.i = c__ * t[i__4].i;
i__5 = j + 1 + jc * t_dim1;
z__3.r = s.r * t[i__5].r - s.i * t[i__5].i, z__3.i = s.r * t[
i__5].i + s.i * t[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp2.r = z__1.r, ctemp2.i = z__1.i;
i__4 = j + 1 + jc * t_dim1;
d_cnjg(&z__4, &s);
z__3.r = -z__4.r, z__3.i = -z__4.i;
i__5 = j + jc * t_dim1;
z__2.r = z__3.r * t[i__5].r - z__3.i * t[i__5].i, z__2.i =
z__3.r * t[i__5].i + z__3.i * t[i__5].r;
i__6 = j + 1 + jc * t_dim1;
z__5.r = c__ * t[i__6].r, z__5.i = c__ * t[i__6].i;
z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
t[i__4].r = z__1.r, t[i__4].i = z__1.i;
i__4 = j + jc * t_dim1;
t[i__4].r = ctemp2.r, t[i__4].i = ctemp2.i;
/* L100: */
}
if (ilq) {
i__3 = *n;
for (jr = 1; jr <= i__3; ++jr) {
i__4 = jr + j * q_dim1;
z__2.r = c__ * q[i__4].r, z__2.i = c__ * q[i__4].i;
d_cnjg(&z__4, &s);
i__5 = jr + (j + 1) * q_dim1;
z__3.r = z__4.r * q[i__5].r - z__4.i * q[i__5].i, z__3.i =
z__4.r * q[i__5].i + z__4.i * q[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__4 = jr + (j + 1) * q_dim1;
z__3.r = -s.r, z__3.i = -s.i;
i__5 = jr + j * q_dim1;
z__2.r = z__3.r * q[i__5].r - z__3.i * q[i__5].i, z__2.i =
z__3.r * q[i__5].i + z__3.i * q[i__5].r;
i__6 = jr + (j + 1) * q_dim1;
z__4.r = c__ * q[i__6].r, z__4.i = c__ * q[i__6].i;
z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
q[i__4].r = z__1.r, q[i__4].i = z__1.i;
i__4 = jr + j * q_dim1;
q[i__4].r = ctemp.r, q[i__4].i = ctemp.i;
/* L110: */
}
}
i__3 = j + 1 + (j + 1) * t_dim1;
ctemp.r = t[i__3].r, ctemp.i = t[i__3].i;
zlartg_(&ctemp, &t[j + 1 + j * t_dim1], &c__, &s, &t[j + 1 + (j +
1) * t_dim1]);
i__3 = j + 1 + j * t_dim1;
t[i__3].r = 0., t[i__3].i = 0.;
/* Computing MIN */
i__4 = j + 2;
i__3 = min(i__4,ilast);
for (jr = ifrstm; jr <= i__3; ++jr) {
i__4 = jr + (j + 1) * h_dim1;
z__2.r = c__ * h__[i__4].r, z__2.i = c__ * h__[i__4].i;
i__5 = jr + j * h_dim1;
z__3.r = s.r * h__[i__5].r - s.i * h__[i__5].i, z__3.i = s.r *
h__[i__5].i + s.i * h__[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__4 = jr + j * h_dim1;
d_cnjg(&z__4, &s);
z__3.r = -z__4.r, z__3.i = -z__4.i;
i__5 = jr + (j + 1) * h_dim1;
z__2.r = z__3.r * h__[i__5].r - z__3.i * h__[i__5].i, z__2.i =
z__3.r * h__[i__5].i + z__3.i * h__[i__5].r;
i__6 = jr + j * h_dim1;
z__5.r = c__ * h__[i__6].r, z__5.i = c__ * h__[i__6].i;
z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
i__4 = jr + (j + 1) * h_dim1;
h__[i__4].r = ctemp.r, h__[i__4].i = ctemp.i;
/* L120: */
}
i__3 = j;
for (jr = ifrstm; jr <= i__3; ++jr) {
i__4 = jr + (j + 1) * t_dim1;
z__2.r = c__ * t[i__4].r, z__2.i = c__ * t[i__4].i;
i__5 = jr + j * t_dim1;
z__3.r = s.r * t[i__5].r - s.i * t[i__5].i, z__3.i = s.r * t[
i__5].i + s.i * t[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__4 = jr + j * t_dim1;
d_cnjg(&z__4, &s);
z__3.r = -z__4.r, z__3.i = -z__4.i;
i__5 = jr + (j + 1) * t_dim1;
z__2.r = z__3.r * t[i__5].r - z__3.i * t[i__5].i, z__2.i =
z__3.r * t[i__5].i + z__3.i * t[i__5].r;
i__6 = jr + j * t_dim1;
z__5.r = c__ * t[i__6].r, z__5.i = c__ * t[i__6].i;
z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
t[i__4].r = z__1.r, t[i__4].i = z__1.i;
i__4 = jr + (j + 1) * t_dim1;
t[i__4].r = ctemp.r, t[i__4].i = ctemp.i;
/* L130: */
}
if (ilz) {
i__3 = *n;
for (jr = 1; jr <= i__3; ++jr) {
i__4 = jr + (j + 1) * z_dim1;
z__2.r = c__ * z__[i__4].r, z__2.i = c__ * z__[i__4].i;
i__5 = jr + j * z_dim1;
z__3.r = s.r * z__[i__5].r - s.i * z__[i__5].i, z__3.i =
s.r * z__[i__5].i + s.i * z__[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__4 = jr + j * z_dim1;
d_cnjg(&z__4, &s);
z__3.r = -z__4.r, z__3.i = -z__4.i;
i__5 = jr + (j + 1) * z_dim1;
z__2.r = z__3.r * z__[i__5].r - z__3.i * z__[i__5].i,
z__2.i = z__3.r * z__[i__5].i + z__3.i * z__[i__5]
.r;
i__6 = jr + j * z_dim1;
z__5.r = c__ * z__[i__6].r, z__5.i = c__ * z__[i__6].i;
z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
z__[i__4].r = z__1.r, z__[i__4].i = z__1.i;
i__4 = jr + (j + 1) * z_dim1;
z__[i__4].r = ctemp.r, z__[i__4].i = ctemp.i;
/* L140: */
}
}
/* L150: */
}
L160:
/* L170: */
;
}
/* Drop-through = non-convergence */
L180:
*info = ilast;
goto L210;
/* Successful completion of all QZ steps */
L190:
/* Set Eigenvalues 1:ILO-1 */
i__1 = *ilo - 1;
for (j = 1; j <= i__1; ++j) {
absb = z_abs(&t[j + j * t_dim1]);
if (absb > safmin) {
i__2 = j + j * t_dim1;
z__2.r = t[i__2].r / absb, z__2.i = t[i__2].i / absb;
d_cnjg(&z__1, &z__2);
signbc.r = z__1.r, signbc.i = z__1.i;
i__2 = j + j * t_dim1;
t[i__2].r = absb, t[i__2].i = 0.;
if (ilschr) {
i__2 = j - 1;
zscal_(&i__2, &signbc, &t[j * t_dim1 + 1], &c__1);
zscal_(&j, &signbc, &h__[j * h_dim1 + 1], &c__1);
} else {
i__2 = j + j * h_dim1;
i__3 = j + j * h_dim1;
z__1.r = h__[i__3].r * signbc.r - h__[i__3].i * signbc.i,
z__1.i = h__[i__3].r * signbc.i + h__[i__3].i *
signbc.r;
h__[i__2].r = z__1.r, h__[i__2].i = z__1.i;
}
if (ilz) {
zscal_(n, &signbc, &z__[j * z_dim1 + 1], &c__1);
}
} else {
i__2 = j + j * t_dim1;
t[i__2].r = 0., t[i__2].i = 0.;
}
i__2 = j;
i__3 = j + j * h_dim1;
alpha[i__2].r = h__[i__3].r, alpha[i__2].i = h__[i__3].i;
i__2 = j;
i__3 = j + j * t_dim1;
beta[i__2].r = t[i__3].r, beta[i__2].i = t[i__3].i;
/* L200: */
}
/* Normal Termination */
*info = 0;
/* Exit (other than argument error) -- return optimal workspace size */
L210:
z__1.r = (doublereal) (*n), z__1.i = 0.;
work[1].r = z__1.r, work[1].i = z__1.i;
return 0;
/* End of ZHGEQZ */
} /* zhgeqz_ */
