/* qzval.f -- translated by f2c (version 19961017).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Subroutine */ int qzval_(integer *nm, integer *n, doublereal *a,
doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta,
logical *matz, doublereal *z__)
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i__1, i__2;
doublereal d__1, d__2, d__3, d__4;
/* Builtin functions */
double sqrt(doublereal), d_sign(doublereal *, doublereal *);
/* Local variables */
static doublereal epsb, c__, d__, e;
static integer i__, j;
static doublereal r__, s, t, a1, a2, u1, u2, v1, v2, a11, a12, a21, a22,
b11, b12, b22, di, ei;
static integer na;
static doublereal an, bn;
static integer en;
static doublereal cq, dr;
static integer nn;
static doublereal cz, ti, tr, a1i, a2i, a11i, a12i, a22i, a11r, a12r,
a22r, sqi, ssi;
static integer isw;
static doublereal sqr, szi, ssr, szr;
/* THIS SUBROUTINE IS THE THIRD STEP OF THE QZ ALGORITHM */
/* FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */
/* SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. */
/* THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM */
/* IN QUASI-TRIANGULAR FORM AND THE OTHER IN UPPER TRIANGULAR FORM. */
/* IT REDUCES THE QUASI-TRIANGULAR MATRIX FURTHER, SO THAT ANY */
/* REMAINING 2-BY-2 BLOCKS CORRESPOND TO PAIRS OF COMPLEX */
/* EIGENVALUES, AND RETURNS QUANTITIES WHOSE RATIOS GIVE THE */
/* GENERALIZED EIGENVALUES. IT IS USUALLY PRECEDED BY QZHES */
/* AND QZIT AND MAY BE FOLLOWED BY QZVEC. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRICES. */
/* A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX. */
/* B CONTAINS A REAL UPPER TRIANGULAR MATRIX. IN ADDITION, */
/* LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) */
/* COMPUTED AND SAVED IN QZIT. */
/* MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS
*/
/* ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING */
/* EIGENVECTORS, AND TO .FALSE. OTHERWISE. */
/* Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE */
/* TRANSFORMATION MATRIX PRODUCED IN THE REDUCTIONS BY QZHES */
/* AND QZIT, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. */
/* IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. */
/* ON OUTPUT */
/* A HAS BEEN REDUCED FURTHER TO A QUASI-TRIANGULAR MATRIX */
/* IN WHICH ALL NONZERO SUBDIAGONAL ELEMENTS CORRESPOND TO */
/* PAIRS OF COMPLEX EIGENVALUES. */
/* B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS */
/* HAVE BEEN ALTERED. B(N,1) IS UNALTERED. */
/* ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS OF THE */
/* DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX THAT WOULD BE */
/* OBTAINED IF A WERE REDUCED COMPLETELY TO TRIANGULAR FORM */
/* BY UNITARY TRANSFORMATIONS. NON-ZERO VALUES OF ALFI OCCUR */
/* IN PAIRS, THE FIRST MEMBER POSITIVE AND THE SECOND NEGATIVE.
*/
/* BETA CONTAINS THE DIAGONAL ELEMENTS OF THE CORRESPONDING B, */
/* NORMALIZED TO BE REAL AND NON-NEGATIVE. THE GENERALIZED */
/* EIGENVALUES ARE THEN THE RATIOS ((ALFR+I*ALFI)/BETA). */
/* Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS */
/* (FOR ALL THREE STEPS) IF MATZ HAS BEEN SET TO .TRUE. */
/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
*/
/* THIS VERSION DATED AUGUST 1983. */
/* ------------------------------------------------------------------
*/
/* Parameter adjustments */
z_dim1 = *nm;
z_offset = z_dim1 + 1;
z__ -= z_offset;
--beta;
--alfi;
--alfr;
b_dim1 = *nm;
b_offset = b_dim1 + 1;
b -= b_offset;
a_dim1 = *nm;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
epsb = b[*n + b_dim1];
isw = 1;
/* .......... FIND EIGENVALUES OF QUASI-TRIANGULAR MATRICES. */
/* FOR EN=N STEP -1 UNTIL 1 DO -- .......... */
i__1 = *n;
for (nn = 1; nn <= i__1; ++nn) {
en = *n + 1 - nn;
na = en - 1;
if (isw == 2) {
goto L505;
}
if (en == 1) {
goto L410;
}
if (a[en + na * a_dim1] != 0.) {
goto L420;
}
/* .......... 1-BY-1 BLOCK, ONE REAL ROOT .......... */
L410:
alfr[en] = a[en + en * a_dim1];
if (b[en + en * b_dim1] < 0.) {
alfr[en] = -alfr[en];
}
beta[en] = (d__1 = b[en + en * b_dim1], abs(d__1));
alfi[en] = 0.;
goto L510;
/* .......... 2-BY-2 BLOCK .......... */
L420:
if ((d__1 = b[na + na * b_dim1], abs(d__1)) <= epsb) {
goto L455;
}
if ((d__1 = b[en + en * b_dim1], abs(d__1)) > epsb) {
goto L430;
}
a1 = a[en + en * a_dim1];
a2 = a[en + na * a_dim1];
bn = 0.;
goto L435;
L430:
an = (d__1 = a[na + na * a_dim1], abs(d__1)) + (d__2 = a[na + en *
a_dim1], abs(d__2)) + (d__3 = a[en + na * a_dim1], abs(d__3))
+ (d__4 = a[en + en * a_dim1], abs(d__4));
bn = (d__1 = b[na + na * b_dim1], abs(d__1)) + (d__2 = b[na + en *
b_dim1], abs(d__2)) + (d__3 = b[en + en * b_dim1], abs(d__3));
a11 = a[na + na * a_dim1] / an;
a12 = a[na + en * a_dim1] / an;
a21 = a[en + na * a_dim1] / an;
a22 = a[en + en * a_dim1] / an;
b11 = b[na + na * b_dim1] / bn;
b12 = b[na + en * b_dim1] / bn;
b22 = b[en + en * b_dim1] / bn;
e = a11 / b11;
ei = a22 / b22;
s = a21 / (b11 * b22);
t = (a22 - e * b22) / b22;
if (abs(e) <= abs(ei)) {
goto L431;
}
e = ei;
t = (a11 - e * b11) / b11;
L431:
c__ = (t - s * b12) * .5;
d__ = c__ * c__ + s * (a12 - e * b12);
if (d__ < 0.) {
goto L480;
}
/* .......... TWO REAL ROOTS. */
/* ZERO BOTH A(EN,NA) AND B(EN,NA) .......... */
d__1 = sqrt(d__);
e += c__ + d_sign(&d__1, &c__);
a11 -= e * b11;
a12 -= e * b12;
a22 -= e * b22;
if (abs(a11) + abs(a12) < abs(a21) + abs(a22)) {
goto L432;
}
a1 = a12;
a2 = a11;
goto L435;
L432:
a1 = a22;
a2 = a21;
/* .......... CHOOSE AND APPLY REAL Z .......... */
L435:
s = abs(a1) + abs(a2);
u1 = a1 / s;
u2 = a2 / s;
d__1 = sqrt(u1 * u1 + u2 * u2);
r__ = d_sign(&d__1, &u1);
v1 = -(u1 + r__) / r__;
v2 = -u2 / r__;
u2 = v2 / v1;
i__2 = en;
for (i__ = 1; i__ <= i__2; ++i__) {
t = a[i__ + en * a_dim1] + u2 * a[i__ + na * a_dim1];
a[i__ + en * a_dim1] += t * v1;
a[i__ + na * a_dim1] += t * v2;
t = b[i__ + en * b_dim1] + u2 * b[i__ + na * b_dim1];
b[i__ + en * b_dim1] += t * v1;
b[i__ + na * b_dim1] += t * v2;
/* L440: */
}
if (! (*matz)) {
goto L450;
}
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
t = z__[i__ + en * z_dim1] + u2 * z__[i__ + na * z_dim1];
z__[i__ + en * z_dim1] += t * v1;
z__[i__ + na * z_dim1] += t * v2;
/* L445: */
}
L450:
if (bn == 0.) {
goto L475;
}
if (an < abs(e) * bn) {
goto L455;
}
a1 = b[na + na * b_dim1];
a2 = b[en + na * b_dim1];
goto L460;
L455:
a1 = a[na + na * a_dim1];
a2 = a[en + na * a_dim1];
/* .......... CHOOSE AND APPLY REAL Q .......... */
L460:
s = abs(a1) + abs(a2);
if (s == 0.) {
goto L475;
}
u1 = a1 / s;
u2 = a2 / s;
d__1 = sqrt(u1 * u1 + u2 * u2);
r__ = d_sign(&d__1, &u1);
v1 = -(u1 + r__) / r__;
v2 = -u2 / r__;
u2 = v2 / v1;
i__2 = *n;
for (j = na; j <= i__2; ++j) {
t = a[na + j * a_dim1] + u2 * a[en + j * a_dim1];
a[na + j * a_dim1] += t * v1;
a[en + j * a_dim1] += t * v2;
t = b[na + j * b_dim1] + u2 * b[en + j * b_dim1];
b[na + j * b_dim1] += t * v1;
b[en + j * b_dim1] += t * v2;
/* L470: */
}
L475:
a[en + na * a_dim1] = 0.;
b[en + na * b_dim1] = 0.;
alfr[na] = a[na + na * a_dim1];
alfr[en] = a[en + en * a_dim1];
if (b[na + na * b_dim1] < 0.) {
alfr[na] = -alfr[na];
}
if (b[en + en * b_dim1] < 0.) {
alfr[en] = -alfr[en];
}
beta[na] = (d__1 = b[na + na * b_dim1], abs(d__1));
beta[en] = (d__1 = b[en + en * b_dim1], abs(d__1));
alfi[en] = 0.;
alfi[na] = 0.;
goto L505;
/* .......... TWO COMPLEX ROOTS .......... */
L480:
e += c__;
ei = sqrt(-d__);
a11r = a11 - e * b11;
a11i = ei * b11;
a12r = a12 - e * b12;
a12i = ei * b12;
a22r = a22 - e * b22;
a22i = ei * b22;
if (abs(a11r) + abs(a11i) + abs(a12r) + abs(a12i) < abs(a21) + abs(
a22r) + abs(a22i)) {
goto L482;
}
a1 = a12r;
a1i = a12i;
a2 = -a11r;
a2i = -a11i;
goto L485;
L482:
a1 = a22r;
a1i = a22i;
a2 = -a21;
a2i = 0.;
/* .......... CHOOSE COMPLEX Z .......... */
L485:
cz = sqrt(a1 * a1 + a1i * a1i);
if (cz == 0.) {
goto L487;
}
szr = (a1 * a2 + a1i * a2i) / cz;
szi = (a1 * a2i - a1i * a2) / cz;
r__ = sqrt(cz * cz + szr * szr + szi * szi);
cz /= r__;
szr /= r__;
szi /= r__;
goto L490;
L487:
szr = 1.;
szi = 0.;
L490:
if (an < (abs(e) + ei) * bn) {
goto L492;
}
a1 = cz * b11 + szr * b12;
a1i = szi * b12;
a2 = szr * b22;
a2i = szi * b22;
goto L495;
L492:
a1 = cz * a11 + szr * a12;
a1i = szi * a12;
a2 = cz * a21 + szr * a22;
a2i = szi * a22;
/* .......... CHOOSE COMPLEX Q .......... */
L495:
cq = sqrt(a1 * a1 + a1i * a1i);
if (cq == 0.) {
goto L497;
}
sqr = (a1 * a2 + a1i * a2i) / cq;
sqi = (a1 * a2i - a1i * a2) / cq;
r__ = sqrt(cq * cq + sqr * sqr + sqi * sqi);
cq /= r__;
sqr /= r__;
sqi /= r__;
goto L500;
L497:
sqr = 1.;
sqi = 0.;
/* .......... COMPUTE DIAGONAL ELEMENTS THAT WOULD RESULT */
/* IF TRANSFORMATIONS WERE APPLIED .......... */
L500:
ssr = sqr * szr + sqi * szi;
ssi = sqr * szi - sqi * szr;
i__ = 1;
tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21 + ssr * a22;
ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22;
dr = cq * cz * b11 + cq * szr * b12 + ssr * b22;
di = cq * szi * b12 + ssi * b22;
goto L503;
L502:
i__ = 2;
tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21 + cq * cz * a22;
ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21;
dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22;
di = -ssi * b11 - sqi * cz * b12;
L503:
t = ti * dr - tr * di;
j = na;
if (t < 0.) {
j = en;
}
r__ = sqrt(dr * dr + di * di);
beta[j] = bn * r__;
alfr[j] = an * (tr * dr + ti * di) / r__;
alfi[j] = an * t / r__;
if (i__ == 1) {
goto L502;
}
L505:
isw = 3 - isw;
L510:
;
}
b[*n + b_dim1] = epsb;
return 0;
} /* qzval_ */
syntax highlighted by Code2HTML, v. 0.9.1