/* qzvec.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 qzvec_(integer *nm, integer *n, doublereal *a,
doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta,
doublereal *z__)
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i__1, i__2,
i__3;
doublereal d__1, d__2;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static doublereal alfm, almi, betm, epsb, almr, d__;
static integer i__, j, k, m;
static doublereal q, r__, s, t, w, x, y, t1, t2, w1, x1, z1, di;
static integer na, ii, en, jj;
static doublereal ra, dr, sa;
static integer nn;
static doublereal ti, rr, tr, zz;
static integer isw, enm2;
/* THIS SUBROUTINE IS THE OPTIONAL FOURTH 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 (IN WHICH EACH 2-BY-2 BLOCK CORRESPONDS TO */
/* A PAIR OF COMPLEX EIGENVALUES) AND THE OTHER IN UPPER TRIANGULAR */
/* FORM. IT COMPUTES THE EIGENVECTORS OF THE TRIANGULAR PROBLEM AND
*/
/* TRANSFORMS THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM. */
/* IT IS USUALLY PRECEDED BY QZHES, QZIT, AND QZVAL. */
/* 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. */
/* ALFR, ALFI, AND BETA ARE VECTORS WITH COMPONENTS WHOSE */
/* RATIOS ((ALFR+I*ALFI)/BETA) ARE THE GENERALIZED */
/* EIGENVALUES. THEY ARE USUALLY OBTAINED FROM QZVAL. */
/* Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE */
/* REDUCTIONS BY QZHES, QZIT, AND QZVAL, IF PERFORMED. */
/* IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM ARE */
/* DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX. */
/* ON OUTPUT */
/* A IS UNALTERED. ITS SUBDIAGONAL ELEMENTS PROVIDE INFORMATION */
/* ABOUT THE STORAGE OF THE COMPLEX EIGENVECTORS. */
/* B HAS BEEN DESTROYED. */
/* ALFR, ALFI, AND BETA ARE UNALTERED. */
/* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. */
/* IF ALFI(I) .EQ. 0.0, THE I-TH EIGENVALUE IS REAL AND */
/* THE I-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR. */
/* IF ALFI(I) .NE. 0.0, THE I-TH EIGENVALUE IS COMPLEX. */
/* IF ALFI(I) .GT. 0.0, THE EIGENVALUE IS THE FIRST OF */
/* A COMPLEX PAIR AND THE I-TH AND (I+1)-TH COLUMNS */
/* OF Z CONTAIN ITS EIGENVECTOR. */
/* IF ALFI(I) .LT. 0.0, THE EIGENVALUE IS THE SECOND OF */
/* A COMPLEX PAIR AND THE (I-1)-TH AND I-TH COLUMNS */
/* OF Z CONTAIN THE CONJUGATE OF ITS EIGENVECTOR. */
/* EACH EIGENVECTOR IS NORMALIZED SO THAT THE MODULUS */
/* OF ITS LARGEST COMPONENT IS 1.0 . */
/* 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;
/* .......... 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 L795;
}
if (alfi[en] != 0.) {
goto L710;
}
/* .......... REAL VECTOR .......... */
m = en;
b[en + en * b_dim1] = 1.;
if (na == 0) {
goto L800;
}
alfm = alfr[m];
betm = beta[m];
/* .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */
i__2 = na;
for (ii = 1; ii <= i__2; ++ii) {
i__ = en - ii;
w = betm * a[i__ + i__ * a_dim1] - alfm * b[i__ + i__ * b_dim1];
r__ = 0.;
i__3 = en;
for (j = m; j <= i__3; ++j) {
/* L610: */
r__ += (betm * a[i__ + j * a_dim1] - alfm * b[i__ + j *
b_dim1]) * b[j + en * b_dim1];
}
if (i__ == 1 || isw == 2) {
goto L630;
}
if (betm * a[i__ + (i__ - 1) * a_dim1] == 0.) {
goto L630;
}
zz = w;
s = r__;
goto L690;
L630:
m = i__;
if (isw == 2) {
goto L640;
}
/* .......... REAL 1-BY-1 BLOCK .......... */
t = w;
if (w == 0.) {
t = epsb;
}
b[i__ + en * b_dim1] = -r__ / t;
goto L700;
/* .......... REAL 2-BY-2 BLOCK .......... */
L640:
x = betm * a[i__ + (i__ + 1) * a_dim1] - alfm * b[i__ + (i__ + 1)
* b_dim1];
y = betm * a[i__ + 1 + i__ * a_dim1];
q = w * zz - x * y;
t = (x * s - zz * r__) / q;
b[i__ + en * b_dim1] = t;
if (abs(x) <= abs(zz)) {
goto L650;
}
b[i__ + 1 + en * b_dim1] = (-r__ - w * t) / x;
goto L690;
L650:
b[i__ + 1 + en * b_dim1] = (-s - y * t) / zz;
L690:
isw = 3 - isw;
L700:
;
}
/* .......... END REAL VECTOR .......... */
goto L800;
/* .......... COMPLEX VECTOR .......... */
L710:
m = na;
almr = alfr[m];
almi = alfi[m];
betm = beta[m];
/* .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT */
/* EIGENVECTOR MATRIX IS TRIANGULAR .......... */
y = betm * a[en + na * a_dim1];
b[na + na * b_dim1] = -almi * b[en + en * b_dim1] / y;
b[na + en * b_dim1] = (almr * b[en + en * b_dim1] - betm * a[en + en *
a_dim1]) / y;
b[en + na * b_dim1] = 0.;
b[en + en * b_dim1] = 1.;
enm2 = na - 1;
if (enm2 == 0) {
goto L795;
}
/* .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- .......... */
i__2 = enm2;
for (ii = 1; ii <= i__2; ++ii) {
i__ = na - ii;
w = betm * a[i__ + i__ * a_dim1] - almr * b[i__ + i__ * b_dim1];
w1 = -almi * b[i__ + i__ * b_dim1];
ra = 0.;
sa = 0.;
i__3 = en;
for (j = m; j <= i__3; ++j) {
x = betm * a[i__ + j * a_dim1] - almr * b[i__ + j * b_dim1];
x1 = -almi * b[i__ + j * b_dim1];
ra = ra + x * b[j + na * b_dim1] - x1 * b[j + en * b_dim1];
sa = sa + x * b[j + en * b_dim1] + x1 * b[j + na * b_dim1];
/* L760: */
}
if (i__ == 1 || isw == 2) {
goto L770;
}
if (betm * a[i__ + (i__ - 1) * a_dim1] == 0.) {
goto L770;
}
zz = w;
z1 = w1;
r__ = ra;
s = sa;
isw = 2;
goto L790;
L770:
m = i__;
if (isw == 2) {
goto L780;
}
/* .......... COMPLEX 1-BY-1 BLOCK .......... */
tr = -ra;
ti = -sa;
L773:
dr = w;
di = w1;
/* .......... COMPLEX DIVIDE (T1,T2) = (TR,TI) / (DR,DI) .....
..... */
L775:
if (abs(di) > abs(dr)) {
goto L777;
}
rr = di / dr;
d__ = dr + di * rr;
t1 = (tr + ti * rr) / d__;
t2 = (ti - tr * rr) / d__;
switch (isw) {
case 1: goto L787;
case 2: goto L782;
}
L777:
rr = dr / di;
d__ = dr * rr + di;
t1 = (tr * rr + ti) / d__;
t2 = (ti * rr - tr) / d__;
switch (isw) {
case 1: goto L787;
case 2: goto L782;
}
/* .......... COMPLEX 2-BY-2 BLOCK .......... */
L780:
x = betm * a[i__ + (i__ + 1) * a_dim1] - almr * b[i__ + (i__ + 1)
* b_dim1];
x1 = -almi * b[i__ + (i__ + 1) * b_dim1];
y = betm * a[i__ + 1 + i__ * a_dim1];
tr = y * ra - w * r__ + w1 * s;
ti = y * sa - w * s - w1 * r__;
dr = w * zz - w1 * z1 - x * y;
di = w * z1 + w1 * zz - x1 * y;
if (dr == 0. && di == 0.) {
dr = epsb;
}
goto L775;
L782:
b[i__ + 1 + na * b_dim1] = t1;
b[i__ + 1 + en * b_dim1] = t2;
isw = 1;
if (abs(y) > abs(w) + abs(w1)) {
goto L785;
}
tr = -ra - x * b[i__ + 1 + na * b_dim1] + x1 * b[i__ + 1 + en *
b_dim1];
ti = -sa - x * b[i__ + 1 + en * b_dim1] - x1 * b[i__ + 1 + na *
b_dim1];
goto L773;
L785:
t1 = (-r__ - zz * b[i__ + 1 + na * b_dim1] + z1 * b[i__ + 1 + en *
b_dim1]) / y;
t2 = (-s - zz * b[i__ + 1 + en * b_dim1] - z1 * b[i__ + 1 + na *
b_dim1]) / y;
L787:
b[i__ + na * b_dim1] = t1;
b[i__ + en * b_dim1] = t2;
L790:
;
}
/* .......... END COMPLEX VECTOR .......... */
L795:
isw = 3 - isw;
L800:
;
}
/* .......... END BACK SUBSTITUTION. */
/* TRANSFORM TO ORIGINAL COORDINATE SYSTEM. */
/* FOR J=N STEP -1 UNTIL 1 DO -- .......... */
i__1 = *n;
for (jj = 1; jj <= i__1; ++jj) {
j = *n + 1 - jj;
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
zz = 0.;
i__3 = j;
for (k = 1; k <= i__3; ++k) {
/* L860: */
zz += z__[i__ + k * z_dim1] * b[k + j * b_dim1];
}
z__[i__ + j * z_dim1] = zz;
/* L880: */
}
}
/* .......... NORMALIZE SO THAT MODULUS OF LARGEST */
/* COMPONENT OF EACH VECTOR IS 1. */
/* (ISW IS 1 INITIALLY FROM BEFORE) .......... */
i__2 = *n;
for (j = 1; j <= i__2; ++j) {
d__ = 0.;
if (isw == 2) {
goto L920;
}
if (alfi[j] != 0.) {
goto L945;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
if ((d__1 = z__[i__ + j * z_dim1], abs(d__1)) > d__) {
d__ = (d__2 = z__[i__ + j * z_dim1], abs(d__2));
}
/* L890: */
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/* L900: */
z__[i__ + j * z_dim1] /= d__;
}
goto L950;
L920:
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
r__ = (d__1 = z__[i__ + (j - 1) * z_dim1], abs(d__1)) + (d__2 =
z__[i__ + j * z_dim1], abs(d__2));
if (r__ != 0.) {
/* Computing 2nd power */
d__1 = z__[i__ + (j - 1) * z_dim1] / r__;
/* Computing 2nd power */
d__2 = z__[i__ + j * z_dim1] / r__;
r__ *= sqrt(d__1 * d__1 + d__2 * d__2);
}
if (r__ > d__) {
d__ = r__;
}
/* L930: */
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
z__[i__ + (j - 1) * z_dim1] /= d__;
z__[i__ + j * z_dim1] /= d__;
/* L940: */
}
L945:
isw = 3 - isw;
L950:
;
}
return 0;
} /* qzvec_ */
syntax highlighted by Code2HTML, v. 0.9.1