Show eispack.cpp syntax highlighted
/* eispack.f -- translated by f2c (version 19970805).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "Oscill8External.h" // emery added for oscill8
#include "auto_f2c.h"
#include "math.h"
#include "auto_c.h"
BEGIN_AUTO_NAMESPACE;
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Eigenvalue solver from EISPACK */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
rg(integer nm, integer n, doublereal *a, doublereal *wr, doublereal *wi, integer matz, doublereal *z__, integer *iv1, doublereal *fv1, integer *ierr)
{
/* System generated locals */
integer a_dim1, a_offset, z_dim1, z_offset;
/* Local variables */
integer is1, is2;
/* THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF */
/* SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) */
/* TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) */
/* OF A REAL GENERAL MATRIX. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRIX A. */
/* A CONTAINS THE REAL GENERAL MATRIX. */
/* MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF */
/* ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO */
/* ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. */
/* ON OUTPUT */
/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
/* RESPECTIVELY, OF THE EIGENVALUES. COMPLEX CONJUGATE */
/* PAIRS OF EIGENVALUES APPEAR CONSECUTIVELY WITH THE */
/* EIGENVALUE HAVING THE POSITIVE IMAGINARY PART FIRST. */
/* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS */
/* IF MATZ IS NOT ZERO. IF THE J-TH EIGENVALUE IS REAL, THE */
/* J-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR. IF THE J-TH */
/* EIGENVALUE IS COMPLEX WITH POSITIVE IMAGINARY PART, THE */
/* J-TH AND (J+1)-TH COLUMNS OF Z CONTAIN THE REAL AND */
/* IMAGINARY PARTS OF ITS EIGENVECTOR. THE CONJUGATE OF THIS */
/* VECTOR IS THE EIGENVECTOR FOR THE CONJUGATE EIGENVALUE. */
/* IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR */
/* COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR HQR */
/* AND HQR2. THE NORMAL COMPLETION CODE IS ZERO. */
/* IV1 AND FV1 ARE TEMPORARY STORAGE ARRAYS. */
/* 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 */
--fv1;
--iv1;
z_dim1 = nm;
z_offset = z_dim1 + 1;
z__ -= z_offset;
--wi;
--wr;
a_dim1 = nm;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
if (n <= nm) {
goto L10;
}
*ierr = n * 10;
goto L50;
L10:
balanc(&nm, &n, &a[a_offset], &is1, &is2, &fv1[1]);
elmhes(&nm, &n, &is1, &is2, &a[a_offset], &iv1[1]);
if (matz != 0) {
goto L20;
}
/* .......... FIND EIGENVALUES ONLY .......... */
hqr(&nm, &n, &is1, &is2, &a[a_offset], &wr[1], &wi[1], ierr);
goto L50;
/* .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... */
L20:
eltran(&nm, &n, &is1, &is2, &a[a_offset], &iv1[1], &z__[z_offset]);
hqr2(&nm, &n, &is1, &is2, &a[a_offset], &wr[1], &wi[1], &z__[z_offset],
ierr);
if (*ierr != 0) {
goto L50;
}
balbak(&nm, &n, &is1, &is2, &fv1[1], &n, &z__[z_offset]);
L50:
return 0;
} /* rg_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
hqr(integer *nm, integer *n, integer *low, integer *igh, doublereal *h__, doublereal *wr, doublereal *wi, integer *ierr)
{
/* System generated locals */
integer h_dim1, h_offset, i__1, i__2, i__3;
doublereal d__1, d__2;
/* Local variables */
doublereal norm;
integer i__, j, k, l, m;
doublereal p, q, r__, s, t, w, x, y;
integer na, en, ll, mm;
doublereal zz;
logical notlas;
integer mp2, itn, its, enm2;
doublereal tst1, tst2;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, */
/* NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). */
/* THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL */
/* UPPER HESSENBERG MATRIX BY THE QR METHOD. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRIX. */
/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
/* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */
/* SET LOW=1, IGH=N. */
/* H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT */
/* THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG */
/* FORM BY ELMHES OR ORTHES, IF PERFORMED, IS STORED */
/* IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. */
/* ON OUTPUT */
/* H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED */
/* BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND */
/* BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED. */
/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
/* RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES */
/* ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS */
/* OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE */
/* HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN */
/* ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
/* FOR INDICES IERR+1,...,N. */
/* IERR IS SET TO */
/* ZERO FOR NORMAL RETURN, */
/* J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
/* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
/* 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 */
--wi;
--wr;
h_dim1 = *nm;
h_offset = h_dim1 + 1;
h__ -= h_offset;
/* Function Body */
*ierr = 0;
norm = 0.;
k = 1;
/* .......... STORE ROOTS ISOLATED BY BALANC */
/* AND COMPUTE MATRIX NORM .......... */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = *n;
for (j = k; j <= i__2; ++j) {
/* L40: */
norm += (d__1 = h__[i__ + j * h_dim1], fabs(d__1));
}
k = i__;
if (i__ >= *low && i__ <= *igh) {
goto L50;
}
wr[i__] = h__[i__ + i__ * h_dim1];
wi[i__] = 0.;
L50:
;
}
en = *igh;
t = 0.;
itn = *n * 30;
/* .......... SEARCH FOR NEXT EIGENVALUES .......... */
L60:
if (en < *low) {
goto L1001;
}
its = 0;
na = en - 1;
enm2 = na - 1;
/* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
/* FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
L70:
i__1 = en;
for (ll = *low; ll <= i__1; ++ll) {
l = en + *low - ll;
if (l == *low) {
goto L100;
}
s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], fabs(d__1)) + (d__2 = h__[l
+ l * h_dim1], fabs(d__2));
if (s == 0.) {
s = norm;
}
tst1 = s;
tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], fabs(d__1));
if (tst2 == tst1) {
goto L100;
}
/* L80: */
}
/* .......... FORM SHIFT .......... */
L100:
x = h__[en + en * h_dim1];
if (l == en) {
goto L270;
}
y = h__[na + na * h_dim1];
w = h__[en + na * h_dim1] * h__[na + en * h_dim1];
if (l == na) {
goto L280;
}
if (itn == 0) {
goto L1000;
}
if (its != 10 && its != 20) {
goto L130;
}
/* .......... FORM EXCEPTIONAL SHIFT .......... */
t += x;
i__1 = en;
for (i__ = *low; i__ <= i__1; ++i__) {
/* L120: */
h__[i__ + i__ * h_dim1] -= x;
}
s = (d__1 = h__[en + na * h_dim1], fabs(d__1)) + (d__2 = h__[na + enm2 *
h_dim1], fabs(d__2));
x = s * .75;
y = x;
w = s * -.4375 * s;
L130:
++its;
--itn;
/* .......... LOOK FOR TWO CONSECUTIVE SMALL */
/* SUB-DIAGONAL ELEMENTS. */
/* FOR M=EN-2 STEP -1 UNTIL L DO -- .......... */
i__1 = enm2;
for (mm = l; mm <= i__1; ++mm) {
m = enm2 + l - mm;
zz = h__[m + m * h_dim1];
r__ = x - zz;
s = y - zz;
p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) *
h_dim1];
q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s;
r__ = h__[m + 2 + (m + 1) * h_dim1];
s = fabs(p) + fabs(q) + fabs(r__);
p /= s;
q /= s;
r__ /= s;
if (m == l) {
goto L150;
}
tst1 = fabs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], fabs(d__1)) +
fabs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], fabs(d__2)));
tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], fabs(d__1)) * (fabs(q)
+ fabs(r__));
if (tst2 == tst1) {
goto L150;
}
/* L140: */
}
L150:
mp2 = m + 2;
i__1 = en;
for (i__ = mp2; i__ <= i__1; ++i__) {
h__[i__ + (i__ - 2) * h_dim1] = 0.;
if (i__ == mp2) {
goto L160;
}
h__[i__ + (i__ - 3) * h_dim1] = 0.;
L160:
;
}
/* .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND */
/* COLUMNS M TO EN .......... */
i__1 = na;
for (k = m; k <= i__1; ++k) {
notlas = k != na;
if (k == m) {
goto L170;
}
p = h__[k + (k - 1) * h_dim1];
q = h__[k + 1 + (k - 1) * h_dim1];
r__ = 0.;
if (notlas) {
r__ = h__[k + 2 + (k - 1) * h_dim1];
}
x = fabs(p) + fabs(q) + fabs(r__);
if (x == 0.) {
goto L260;
}
p /= x;
q /= x;
r__ /= x;
L170:
d__1 = sqrt(p * p + q * q + r__ * r__);
s = d_sign(d__1, p);
if (k == m) {
goto L180;
}
h__[k + (k - 1) * h_dim1] = -s * x;
goto L190;
L180:
if (l != m) {
h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1];
}
L190:
p += s;
x = p / s;
y = q / s;
zz = r__ / s;
q /= p;
r__ /= p;
if (notlas) {
goto L225;
}
/* .......... ROW MODIFICATION .......... */
i__2 = *n;
for (j = k; j <= i__2; ++j) {
p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1];
h__[k + j * h_dim1] -= p * x;
h__[k + 1 + j * h_dim1] -= p * y;
/* L200: */
}
/* Computing MIN */
i__2 = en, i__3 = k + 3;
j = min(i__2,i__3);
/* .......... COLUMN MODIFICATION .......... */
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__) {
p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1];
h__[i__ + k * h_dim1] -= p;
h__[i__ + (k + 1) * h_dim1] -= p * q;
/* L210: */
}
goto L255;
L225:
/* .......... ROW MODIFICATION .......... */
i__2 = *n;
for (j = k; j <= i__2; ++j) {
p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[
k + 2 + j * h_dim1];
h__[k + j * h_dim1] -= p * x;
h__[k + 1 + j * h_dim1] -= p * y;
h__[k + 2 + j * h_dim1] -= p * zz;
/* L230: */
}
/* Computing MIN */
i__2 = en, i__3 = k + 3;
j = min(i__2,i__3);
/* .......... COLUMN MODIFICATION .......... */
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__) {
p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] +
zz * h__[i__ + (k + 2) * h_dim1];
h__[i__ + k * h_dim1] -= p;
h__[i__ + (k + 1) * h_dim1] -= p * q;
h__[i__ + (k + 2) * h_dim1] -= p * r__;
/* L240: */
}
L255:
L260:
;
}
goto L70;
/* .......... ONE ROOT FOUND .......... */
L270:
wr[en] = x + t;
wi[en] = 0.;
en = na;
goto L60;
/* .......... TWO ROOTS FOUND .......... */
L280:
p = (y - x) / 2.;
q = p * p + w;
zz = sqrt((fabs(q)));
x += t;
if (q < 0.) {
goto L320;
}
/* .......... REAL PAIR .......... */
zz = p + d_sign(zz, p);
wr[na] = x + zz;
wr[en] = wr[na];
if (zz != 0.) {
wr[en] = x - w / zz;
}
wi[na] = 0.;
wi[en] = 0.;
goto L330;
/* .......... COMPLEX PAIR .......... */
L320:
wr[na] = x + p;
wr[en] = x + p;
wi[na] = zz;
wi[en] = -zz;
L330:
en = enm2;
goto L60;
/* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
/* CONVERGED AFTER 30*N ITERATIONS .......... */
L1000:
*ierr = en;
L1001:
return 0;
} /* hqr_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
hqr2(integer *nm, integer *n, integer *low, integer *igh, doublereal *h__, doublereal *wr, doublereal *wi, doublereal *z__, integer *ierr)
{
/* System generated locals */
integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3;
doublereal d__1, d__2, d__3, d__4;
/* Local variables */
doublereal norm;
integer i__, j, k, l, m;
doublereal p, q, r__, s, t, w, x, y;
integer na, ii, en, jj;
doublereal ra, sa;
integer ll, mm, nn;
doublereal vi, vr, zz;
logical notlas;
integer mp2, itn, its, enm2;
doublereal tst1, tst2;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2, */
/* NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
/* THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
/* OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD. THE */
/* EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND */
/* IF ELMHES AND ELTRAN OR ORTHES AND ORTRAN HAVE */
/* BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM */
/* AND TO ACCUMULATE THE SIMILARITY TRANSFORMATIONS. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRIX. */
/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
/* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */
/* SET LOW=1, IGH=N. */
/* H CONTAINS THE UPPER HESSENBERG MATRIX. */
/* Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY ELTRAN */
/* AFTER THE REDUCTION BY ELMHES, OR BY ORTRAN AFTER THE */
/* REDUCTION BY ORTHES, IF PERFORMED. IF THE EIGENVECTORS */
/* OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE */
/* IDENTITY MATRIX. */
/* ON OUTPUT */
/* H HAS BEEN DESTROYED. */
/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
/* RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES */
/* ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS */
/* OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE */
/* HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN */
/* ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
/* FOR INDICES IERR+1,...,N. */
/* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. */
/* IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z */
/* CONTAINS ITS EIGENVECTOR. IF THE I-TH EIGENVALUE IS COMPLEX
*/
/* WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH */
/* COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS */
/* EIGENVECTOR. THE EIGENVECTORS ARE UNNORMALIZED. IF AN */
/* ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND.
*/
/* IERR IS SET TO */
/* ZERO FOR NORMAL RETURN, */
/* J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
/* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
/* CALLS CDIV FOR COMPLEX DIVISION. */
/* 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;
--wi;
--wr;
h_dim1 = *nm;
h_offset = h_dim1 + 1;
h__ -= h_offset;
/* Function Body */
*ierr = 0;
norm = 0.;
k = 1;
/* .......... STORE ROOTS ISOLATED BY BALANC */
/* AND COMPUTE MATRIX NORM .......... */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = *n;
for (j = k; j <= i__2; ++j) {
/* L40: */
norm += (d__1 = h__[i__ + j * h_dim1], fabs(d__1));
}
k = i__;
if (i__ >= *low && i__ <= *igh) {
goto L50;
}
wr[i__] = h__[i__ + i__ * h_dim1];
wi[i__] = 0.;
L50:
;
}
en = *igh;
t = 0.;
itn = *n * 30;
/* .......... SEARCH FOR NEXT EIGENVALUES .......... */
L60:
if (en < *low) {
goto L340;
}
its = 0;
na = en - 1;
enm2 = na - 1;
/* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
/* FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
L70:
i__1 = en;
for (ll = *low; ll <= i__1; ++ll) {
l = en + *low - ll;
if (l == *low) {
goto L100;
}
s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], fabs(d__1)) + (d__2 = h__[l
+ l * h_dim1], fabs(d__2));
if (s == 0.) {
s = norm;
}
tst1 = s;
tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], fabs(d__1));
if (tst2 == tst1) {
goto L100;
}
/* L80: */
}
/* .......... FORM SHIFT .......... */
L100:
x = h__[en + en * h_dim1];
if (l == en) {
goto L270;
}
y = h__[na + na * h_dim1];
w = h__[en + na * h_dim1] * h__[na + en * h_dim1];
if (l == na) {
goto L280;
}
if (itn == 0) {
goto L1000;
}
if (its != 10 && its != 20) {
goto L130;
}
/* .......... FORM EXCEPTIONAL SHIFT .......... */
t += x;
i__1 = en;
for (i__ = *low; i__ <= i__1; ++i__) {
/* L120: */
h__[i__ + i__ * h_dim1] -= x;
}
s = (d__1 = h__[en + na * h_dim1], fabs(d__1)) + (d__2 = h__[na + enm2 *
h_dim1], fabs(d__2));
x = s * .75;
y = x;
w = s * -.4375 * s;
L130:
++its;
--itn;
/* .......... LOOK FOR TWO CONSECUTIVE SMALL */
/* SUB-DIAGONAL ELEMENTS. */
/* FOR M=EN-2 STEP -1 UNTIL L DO -- .......... */
i__1 = enm2;
for (mm = l; mm <= i__1; ++mm) {
m = enm2 + l - mm;
zz = h__[m + m * h_dim1];
r__ = x - zz;
s = y - zz;
p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) *
h_dim1];
q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s;
r__ = h__[m + 2 + (m + 1) * h_dim1];
s = fabs(p) + fabs(q) + fabs(r__);
p /= s;
q /= s;
r__ /= s;
if (m == l) {
goto L150;
}
tst1 = fabs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], fabs(d__1)) +
fabs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], fabs(d__2)));
tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], fabs(d__1)) * (fabs(q)
+ fabs(r__));
if (tst2 == tst1) {
goto L150;
}
/* L140: */
}
L150:
mp2 = m + 2;
i__1 = en;
for (i__ = mp2; i__ <= i__1; ++i__) {
h__[i__ + (i__ - 2) * h_dim1] = 0.;
if (i__ == mp2) {
goto L160;
}
h__[i__ + (i__ - 3) * h_dim1] = 0.;
L160:
;
}
/* .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND */
/* COLUMNS M TO EN .......... */
i__1 = na;
for (k = m; k <= i__1; ++k) {
notlas = k != na;
if (k == m) {
goto L170;
}
p = h__[k + (k - 1) * h_dim1];
q = h__[k + 1 + (k - 1) * h_dim1];
r__ = 0.;
if (notlas) {
r__ = h__[k + 2 + (k - 1) * h_dim1];
}
x = fabs(p) + fabs(q) + fabs(r__);
if (x == 0.) {
goto L260;
}
p /= x;
q /= x;
r__ /= x;
L170:
d__1 = sqrt(p * p + q * q + r__ * r__);
s = d_sign(d__1, p);
if (k == m) {
goto L180;
}
h__[k + (k - 1) * h_dim1] = -s * x;
goto L190;
L180:
if (l != m) {
h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1];
}
L190:
p += s;
x = p / s;
y = q / s;
zz = r__ / s;
q /= p;
r__ /= p;
if (notlas) {
goto L225;
}
/* .......... ROW MODIFICATION .......... */
i__2 = *n;
for (j = k; j <= i__2; ++j) {
p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1];
h__[k + j * h_dim1] -= p * x;
h__[k + 1 + j * h_dim1] -= p * y;
/* L200: */
}
/* Computing MIN */
i__2 = en, i__3 = k + 3;
j = min(i__2,i__3);
/* .......... COLUMN MODIFICATION .......... */
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__) {
p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1];
h__[i__ + k * h_dim1] -= p;
h__[i__ + (k + 1) * h_dim1] -= p * q;
/* L210: */
}
/* .......... ACCUMULATE TRANSFORMATIONS .......... */
i__2 = *igh;
for (i__ = *low; i__ <= i__2; ++i__) {
p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1];
z__[i__ + k * z_dim1] -= p;
z__[i__ + (k + 1) * z_dim1] -= p * q;
/* L220: */
}
goto L255;
L225:
/* .......... ROW MODIFICATION .......... */
i__2 = *n;
for (j = k; j <= i__2; ++j) {
p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[
k + 2 + j * h_dim1];
h__[k + j * h_dim1] -= p * x;
h__[k + 1 + j * h_dim1] -= p * y;
h__[k + 2 + j * h_dim1] -= p * zz;
/* L230: */
}
/* Computing MIN */
i__2 = en, i__3 = k + 3;
j = min(i__2,i__3);
/* .......... COLUMN MODIFICATION .......... */
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__) {
p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] +
zz * h__[i__ + (k + 2) * h_dim1];
h__[i__ + k * h_dim1] -= p;
h__[i__ + (k + 1) * h_dim1] -= p * q;
h__[i__ + (k + 2) * h_dim1] -= p * r__;
/* L240: */
}
/* .......... ACCUMULATE TRANSFORMATIONS .......... */
i__2 = *igh;
for (i__ = *low; i__ <= i__2; ++i__) {
p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1] +
zz * z__[i__ + (k + 2) * z_dim1];
z__[i__ + k * z_dim1] -= p;
z__[i__ + (k + 1) * z_dim1] -= p * q;
z__[i__ + (k + 2) * z_dim1] -= p * r__;
/* L250: */
}
L255:
L260:
;
}
goto L70;
/* .......... ONE ROOT FOUND .......... */
L270:
h__[en + en * h_dim1] = x + t;
wr[en] = h__[en + en * h_dim1];
wi[en] = 0.;
en = na;
goto L60;
/* .......... TWO ROOTS FOUND .......... */
L280:
p = (y - x) / 2.;
q = p * p + w;
zz = sqrt((fabs(q)));
h__[en + en * h_dim1] = x + t;
x = h__[en + en * h_dim1];
h__[na + na * h_dim1] = y + t;
if (q < 0.) {
goto L320;
}
/* .......... REAL PAIR .......... */
zz = p + d_sign(zz, p);
wr[na] = x + zz;
wr[en] = wr[na];
if (zz != 0.) {
wr[en] = x - w / zz;
}
wi[na] = 0.;
wi[en] = 0.;
x = h__[en + na * h_dim1];
s = fabs(x) + fabs(zz);
p = x / s;
q = zz / s;
r__ = sqrt(p * p + q * q);
p /= r__;
q /= r__;
/* .......... ROW MODIFICATION .......... */
i__1 = *n;
for (j = na; j <= i__1; ++j) {
zz = h__[na + j * h_dim1];
h__[na + j * h_dim1] = q * zz + p * h__[en + j * h_dim1];
h__[en + j * h_dim1] = q * h__[en + j * h_dim1] - p * zz;
/* L290: */
}
/* .......... COLUMN MODIFICATION .......... */
i__1 = en;
for (i__ = 1; i__ <= i__1; ++i__) {
zz = h__[i__ + na * h_dim1];
h__[i__ + na * h_dim1] = q * zz + p * h__[i__ + en * h_dim1];
h__[i__ + en * h_dim1] = q * h__[i__ + en * h_dim1] - p * zz;
/* L300: */
}
/* .......... ACCUMULATE TRANSFORMATIONS .......... */
i__1 = *igh;
for (i__ = *low; i__ <= i__1; ++i__) {
zz = z__[i__ + na * z_dim1];
z__[i__ + na * z_dim1] = q * zz + p * z__[i__ + en * z_dim1];
z__[i__ + en * z_dim1] = q * z__[i__ + en * z_dim1] - p * zz;
/* L310: */
}
goto L330;
/* .......... COMPLEX PAIR .......... */
L320:
wr[na] = x + p;
wr[en] = x + p;
wi[na] = zz;
wi[en] = -zz;
L330:
en = enm2;
goto L60;
/* .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND */
/* VECTORS OF UPPER TRIANGULAR FORM .......... */
L340:
if (norm == 0.) {
goto L1001;
}
/* .......... FOR EN=N STEP -1 UNTIL 1 DO -- .......... */
i__1 = *n;
for (nn = 1; nn <= i__1; ++nn) {
en = *n + 1 - nn;
p = wr[en];
q = wi[en];
na = en - 1;
if (q < 0.) {
goto L710;
} else if (q == 0) {
goto L600;
} else {
goto L800;
}
/* .......... REAL VECTOR .......... */
L600:
m = en;
h__[en + en * h_dim1] = 1.;
if (na == 0) {
goto L800;
}
/* .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */
i__2 = na;
for (ii = 1; ii <= i__2; ++ii) {
i__ = en - ii;
w = h__[i__ + i__ * h_dim1] - p;
r__ = 0.;
i__3 = en;
for (j = m; j <= i__3; ++j) {
/* L610: */
r__ += h__[i__ + j * h_dim1] * h__[j + en * h_dim1];
}
if (wi[i__] >= 0.) {
goto L630;
}
zz = w;
s = r__;
goto L700;
L630:
m = i__;
if (wi[i__] != 0.) {
goto L640;
}
t = w;
if (t != 0.) {
goto L635;
}
tst1 = norm;
t = tst1;
L632:
t *= .01;
tst2 = norm + t;
if (tst2 > tst1) {
goto L632;
}
L635:
h__[i__ + en * h_dim1] = -r__ / t;
goto L680;
/* .......... SOLVE REAL EQUATIONS .......... */
L640:
x = h__[i__ + (i__ + 1) * h_dim1];
y = h__[i__ + 1 + i__ * h_dim1];
q = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__];
t = (x * s - zz * r__) / q;
h__[i__ + en * h_dim1] = t;
if (fabs(x) <= fabs(zz)) {
goto L650;
}
h__[i__ + 1 + en * h_dim1] = (-r__ - w * t) / x;
goto L680;
L650:
h__[i__ + 1 + en * h_dim1] = (-s - y * t) / zz;
/* .......... OVERFLOW CONTROL .......... */
L680:
t = (d__1 = h__[i__ + en * h_dim1], fabs(d__1));
if (t == 0.) {
goto L700;
}
tst1 = t;
tst2 = tst1 + 1. / tst1;
if (tst2 > tst1) {
goto L700;
}
i__3 = en;
for (j = i__; j <= i__3; ++j) {
h__[j + en * h_dim1] /= t;
/* L690: */
}
L700:
;
}
/* .......... END REAL VECTOR .......... */
goto L800;
/* .......... COMPLEX VECTOR .......... */
L710:
m = na;
/* .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT */
/* EIGENVECTOR MATRIX IS TRIANGULAR .......... */
if ((d__1 = h__[en + na * h_dim1], fabs(d__1)) <= (d__2 = h__[na + en *
h_dim1], fabs(d__2))) {
goto L720;
}
h__[na + na * h_dim1] = q / h__[en + na * h_dim1];
h__[na + en * h_dim1] = -(h__[en + en * h_dim1] - p) / h__[en + na *
h_dim1];
goto L730;
L720:
d__1 = -h__[na + en * h_dim1];
d__2 = h__[na + na * h_dim1] - p;
{
doublereal c_b81 = 0.;
cdiv(&c_b81, &d__1, &d__2, &q, &h__[na + na * h_dim1], &h__[na + en *
h_dim1]);
}
L730:
h__[en + na * h_dim1] = 0.;
h__[en + en * h_dim1] = 1.;
enm2 = na - 1;
if (enm2 == 0) {
goto L800;
}
/* .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- .......... */
i__2 = enm2;
for (ii = 1; ii <= i__2; ++ii) {
i__ = na - ii;
w = h__[i__ + i__ * h_dim1] - p;
ra = 0.;
sa = 0.;
i__3 = en;
for (j = m; j <= i__3; ++j) {
ra += h__[i__ + j * h_dim1] * h__[j + na * h_dim1];
sa += h__[i__ + j * h_dim1] * h__[j + en * h_dim1];
/* L760: */
}
if (wi[i__] >= 0.) {
goto L770;
}
zz = w;
r__ = ra;
s = sa;
goto L795;
L770:
m = i__;
if (wi[i__] != 0.) {
goto L780;
}
d__1 = -ra;
d__2 = -sa;
cdiv(&d__1, &d__2, &w, &q, &h__[i__ + na * h_dim1], &h__[i__ +
en * h_dim1]);
goto L790;
/* .......... SOLVE COMPLEX EQUATIONS .......... */
L780:
x = h__[i__ + (i__ + 1) * h_dim1];
y = h__[i__ + 1 + i__ * h_dim1];
vr = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__] - q * q;
vi = (wr[i__] - p) * 2. * q;
if (vr != 0. || vi != 0.) {
goto L784;
}
tst1 = norm * (fabs(w) + fabs(q) + fabs(x) + fabs(y) + fabs(zz));
vr = tst1;
L783:
vr *= .01;
tst2 = tst1 + vr;
if (tst2 > tst1) {
goto L783;
}
L784:
d__1 = x * r__ - zz * ra + q * sa;
d__2 = x * s - zz * sa - q * ra;
cdiv(&d__1, &d__2, &vr, &vi, &h__[i__ + na * h_dim1], &h__[i__ +
en * h_dim1]);
if (fabs(x) <= fabs(zz) + fabs(q)) {
goto L785;
}
h__[i__ + 1 + na * h_dim1] = (-ra - w * h__[i__ + na * h_dim1] +
q * h__[i__ + en * h_dim1]) / x;
h__[i__ + 1 + en * h_dim1] = (-sa - w * h__[i__ + en * h_dim1] -
q * h__[i__ + na * h_dim1]) / x;
goto L790;
L785:
d__1 = -r__ - y * h__[i__ + na * h_dim1];
d__2 = -s - y * h__[i__ + en * h_dim1];
cdiv(&d__1, &d__2, &zz, &q, &h__[i__ + 1 + na * h_dim1], &h__[
i__ + 1 + en * h_dim1]);
/* .......... OVERFLOW CONTROL .......... */
L790:
/* Computing MAX */
d__3 = (d__1 = h__[i__ + na * h_dim1], fabs(d__1)), d__4 = (d__2 =
h__[i__ + en * h_dim1], fabs(d__2));
t = max(d__3,d__4);
if (t == 0.) {
goto L795;
}
tst1 = t;
tst2 = tst1 + 1. / tst1;
if (tst2 > tst1) {
goto L795;
}
i__3 = en;
for (j = i__; j <= i__3; ++j) {
h__[j + na * h_dim1] /= t;
h__[j + en * h_dim1] /= t;
/* L792: */
}
L795:
;
}
/* .......... END COMPLEX VECTOR .......... */
L800:
;
}
/* .......... END BACK SUBSTITUTION. */
/* VECTORS OF ISOLATED ROOTS .......... */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
if (i__ >= *low && i__ <= *igh) {
goto L840;
}
i__2 = *n;
for (j = i__; j <= i__2; ++j) {
/* L820: */
z__[i__ + j * z_dim1] = h__[i__ + j * h_dim1];
}
L840:
;
}
/* .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */
/* VECTORS OF ORIGINAL FULL MATRIX. */
/* FOR J=N STEP -1 UNTIL LOW DO -- .......... */
i__1 = *n;
for (jj = *low; jj <= i__1; ++jj) {
j = *n + *low - jj;
m = min(j,*igh);
i__2 = *igh;
for (i__ = *low; i__ <= i__2; ++i__) {
zz = 0.;
i__3 = m;
for (k = *low; k <= i__3; ++k) {
/* L860: */
zz += z__[i__ + k * z_dim1] * h__[k + j * h_dim1];
}
z__[i__ + j * z_dim1] = zz;
/* L880: */
}
}
goto L1001;
/* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
/* CONVERGED AFTER 30*N ITERATIONS .......... */
L1000:
*ierr = en;
L1001:
return 0;
} /* hqr2_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
cdiv(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci)
{
/* System generated locals */
doublereal d__1, d__2;
/* Local variables */
doublereal s, ais, bis, ars, brs;
/* COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI) */
s = fabs(*br) + fabs(*bi);
ars = *ar / s;
ais = *ai / s;
brs = *br / s;
bis = *bi / s;
/* Computing 2nd power */
d__1 = brs;
/* Computing 2nd power */
d__2 = bis;
s = d__1 * d__1 + d__2 * d__2;
*cr = (ars * brs + ais * bis) / s;
*ci = (ais * brs - ars * bis) / s;
return 0;
} /* cdiv_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
balanc(integer *nm, integer *n, doublereal *a, integer *low, integer *igh, doublereal *scale)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer iexc;
doublereal c__, f, g;
integer i__, j, k, l, m;
doublereal r__, s, radix, b2;
integer jj;
logical noconv;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE, */
/* NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
/* THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES */
/* EIGENVALUES WHENEVER POSSIBLE. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRIX. */
/* A CONTAINS THE INPUT MATRIX TO BE BALANCED. */
/* ON OUTPUT */
/* A CONTAINS THE BALANCED MATRIX. */
/* LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J) */
/* IS EQUAL TO ZERO IF */
/* (1) I IS GREATER THAN J AND */
/* (2) J=1,...,LOW-1 OR I=IGH+1,...,N. */
/* SCALE CONTAINS INFORMATION DETERMINING THE */
/* PERMUTATIONS AND SCALING FACTORS USED. */
/* SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH */
/* HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED */
/* WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS */
/* OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN */
/* SCALE(J) = P(J), FOR J = 1,...,LOW-1 */
/* = D(J,J), J = LOW,...,IGH */
/* = P(J) J = IGH+1,...,N. */
/* THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, */
/* THEN 1 TO LOW-1. */
/* NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. */
/* THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN */
/* BALANC IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS */
/* K,L HAVE BEEN REVERSED.) */
/* 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 */
--scale;
a_dim1 = *nm;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
radix = 16.;
b2 = radix * radix;
k = 1;
l = *n;
goto L100;
/* .......... IN-LINE PROCEDURE FOR ROW AND */
/* COLUMN EXCHANGE .......... */
L20:
scale[m] = (doublereal) j;
if (j == m) {
goto L50;
}
i__1 = l;
for (i__ = 1; i__ <= i__1; ++i__) {
f = a[i__ + j * a_dim1];
a[i__ + j * a_dim1] = a[i__ + m * a_dim1];
a[i__ + m * a_dim1] = f;
/* L30: */
}
i__1 = *n;
for (i__ = k; i__ <= i__1; ++i__) {
f = a[j + i__ * a_dim1];
a[j + i__ * a_dim1] = a[m + i__ * a_dim1];
a[m + i__ * a_dim1] = f;
/* L40: */
}
L50:
switch ((int)iexc) {
case 1: goto L80;
case 2: goto L130;
}
/* .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE */
/* AND PUSH THEM DOWN .......... */
L80:
if (l == 1) {
goto L280;
}
--l;
/* .......... FOR J=L STEP -1 UNTIL 1 DO -- .......... */
L100:
i__1 = l;
for (jj = 1; jj <= i__1; ++jj) {
j = l + 1 - jj;
i__2 = l;
for (i__ = 1; i__ <= i__2; ++i__) {
if (i__ == j) {
goto L110;
}
if (a[j + i__ * a_dim1] != 0.) {
goto L120;
}
L110:
;
}
m = l;
iexc = 1;
goto L20;
L120:
;
}
goto L140;
/* .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE */
/* AND PUSH THEM LEFT .......... */
L130:
++k;
L140:
i__1 = l;
for (j = k; j <= i__1; ++j) {
i__2 = l;
for (i__ = k; i__ <= i__2; ++i__) {
if (i__ == j) {
goto L150;
}
if (a[i__ + j * a_dim1] != 0.) {
goto L170;
}
L150:
;
}
m = k;
iexc = 2;
goto L20;
L170:
;
}
/* .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L .......... */
i__1 = l;
for (i__ = k; i__ <= i__1; ++i__) {
/* L180: */
scale[i__] = 1.;
}
/* .......... ITERATIVE LOOP FOR NORM REDUCTION .......... */
L190:
noconv = FALSE_;
i__1 = l;
for (i__ = k; i__ <= i__1; ++i__) {
c__ = 0.;
r__ = 0.;
i__2 = l;
for (j = k; j <= i__2; ++j) {
if (j == i__) {
goto L200;
}
c__ += (d__1 = a[j + i__ * a_dim1], fabs(d__1));
r__ += (d__1 = a[i__ + j * a_dim1], fabs(d__1));
L200:
;
}
/* .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW .........
. */
if (c__ == 0. || r__ == 0.) {
goto L270;
}
g = r__ / radix;
f = 1.;
s = c__ + r__;
L210:
if (c__ >= g) {
goto L220;
}
f *= radix;
c__ *= b2;
goto L210;
L220:
g = r__ * radix;
L230:
if (c__ < g) {
goto L240;
}
f /= radix;
c__ /= b2;
goto L230;
/* .......... NOW BALANCE .......... */
L240:
if ((c__ + r__) / f >= s * .95) {
goto L270;
}
g = 1. / f;
scale[i__] *= f;
noconv = TRUE_;
i__2 = *n;
for (j = k; j <= i__2; ++j) {
/* L250: */
a[i__ + j * a_dim1] *= g;
}
i__2 = l;
for (j = 1; j <= i__2; ++j) {
/* L260: */
a[j + i__ * a_dim1] *= f;
}
L270:
;
}
if (noconv) {
goto L190;
}
L280:
*low = k;
*igh = l;
return 0;
} /* balanc_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
balbak(integer *nm, integer *n, integer *low, integer *igh, doublereal *scale, integer *m, doublereal *z__)
{
/* System generated locals */
integer z_dim1, z_offset, i__1, i__2;
/* Local variables */
integer i__, j, k;
doublereal s;
integer ii;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALBAK, */
/* NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
/* THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL */
/* MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
/* BALANCED MATRIX DETERMINED BY BALANC. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRIX. */
/* LOW AND IGH ARE INTEGERS DETERMINED BY BALANC. */
/* SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS */
/* AND SCALING FACTORS USED BY BALANC. */
/* M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. */
/* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- */
/* VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. */
/* ON OUTPUT */
/* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE */
/* TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. */
/* 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 */
--scale;
z_dim1 = *nm;
z_offset = z_dim1 + 1;
z__ -= z_offset;
/* Function Body */
if (*m == 0) {
goto L200;
}
if (*igh == *low) {
goto L120;
}
i__1 = *igh;
for (i__ = *low; i__ <= i__1; ++i__) {
s = scale[i__];
/* .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED */
/* IF THE FOREGOING STATEMENT IS REPLACED BY */
/* S=1.0D0/SCALE(I). .......... */
i__2 = *m;
for (j = 1; j <= i__2; ++j) {
/* L100: */
z__[i__ + j * z_dim1] *= s;
}
/* L110: */
}
/* ......... FOR I=LOW-1 STEP -1 UNTIL 1, */
/* IGH+1 STEP 1 UNTIL N DO -- .......... */
L120:
i__1 = *n;
for (ii = 1; ii <= i__1; ++ii) {
i__ = ii;
if (i__ >= *low && i__ <= *igh) {
goto L140;
}
if (i__ < *low) {
i__ = *low - ii;
}
k = (integer) scale[i__];
if (k == i__) {
goto L140;
}
i__2 = *m;
for (j = 1; j <= i__2; ++j) {
s = z__[i__ + j * z_dim1];
z__[i__ + j * z_dim1] = z__[k + j * z_dim1];
z__[k + j * z_dim1] = s;
/* L130: */
}
L140:
;
}
L200:
return 0;
} /* balbak_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
elmhes(integer *nm, integer *n, integer *low, integer *igh, doublereal *a, integer *int__)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3;
doublereal d__1;
/* Local variables */
integer i__, j, m;
doublereal x, y;
integer la, mm1, kp1, mp1;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES, */
/* NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
/* GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE */
/* REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS */
/* LOW THROUGH IGH TO UPPER HESSENBERG FORM BY */
/* STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRIX. */
/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
/* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */
/* SET LOW=1, IGH=N. */
/* A CONTAINS THE INPUT MATRIX. */
/* ON OUTPUT */
/* A CONTAINS THE HESSENBERG MATRIX. THE MULTIPLIERS */
/* WHICH WERE USED IN THE REDUCTION ARE STORED IN THE */
/* REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. */
/* INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
/* INTERCHANGED IN THE REDUCTION. */
/* ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
/* 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 */
a_dim1 = *nm;
a_offset = a_dim1 + 1;
a -= a_offset;
--int__;
/* Function Body */
la = *igh - 1;
kp1 = *low + 1;
if (la < kp1) {
goto L200;
}
i__1 = la;
for (m = kp1; m <= i__1; ++m) {
mm1 = m - 1;
x = 0.;
i__ = m;
i__2 = *igh;
for (j = m; j <= i__2; ++j) {
if ((d__1 = a[j + mm1 * a_dim1], fabs(d__1)) <= fabs(x)) {
goto L100;
}
x = a[j + mm1 * a_dim1];
i__ = j;
L100:
;
}
int__[m] = i__;
if (i__ == m) {
goto L130;
}
/* .......... INTERCHANGE ROWS AND COLUMNS OF A .......... */
i__2 = *n;
for (j = mm1; j <= i__2; ++j) {
y = a[i__ + j * a_dim1];
a[i__ + j * a_dim1] = a[m + j * a_dim1];
a[m + j * a_dim1] = y;
/* L110: */
}
i__2 = *igh;
for (j = 1; j <= i__2; ++j) {
y = a[j + i__ * a_dim1];
a[j + i__ * a_dim1] = a[j + m * a_dim1];
a[j + m * a_dim1] = y;
/* L120: */
}
/* .......... END INTERCHANGE .......... */
L130:
if (x == 0.) {
goto L180;
}
mp1 = m + 1;
i__2 = *igh;
for (i__ = mp1; i__ <= i__2; ++i__) {
y = a[i__ + mm1 * a_dim1];
if (y == 0.) {
goto L160;
}
y /= x;
a[i__ + mm1 * a_dim1] = y;
i__3 = *n;
for (j = m; j <= i__3; ++j) {
/* L140: */
a[i__ + j * a_dim1] -= y * a[m + j * a_dim1];
}
i__3 = *igh;
for (j = 1; j <= i__3; ++j) {
/* L150: */
a[j + m * a_dim1] += y * a[j + i__ * a_dim1];
}
L160:
;
}
L180:
;
}
L200:
return 0;
} /* elmhes_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
eltran(integer *nm, integer *n, integer *low, integer *igh, doublereal *a, integer *int__, doublereal *z__)
{
/* System generated locals */
integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
/* Local variables */
integer i__, j, kl, mm, mp, mp1;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMTRANS,
*/
/* NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
/* THIS SUBROUTINE ACCUMULATES THE STABILIZED ELEMENTARY */
/* SIMILARITY TRANSFORMATIONS USED IN THE REDUCTION OF A */
/* REAL GENERAL MATRIX TO UPPER HESSENBERG FORM BY ELMHES. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRIX. */
/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
/* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */
/* SET LOW=1, IGH=N. */
/* A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE */
/* REDUCTION BY ELMHES IN ITS LOWER TRIANGLE */
/* BELOW THE SUBDIAGONAL. */
/* INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
/* INTERCHANGED IN THE REDUCTION BY ELMHES. */
/* ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
/* ON OUTPUT */
/* Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE */
/* REDUCTION BY ELMHES. */
/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
*/
/* THIS VERSION DATED AUGUST 1983. */
/* ------------------------------------------------------------------
*/
/* .......... INITIALIZE Z TO IDENTITY MATRIX .......... */
/* Parameter adjustments */
z_dim1 = *nm;
z_offset = z_dim1 + 1;
z__ -= z_offset;
--int__;
a_dim1 = *nm;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
/* L60: */
z__[i__ + j * z_dim1] = 0.;
}
z__[j + j * z_dim1] = 1.;
/* L80: */
}
kl = *igh - *low - 1;
if (kl < 1) {
goto L200;
}
/* .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
i__1 = kl;
for (mm = 1; mm <= i__1; ++mm) {
mp = *igh - mm;
mp1 = mp + 1;
i__2 = *igh;
for (i__ = mp1; i__ <= i__2; ++i__) {
/* L100: */
z__[i__ + mp * z_dim1] = a[i__ + (mp - 1) * a_dim1];
}
i__ = int__[mp];
if (i__ == mp) {
goto L140;
}
i__2 = *igh;
for (j = mp; j <= i__2; ++j) {
z__[mp + j * z_dim1] = z__[i__ + j * z_dim1];
z__[i__ + j * z_dim1] = 0.;
/* L130: */
}
z__[i__ + mp * z_dim1] = 1.;
L140:
;
}
L200:
return 0;
} /* eltran_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* EISPACK routines needed in the computation of Floquet multipliers */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
qzhes(integer nm, integer n, doublereal *a, doublereal *b, logical matz, 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;
/* Local variables */
integer i__, j, k, l;
doublereal r__, s, t;
integer l1;
doublereal u1, u2, v1, v2;
integer lb, nk1, nm1, nm2;
doublereal rho;
/* THIS SUBROUTINE IS THE FIRST 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 GENERAL MATRICES AND */
/* REDUCES ONE OF THEM TO UPPER HESSENBERG FORM AND THE OTHER */
/* TO UPPER TRIANGULAR FORM USING ORTHOGONAL TRANSFORMATIONS. */
/* IT IS USUALLY FOLLOWED BY QZIT, QZVAL AND, POSSIBLY, QZVEC. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRICES. */
/* A CONTAINS A REAL GENERAL MATRIX. */
/* B CONTAINS A REAL GENERAL MATRIX. */
/* 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. */
/* ON OUTPUT */
/* A HAS BEEN REDUCED TO UPPER HESSENBERG FORM. THE ELEMENTS */
/* BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO. */
/* B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS */
/* BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO. */
/* Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS IF */
/* MATZ HAS BEEN SET TO .TRUE. OTHERWISE, Z IS NOT REFERENCED.
*/
/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
*/
/* THIS VERSION DATED AUGUST 1983. */
/* ------------------------------------------------------------------
*/
/* .......... INITIALIZE Z .......... */
/* Parameter adjustments */
z_dim1 = nm;
z_offset = z_dim1 + 1;
z__ -= z_offset;
b_dim1 = nm;
b_offset = b_dim1 + 1;
b -= b_offset;
a_dim1 = nm;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
if (! (matz)) {
goto L10;
}
i__1 = n;
for (j = 1; j <= i__1; ++j) {
i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
z__[i__ + j * z_dim1] = 0.;
/* L2: */
}
z__[j + j * z_dim1] = 1.;
/* L3: */
}
/* .......... REDUCE B TO UPPER TRIANGULAR FORM .......... */
L10:
if (n <= 1) {
goto L170;
}
nm1 = n - 1;
i__1 = nm1;
for (l = 1; l <= i__1; ++l) {
l1 = l + 1;
s = 0.;
i__2 = n;
for (i__ = l1; i__ <= i__2; ++i__) {
s += (d__1 = b[i__ + l * b_dim1], fabs(d__1));
/* L20: */
}
if (s == 0.) {
goto L100;
}
s += (d__1 = b[l + l * b_dim1], fabs(d__1));
r__ = 0.;
i__2 = n;
for (i__ = l; i__ <= i__2; ++i__) {
b[i__ + l * b_dim1] /= s;
/* Computing 2nd power */
d__1 = b[i__ + l * b_dim1];
r__ += d__1 * d__1;
/* L25: */
}
d__1 = sqrt(r__);
r__ = d_sign(d__1, b[l + l * b_dim1]);
b[l + l * b_dim1] += r__;
rho = r__ * b[l + l * b_dim1];
i__2 = n;
for (j = l1; j <= i__2; ++j) {
t = 0.;
i__3 = n;
for (i__ = l; i__ <= i__3; ++i__) {
t += b[i__ + l * b_dim1] * b[i__ + j * b_dim1];
/* L30: */
}
t = -t / rho;
i__3 = n;
for (i__ = l; i__ <= i__3; ++i__) {
b[i__ + j * b_dim1] += t * b[i__ + l * b_dim1];
/* L40: */
}
/* L50: */
}
i__2 = n;
for (j = 1; j <= i__2; ++j) {
t = 0.;
i__3 = n;
for (i__ = l; i__ <= i__3; ++i__) {
t += b[i__ + l * b_dim1] * a[i__ + j * a_dim1];
/* L60: */
}
t = -t / rho;
i__3 = n;
for (i__ = l; i__ <= i__3; ++i__) {
a[i__ + j * a_dim1] += t * b[i__ + l * b_dim1];
/* L70: */
}
/* L80: */
}
b[l + l * b_dim1] = -s * r__;
i__2 = n;
for (i__ = l1; i__ <= i__2; ++i__) {
b[i__ + l * b_dim1] = 0.;
/* L90: */
}
L100:
;
}
/* .......... REDUCE A TO UPPER HESSENBERG FORM, WHILE */
/* KEEPING B TRIANGULAR .......... */
if (n == 2) {
goto L170;
}
nm2 = n - 2;
i__1 = nm2;
for (k = 1; k <= i__1; ++k) {
nk1 = nm1 - k;
/* .......... FOR L=N-1 STEP -1 UNTIL K+1 DO -- .......... */
i__2 = nk1;
for (lb = 1; lb <= i__2; ++lb) {
l = n - lb;
l1 = l + 1;
/* .......... ZERO A(L+1,K) .......... */
s = (d__1 = a[l + k * a_dim1], fabs(d__1)) + (d__2 = a[l1 + k *
a_dim1], fabs(d__2));
if (s == 0.) {
goto L150;
}
u1 = a[l + k * a_dim1] / s;
u2 = a[l1 + k * a_dim1] / 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__3 = n;
for (j = k; j <= i__3; ++j) {
t = a[l + j * a_dim1] + u2 * a[l1 + j * a_dim1];
a[l + j * a_dim1] += t * v1;
a[l1 + j * a_dim1] += t * v2;
/* L110: */
}
a[l1 + k * a_dim1] = 0.;
i__3 = n;
for (j = l; j <= i__3; ++j) {
t = b[l + j * b_dim1] + u2 * b[l1 + j * b_dim1];
b[l + j * b_dim1] += t * v1;
b[l1 + j * b_dim1] += t * v2;
/* L120: */
}
/* .......... ZERO B(L+1,L) .......... */
s = (d__1 = b[l1 + l1 * b_dim1], fabs(d__1)) + (d__2 = b[l1 + l *
b_dim1], fabs(d__2));
if (s == 0.) {
goto L150;
}
u1 = b[l1 + l1 * b_dim1] / s;
u2 = b[l1 + l * b_dim1] / 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__3 = l1;
for (i__ = 1; i__ <= i__3; ++i__) {
t = b[i__ + l1 * b_dim1] + u2 * b[i__ + l * b_dim1];
b[i__ + l1 * b_dim1] += t * v1;
b[i__ + l * b_dim1] += t * v2;
/* L130: */
}
b[l1 + l * b_dim1] = 0.;
i__3 = n;
for (i__ = 1; i__ <= i__3; ++i__) {
t = a[i__ + l1 * a_dim1] + u2 * a[i__ + l * a_dim1];
a[i__ + l1 * a_dim1] += t * v1;
a[i__ + l * a_dim1] += t * v2;
/* L140: */
}
if (! (matz)) {
goto L150;
}
i__3 = n;
for (i__ = 1; i__ <= i__3; ++i__) {
t = z__[i__ + l1 * z_dim1] + u2 * z__[i__ + l * z_dim1];
z__[i__ + l1 * z_dim1] += t * v1;
z__[i__ + l * z_dim1] += t * v2;
/* L145: */
}
L150:
;
}
/* L160: */
}
L170:
return 0;
} /* qzhes_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
qzit(integer nm, integer n, doublereal *a, doublereal *b, doublereal eps1, logical matz, doublereal *z__, integer *ierr)
{
/* 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, d__3;
/* Local variables */
doublereal epsa, epsb;
integer i__, j, k, l;
doublereal r__, s, t, anorm, bnorm;
integer enorn;
doublereal a1, a2, a3;
integer k1, k2, l1;
doublereal u1, u2, u3, v1, v2, v3, a11, a12, a21, a22, a33, a34,
a43, a44, b11, b12, b22, b33;
integer na, ld;
doublereal b34, b44;
integer en;
doublereal ep;
integer ll;
doublereal sh;
logical notlas;
integer km1, lm1;
doublereal ani, bni;
integer ish, itn, its, enm2, lor1;
/* THIS SUBROUTINE IS THE SECOND STEP OF THE QZ ALGORITHM */
/* FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */
/* SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART, */
/* AS MODIFIED IN TECHNICAL NOTE NASA TN D-7305(1973) BY WARD. */
/* THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM */
/* IN UPPER HESSENBERG FORM AND THE OTHER IN UPPER TRIANGULAR FORM. */
/* IT REDUCES THE HESSENBERG MATRIX TO QUASI-TRIANGULAR FORM USING */
/* ORTHOGONAL TRANSFORMATIONS WHILE MAINTAINING THE TRIANGULAR FORM */
/* OF THE OTHER MATRIX. IT IS USUALLY PRECEDED BY QZHES AND */
/* FOLLOWED BY QZVAL AND, POSSIBLY, QZVEC. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRICES. */
/* A CONTAINS A REAL UPPER HESSENBERG MATRIX. */
/* B CONTAINS A REAL UPPER TRIANGULAR MATRIX. */
/* EPS1 IS A TOLERANCE USED TO DETERMINE NEGLIGIBLE ELEMENTS. */
/* EPS1 = 0.0 (OR NEGATIVE) MAY BE INPUT, IN WHICH CASE AN */
/* ELEMENT WILL BE NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF */
/* ERROR TIMES THE NORM OF ITS MATRIX. IF THE INPUT EPS1 IS */
/* POSITIVE, THEN AN ELEMENT WILL BE CONSIDERED NEGLIGIBLE */
/* IF IT IS LESS THAN EPS1 TIMES THE NORM OF ITS MATRIX. A */
/* POSITIVE VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION, */
/* BUT LESS ACCURATE RESULTS. */
/* 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 REDUCTION */
/* BY QZHES, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. */
/* IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. */
/* ON OUTPUT */
/* A HAS BEEN REDUCED TO QUASI-TRIANGULAR FORM. THE ELEMENTS */
/* BELOW THE FIRST SUBDIAGONAL ARE STILL ZERO AND NO TWO */
/* CONSECUTIVE SUBDIAGONAL ELEMENTS ARE NONZERO. */
/* B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS */
/* HAVE BEEN ALTERED. THE LOCATION B(N,1) IS USED TO STORE */
/* EPS1 TIMES THE NORM OF B FOR LATER USE BY QZVAL AND QZVEC.
*/
/* Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS */
/* (FOR BOTH STEPS) IF MATZ HAS BEEN SET TO .TRUE.. */
/* IERR IS SET TO */
/* ZERO FOR NORMAL RETURN, */
/* J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
/* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
/* 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;
b_dim1 = nm;
b_offset = b_dim1 + 1;
b -= b_offset;
a_dim1 = nm;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
*ierr = 0;
/* .......... COMPUTE EPSA,EPSB .......... */
anorm = 0.;
bnorm = 0.;
i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
ani = 0.;
if (i__ != 1) {
ani = (d__1 = a[i__ + (i__ - 1) * a_dim1], fabs(d__1));
}
bni = 0.;
i__2 = n;
for (j = i__; j <= i__2; ++j) {
ani += (d__1 = a[i__ + j * a_dim1], fabs(d__1));
bni += (d__1 = b[i__ + j * b_dim1], fabs(d__1));
/* L20: */
}
if (ani > anorm) {
anorm = ani;
}
if (bni > bnorm) {
bnorm = bni;
}
/* L30: */
}
if (anorm == 0.) {
anorm = 1.;
}
if (bnorm == 0.) {
bnorm = 1.;
}
ep = eps1;
if (ep > 0.) {
goto L50;
}
/* .......... USE ROUNDOFF LEVEL IF EPS1 IS ZERO .......... */
ep = epslon(1.0);
L50:
epsa = ep * anorm;
epsb = ep * bnorm;
/* .......... REDUCE A TO QUASI-TRIANGULAR FORM, WHILE */
/* KEEPING B TRIANGULAR .......... */
lor1 = 1;
enorn = n;
en = n;
itn = n * 30;
/* .......... BEGIN QZ STEP .......... */
L60:
if (en <= 2) {
goto L1001;
}
if (! (matz)) {
enorn = en;
}
its = 0;
na = en - 1;
enm2 = na - 1;
L70:
ish = 2;
/* .......... CHECK FOR CONVERGENCE OR REDUCIBILITY. */
/* FOR L=EN STEP -1 UNTIL 1 DO -- .......... */
i__1 = en;
for (ll = 1; ll <= i__1; ++ll) {
lm1 = en - ll;
l = lm1 + 1;
if (l == 1) {
goto L95;
}
if ((d__1 = a[l + lm1 * a_dim1], fabs(d__1)) <= epsa) {
goto L90;
}
/* L80: */
}
L90:
a[l + lm1 * a_dim1] = 0.;
if (l < na) {
goto L95;
}
/* .......... 1-BY-1 OR 2-BY-2 BLOCK ISOLATED .......... */
en = lm1;
goto L60;
/* .......... CHECK FOR SMALL TOP OF B .......... */
L95:
ld = l;
L100:
l1 = l + 1;
b11 = b[l + l * b_dim1];
if (fabs(b11) > epsb) {
goto L120;
}
b[l + l * b_dim1] = 0.;
s = (d__1 = a[l + l * a_dim1], fabs(d__1)) + (d__2 = a[l1 + l * a_dim1],
fabs(d__2));
u1 = a[l + l * a_dim1] / s;
u2 = a[l1 + l * a_dim1] / 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__1 = enorn;
for (j = l; j <= i__1; ++j) {
t = a[l + j * a_dim1] + u2 * a[l1 + j * a_dim1];
a[l + j * a_dim1] += t * v1;
a[l1 + j * a_dim1] += t * v2;
t = b[l + j * b_dim1] + u2 * b[l1 + j * b_dim1];
b[l + j * b_dim1] += t * v1;
b[l1 + j * b_dim1] += t * v2;
/* L110: */
}
if (l != 1) {
a[l + lm1 * a_dim1] = -a[l + lm1 * a_dim1];
}
lm1 = l;
l = l1;
goto L90;
L120:
a11 = a[l + l * a_dim1] / b11;
a21 = a[l1 + l * a_dim1] / b11;
if (ish == 1) {
goto L140;
}
/* .......... ITERATION STRATEGY .......... */
if (itn == 0) {
goto L1000;
}
if (its == 10) {
goto L155;
}
/* .......... DETERMINE TYPE OF SHIFT .......... */
b22 = b[l1 + l1 * b_dim1];
if (fabs(b22) < epsb) {
b22 = epsb;
}
b33 = b[na + na * b_dim1];
if (fabs(b33) < epsb) {
b33 = epsb;
}
b44 = b[en + en * b_dim1];
if (fabs(b44) < epsb) {
b44 = epsb;
}
a33 = a[na + na * a_dim1] / b33;
a34 = a[na + en * a_dim1] / b44;
a43 = a[en + na * a_dim1] / b33;
a44 = a[en + en * a_dim1] / b44;
b34 = b[na + en * b_dim1] / b44;
t = (a43 * b34 - a33 - a44) * .5;
r__ = t * t + a34 * a43 - a33 * a44;
if (r__ < 0.) {
goto L150;
}
/* .......... DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A .......... */
ish = 1;
r__ = sqrt(r__);
sh = -t + r__;
s = -t - r__;
if ((d__1 = s - a44, fabs(d__1)) < (d__2 = sh - a44, fabs(d__2))) {
sh = s;
}
/* .......... LOOK FOR TWO CONSECUTIVE SMALL */
/* SUB-DIAGONAL ELEMENTS OF A. */
/* FOR L=EN-2 STEP -1 UNTIL LD DO -- .......... */
i__1 = enm2;
for (ll = ld; ll <= i__1; ++ll) {
l = enm2 + ld - ll;
if (l == ld) {
goto L140;
}
lm1 = l - 1;
l1 = l + 1;
t = a[l + l * a_dim1];
if ((d__1 = b[l + l * b_dim1], fabs(d__1)) > epsb) {
t -= sh * b[l + l * b_dim1];
}
if ((d__1 = a[l + lm1 * a_dim1], fabs(d__1)) <= (d__2 = t / a[l1 + l *
a_dim1], fabs(d__2)) * epsa) {
goto L100;
}
/* L130: */
}
L140:
a1 = a11 - sh;
a2 = a21;
if (l != ld) {
a[l + lm1 * a_dim1] = -a[l + lm1 * a_dim1];
}
goto L160;
/* .......... DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A .......... */
L150:
a12 = a[l + l1 * a_dim1] / b22;
a22 = a[l1 + l1 * a_dim1] / b22;
b12 = b[l + l1 * b_dim1] / b22;
a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11) / a21 +
a12 - a11 * b12;
a2 = a22 - a11 - a21 * b12 - (a33 - a11) - (a44 - a11) + a43 * b34;
a3 = a[l1 + 1 + l1 * a_dim1] / b22;
goto L160;
/* .......... AD HOC SHIFT .......... */
L155:
a1 = 0.;
a2 = 1.;
a3 = 1.1605;
L160:
++its;
--itn;
if (! (matz)) {
lor1 = ld;
}
/* .......... MAIN LOOP .......... */
i__1 = na;
for (k = l; k <= i__1; ++k) {
notlas = k != na && ish == 2;
k1 = k + 1;
k2 = k + 2;
/* Computing MAX */
i__2 = k - 1;
km1 = max(i__2,l);
/* Computing MIN */
i__2 = en, i__3 = k1 + ish;
ll = min(i__2,i__3);
if (notlas) {
goto L190;
}
/* .......... ZERO A(K+1,K-1) .......... */
if (k == l) {
goto L170;
}
a1 = a[k + km1 * a_dim1];
a2 = a[k1 + km1 * a_dim1];
L170:
s = fabs(a1) + fabs(a2);
if (s == 0.) {
goto L70;
}
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 = enorn;
for (j = km1; j <= i__2; ++j) {
t = a[k + j * a_dim1] + u2 * a[k1 + j * a_dim1];
a[k + j * a_dim1] += t * v1;
a[k1 + j * a_dim1] += t * v2;
t = b[k + j * b_dim1] + u2 * b[k1 + j * b_dim1];
b[k + j * b_dim1] += t * v1;
b[k1 + j * b_dim1] += t * v2;
/* L180: */
}
if (k != l) {
a[k1 + km1 * a_dim1] = 0.;
}
goto L240;
/* .......... ZERO A(K+1,K-1) AND A(K+2,K-1) .......... */
L190:
if (k == l) {
goto L200;
}
a1 = a[k + km1 * a_dim1];
a2 = a[k1 + km1 * a_dim1];
a3 = a[k2 + km1 * a_dim1];
L200:
s = fabs(a1) + fabs(a2) + fabs(a3);
if (s == 0.) {
goto L260;
}
u1 = a1 / s;
u2 = a2 / s;
u3 = a3 / s;
d__1 = sqrt(u1 * u1 + u2 * u2 + u3 * u3);
r__ = d_sign(d__1, u1);
v1 = -(u1 + r__) / r__;
v2 = -u2 / r__;
v3 = -u3 / r__;
u2 = v2 / v1;
u3 = v3 / v1;
i__2 = enorn;
for (j = km1; j <= i__2; ++j) {
t = a[k + j * a_dim1] + u2 * a[k1 + j * a_dim1] + u3 * a[k2 + j *
a_dim1];
a[k + j * a_dim1] += t * v1;
a[k1 + j * a_dim1] += t * v2;
a[k2 + j * a_dim1] += t * v3;
t = b[k + j * b_dim1] + u2 * b[k1 + j * b_dim1] + u3 * b[k2 + j *
b_dim1];
b[k + j * b_dim1] += t * v1;
b[k1 + j * b_dim1] += t * v2;
b[k2 + j * b_dim1] += t * v3;
/* L210: */
}
if (k == l) {
goto L220;
}
a[k1 + km1 * a_dim1] = 0.;
a[k2 + km1 * a_dim1] = 0.;
/* .......... ZERO B(K+2,K+1) AND B(K+2,K) .......... */
L220:
s = (d__1 = b[k2 + k2 * b_dim1], fabs(d__1)) + (d__2 = b[k2 + k1 *
b_dim1], fabs(d__2)) + (d__3 = b[k2 + k * b_dim1], fabs(d__3));
if (s == 0.) {
goto L240;
}
u1 = b[k2 + k2 * b_dim1] / s;
u2 = b[k2 + k1 * b_dim1] / s;
u3 = b[k2 + k * b_dim1] / s;
d__1 = sqrt(u1 * u1 + u2 * u2 + u3 * u3);
r__ = d_sign(d__1, u1);
v1 = -(u1 + r__) / r__;
v2 = -u2 / r__;
v3 = -u3 / r__;
u2 = v2 / v1;
u3 = v3 / v1;
i__2 = ll;
for (i__ = lor1; i__ <= i__2; ++i__) {
t = a[i__ + k2 * a_dim1] + u2 * a[i__ + k1 * a_dim1] + u3 * a[i__
+ k * a_dim1];
a[i__ + k2 * a_dim1] += t * v1;
a[i__ + k1 * a_dim1] += t * v2;
a[i__ + k * a_dim1] += t * v3;
t = b[i__ + k2 * b_dim1] + u2 * b[i__ + k1 * b_dim1] + u3 * b[i__
+ k * b_dim1];
b[i__ + k2 * b_dim1] += t * v1;
b[i__ + k1 * b_dim1] += t * v2;
b[i__ + k * b_dim1] += t * v3;
/* L230: */
}
b[k2 + k * b_dim1] = 0.;
b[k2 + k1 * b_dim1] = 0.;
if (! (matz)) {
goto L240;
}
i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
t = z__[i__ + k2 * z_dim1] + u2 * z__[i__ + k1 * z_dim1] + u3 *
z__[i__ + k * z_dim1];
z__[i__ + k2 * z_dim1] += t * v1;
z__[i__ + k1 * z_dim1] += t * v2;
z__[i__ + k * z_dim1] += t * v3;
/* L235: */
}
/* .......... ZERO B(K+1,K) .......... */
L240:
s = (d__1 = b[k1 + k1 * b_dim1], fabs(d__1)) + (d__2 = b[k1 + k *
b_dim1], fabs(d__2));
if (s == 0.) {
goto L260;
}
u1 = b[k1 + k1 * b_dim1] / s;
u2 = b[k1 + k * b_dim1] / 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 = ll;
for (i__ = lor1; i__ <= i__2; ++i__) {
t = a[i__ + k1 * a_dim1] + u2 * a[i__ + k * a_dim1];
a[i__ + k1 * a_dim1] += t * v1;
a[i__ + k * a_dim1] += t * v2;
t = b[i__ + k1 * b_dim1] + u2 * b[i__ + k * b_dim1];
b[i__ + k1 * b_dim1] += t * v1;
b[i__ + k * b_dim1] += t * v2;
/* L250: */
}
b[k1 + k * b_dim1] = 0.;
if (! (matz)) {
goto L260;
}
i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
t = z__[i__ + k1 * z_dim1] + u2 * z__[i__ + k * z_dim1];
z__[i__ + k1 * z_dim1] += t * v1;
z__[i__ + k * z_dim1] += t * v2;
/* L255: */
}
L260:
;
}
/* .......... END QZ STEP .......... */
goto L70;
/* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
/* CONVERGED AFTER 30*N ITERATIONS .......... */
L1000:
*ierr = en;
/* .......... SAVE EPSB FOR USE BY QZVAL AND QZVEC .......... */
L1001:
if (n > 1) {
b[n + b_dim1] = epsb;
}
return 0;
} /* qzit_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* 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;
/* Local variables */
doublereal epsb, c__, d__, e;
integer i__, j;
doublereal r__, s, t, a1, a2, u1, u2, v1, v2, a11, a12, a21, a22,
b11, b12, b22, di, ei;
integer na;
doublereal an, bn;
integer en;
doublereal cq, dr;
integer nn;
doublereal cz, ti, tr, a1i, a2i, a11i, a12i, a22i, a11r, a12r,
a22r, sqi, ssi;
integer isw;
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 */
/* 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], fabs(d__1));
alfi[en] = 0.;
goto L510;
/* .......... 2-BY-2 BLOCK .......... */
L420:
if ((d__1 = b[na + na * b_dim1], fabs(d__1)) <= epsb) {
goto L455;
}
if ((d__1 = b[en + en * b_dim1], fabs(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], fabs(d__1)) + (d__2 = a[na + en *
a_dim1], fabs(d__2)) + (d__3 = a[en + na * a_dim1], fabs(d__3))
+ (d__4 = a[en + en * a_dim1], fabs(d__4));
bn = (d__1 = b[na + na * b_dim1], fabs(d__1)) + (d__2 = b[na + en *
b_dim1], fabs(d__2)) + (d__3 = b[en + en * b_dim1], fabs(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 (fabs(e) <= fabs(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 (fabs(a11) + fabs(a12) < fabs(a21) + fabs(a22)) {
goto L432;
}
a1 = a12;
a2 = a11;
goto L435;
L432:
a1 = a22;
a2 = a21;
/* .......... CHOOSE AND APPLY REAL Z .......... */
L435:
s = fabs(a1) + fabs(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 < fabs(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 = fabs(a1) + fabs(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], fabs(d__1));
beta[en] = (d__1 = b[en + en * b_dim1], fabs(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 (fabs(a11r) + fabs(a11i) + fabs(a12r) + fabs(a12i) < fabs(a21) + fabs(
a22r) + fabs(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 < (fabs(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_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
doublereal
epslon(doublereal x)
{
/* System generated locals */
doublereal ret_val, d__1;
/* Local variables */
doublereal a, b, c__, eps;
/* ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. */
/* THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS */
/* SATISFYING THE FOLLOWING TWO ASSUMPTIONS, */
/* 1. THE BASE USED IN REPRESENTING FLOATING POINT */
/* NUMBERS IS NOT A POWER OF THREE. */
/* 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO */
/* THE ACCURACY USED IN FLOATING POINT VARIABLES */
/* THAT ARE STORED IN MEMORY. */
/* THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO */
/* FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING */
/* ASSUMPTION 2. */
/* UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, */
/* A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, */
/* B HAS A ZERO FOR ITS LAST BIT OR DIGIT, */
/* C IS NOT EXACTLY EQUAL TO ONE, */
/* EPS MEASURES THE SEPARATION OF 1.0 FROM */
/* THE NEXT LARGER FLOATING POINT NUMBER. */
/* THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED */
/* ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD. */
/* THIS VERSION DATED 4/6/83. */
a = 1.3333333333333333;
L10:
b = a - 1.;
c__ = b + b + b;
eps = (d__1 = c__ - 1., fabs(d__1));
if (eps == 0.) {
goto L10;
}
ret_val = eps * fabs(x);
return ret_val;
} /* epslon_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* BLAS-1 routines needed in the computation of Floquet multipliers */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
doublereal
dnrm2(integer *n, doublereal *dx, integer *incx)
{
/* Initialized data */
#define zero 0.
#define one 1.
#define cutlo 8.232e-11
#define cuthi 1.304e19
/* Format strings */
/* System generated locals */
integer i__1, i__2;
doublereal ret_val, d__1;
/* Local variables */
doublereal xmax;
integer next, i__, j, nn;
doublereal hitest, sum;
/* Parameter adjustments */
--dx;
/* Function Body */
/* EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE */
/* INCREMENT INCX . */
/* IF N .LE. 0 RETURN WITH RESULT = 0. */
/* IF N .GE. 1 THEN INCX MUST BE .GE. 1 */
/* C.L.LAWSON, 1978 JAN 08 */
/* FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE */
/* HOPEFULLY APPLICABLE TO ALL MACHINES. */
/* CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. */
/* CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES. */
/* WHERE */
/* EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. */
/* U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) */
/* V = LARGEST NO. (OVERFLOW LIMIT) */
/* BRIEF OUTLINE OF ALGORITHM.. */
/* PHASE 1 SCANS ZERO COMPONENTS. */
/* MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO */
/* MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO */
/* MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M */
/* WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. */
/* VALUES FOR CUTLO AND CUTHI.. */
/* DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. */
/* CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE
*/
/* UNIVAC AND DEC AT 2**(-103) */
/* THUS CUTLO = 2**(-51) = 4.44089E-16 */
/* CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. */
/* THUS CUTHI = 2**(63.5) = 1.30438E19 */
/* CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. */
/* THUS CUTLO = 2**(-33.5) = 8.23181D-11 */
/* CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 */
/* DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / */
/* DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / */
if (*n > 0) {
goto L10;
}
ret_val = zero;
goto L300;
L10:
next = 0;
sum = zero;
nn = *n * *incx;
/* BEGIN MAIN LOOP */
i__ = 1;
L20:
switch ((int)next) {
case 0: goto L30;
case 1: goto L50;
case 2: goto L70;
case 3: goto L110;
}
L30:
if ((d__1 = dx[i__], fabs(d__1)) > cutlo) {
goto L85;
}
next = 1;
xmax = zero;
/* PHASE 1. SUM IS ZERO */
L50:
if (dx[i__] == zero) {
goto L200;
}
if ((d__1 = dx[i__], fabs(d__1)) > cutlo) {
goto L85;
}
/* PREPARE FOR PHASE 2. */
next = 2;
goto L105;
/* PREPARE FOR PHASE 4. */
L100:
i__ = j;
next = 3;
sum = sum / dx[i__] / dx[i__];
L105:
xmax = (d__1 = dx[i__], fabs(d__1));
goto L115;
/* PHASE 2. SUM IS SMALL. */
/* SCALE TO AVOID DESTRUCTIVE UNDERFLOW. */
L70:
if ((d__1 = dx[i__], fabs(d__1)) > cutlo) {
goto L75;
}
/* COMMON CODE FOR PHASES 2 AND 4. */
/* IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
*/
L110:
if ((d__1 = dx[i__], fabs(d__1)) <= xmax) {
goto L115;
}
/* Computing 2nd power */
d__1 = xmax / dx[i__];
sum = one + sum * (d__1 * d__1);
xmax = (d__1 = dx[i__], fabs(d__1));
goto L200;
L115:
/* Computing 2nd power */
d__1 = dx[i__] / xmax;
sum += d__1 * d__1;
goto L200;
/* PREPARE FOR PHASE 3. */
L75:
sum = sum * xmax * xmax;
/* FOR REAL OR D.P. SET HITEST = CUTHI/N */
/* FOR COMPLEX SET HITEST = CUTHI/(2*N) */
L85:
hitest = cuthi / (real) (*n);
/* PHASE 3. SUM IS MID-RANGE. NO SCALING. */
i__1 = nn;
i__2 = *incx;
for (j = i__; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
if ((d__1 = dx[j], fabs(d__1)) >= hitest) {
goto L100;
}
/* L95: */
/* Computing 2nd power */
d__1 = dx[j];
sum += d__1 * d__1;
}
ret_val = sqrt(sum);
goto L300;
L200:
i__ += *incx;
if (i__ <= nn) {
goto L20;
}
/* END OF MAIN LOOP. */
/* COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. */
ret_val = xmax * sqrt(sum);
L300:
return ret_val;
} /* dnrm2_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
doublereal
ddot(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
{
/* System generated locals */
integer i__1;
doublereal ret_val;
/* Local variables */
integer i__, m;
doublereal dtemp;
integer ix, iy, mp1;
/* FORMS THE DOT PRODUCT OF TWO VECTORS. */
/* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
/* JACK DONGARRA, LINPACK, 3/11/78. */
/* Parameter adjustments */
--dy;
--dx;
/* Function Body */
ret_val = 0.;
dtemp = 0.;
if (*n <= 0) {
return ret_val;
}
if (*incx == 1 && *incy == 1) {
goto L20;
}
/* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
/* NOT EQUAL TO 1 */
ix = 1;
iy = 1;
if (*incx < 0) {
ix = (-(*n) + 1) * *incx + 1;
}
if (*incy < 0) {
iy = (-(*n) + 1) * *incy + 1;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
dtemp += dx[ix] * dy[iy];
ix += *incx;
iy += *incy;
/* L10: */
}
ret_val = dtemp;
return ret_val;
/* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
/* CLEAN-UP LOOP */
L20:
m = *n % 5;
if (m == 0) {
goto L40;
}
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__) {
dtemp += dx[i__] * dy[i__];
/* L30: */
}
if (*n < 5) {
goto L60;
}
L40:
mp1 = m + 1;
i__1 = *n;
for (i__ = mp1; i__ <= i__1; i__ += 5) {
dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ +
4] * dy[i__ + 4];
/* L50: */
}
L60:
ret_val = dtemp;
return ret_val;
} /* ddot_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
dscal(integer *n, doublereal *da, doublereal *dx, integer *incx)
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
integer i__, m, nincx, mp1;
/* SCALES A VECTOR BY A CONSTANT. */
/* USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. */
/* JACK DONGARRA, LINPACK, 3/11/78. */
/* Parameter adjustments */
--dx;
/* Function Body */
if (*n <= 0) {
return 0;
}
if (*incx == 1) {
goto L20;
}
/* CODE FOR INCREMENT NOT EQUAL TO 1 */
nincx = *n * *incx;
i__1 = nincx;
i__2 = *incx;
for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
dx[i__] = *da * dx[i__];
/* L10: */
}
return 0;
/* CODE FOR INCREMENT EQUAL TO 1 */
/* CLEAN-UP LOOP */
L20:
m = *n % 5;
if (m == 0) {
goto L40;
}
i__2 = m;
for (i__ = 1; i__ <= i__2; ++i__) {
dx[i__] = *da * dx[i__];
/* L30: */
}
if (*n < 5) {
return 0;
}
L40:
mp1 = m + 1;
i__2 = *n;
for (i__ = mp1; i__ <= i__2; i__ += 5) {
dx[i__] = *da * dx[i__];
dx[i__ + 1] = *da * dx[i__ + 1];
dx[i__ + 2] = *da * dx[i__ + 2];
dx[i__ + 3] = *da * dx[i__ + 3];
dx[i__ + 4] = *da * dx[i__ + 4];
/* L50: */
}
return 0;
} /* dscal_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
integer
idamax(integer *n, doublereal *dx, integer *incx)
{
/* System generated locals */
integer ret_val, i__1;
doublereal d__1;
/* Local variables */
doublereal dmax__;
integer i__, ix;
/* FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. */
/* JACK DONGARRA, LINPACK, 3/11/78. */
/* Parameter adjustments */
--dx;
/* Function Body */
ret_val = 0;
if (*n < 1) {
return ret_val;
}
ret_val = 1;
if (*n == 1) {
return ret_val;
}
if (*incx == 1) {
goto L20;
}
/* CODE FOR INCREMENT NOT EQUAL TO 1 */
ix = 1;
dmax__ = fabs(dx[1]);
ix += *incx;
i__1 = *n;
for (i__ = 2; i__ <= i__1; ++i__) {
if ((d__1 = dx[ix], fabs(d__1)) <= dmax__) {
goto L5;
}
ret_val = i__;
dmax__ = (d__1 = dx[ix], fabs(d__1));
L5:
ix += *incx;
/* L10: */
}
return ret_val;
/* CODE FOR INCREMENT EQUAL TO 1 */
L20:
dmax__ = fabs(dx[1]);
i__1 = *n;
for (i__ = 2; i__ <= i__1; ++i__) {
if ((d__1 = dx[i__], fabs(d__1)) <= dmax__) {
goto L30;
}
ret_val = i__;
dmax__ = (d__1 = dx[i__], fabs(d__1));
L30:
;
}
return ret_val;
} /* idamax_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
daxpy(integer *n, doublereal *da, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
{
/* System generated locals */
integer i__1;
/* Local variables */
integer i__, m, ix, iy, mp1;
/* CONSTANT TIMES A VECTOR PLUS A VECTOR. */
/* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
/* JACK DONGARRA, LINPACK, 3/11/78. */
/* Parameter adjustments */
--dy;
--dx;
/* Function Body */
if (*n <= 0) {
return 0;
}
if (*da == 0.) {
return 0;
}
if (*incx == 1 && *incy == 1) {
goto L20;
}
/* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
/* NOT EQUAL TO 1 */
ix = 1;
iy = 1;
if (*incx < 0) {
ix = (-(*n) + 1) * *incx + 1;
}
if (*incy < 0) {
iy = (-(*n) + 1) * *incy + 1;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
dy[iy] += *da * dx[ix];
ix += *incx;
iy += *incy;
/* L10: */
}
return 0;
/* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
/* CLEAN-UP LOOP */
L20:
m = *n % 4;
if (m == 0) {
goto L40;
}
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__) {
dy[i__] += *da * dx[i__];
/* L30: */
}
if (*n < 4) {
return 0;
}
L40:
mp1 = m + 1;
i__1 = *n;
for (i__ = mp1; i__ <= i__1; i__ += 4) {
dy[i__] += *da * dx[i__];
dy[i__ + 1] += *da * dx[i__ + 1];
dy[i__ + 2] += *da * dx[i__ + 2];
dy[i__ + 3] += *da * dx[i__ + 3];
/* L50: */
}
return 0;
} /* daxpy_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
drot(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy, doublereal *c__, doublereal *s)
{
/* System generated locals */
integer i__1;
/* Local variables */
integer i__;
doublereal dtemp;
integer ix, iy;
/* APPLIES A PLANE ROTATION. */
/* JACK DONGARRA, LINPACK, 3/11/78. */
/* Parameter adjustments */
--dy;
--dx;
/* Function Body */
if (*n <= 0) {
return 0;
}
if (*incx == 1 && *incy == 1) {
goto L20;
}
/* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL */
/* TO 1 */
ix = 1;
iy = 1;
if (*incx < 0) {
ix = (-(*n) + 1) * *incx + 1;
}
if (*incy < 0) {
iy = (-(*n) + 1) * *incy + 1;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
dtemp = *c__ * dx[ix] + *s * dy[iy];
dy[iy] = *c__ * dy[iy] - *s * dx[ix];
dx[ix] = dtemp;
ix += *incx;
iy += *incy;
/* L10: */
}
return 0;
/* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
L20:
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
dtemp = *c__ * dx[i__] + *s * dy[i__];
dy[i__] = *c__ * dy[i__] - *s * dx[i__];
dx[i__] = dtemp;
/* L30: */
}
return 0;
} /* drot_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
dswap(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
{
/* System generated locals */
integer i__1;
/* Local variables */
integer i__, m;
doublereal dtemp;
integer ix, iy, mp1;
/* INTERCHANGES TWO VECTORS. */
/* USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE. */
/* JACK DONGARRA, LINPACK, 3/11/78. */
/* Parameter adjustments */
--dy;
--dx;
/* Function Body */
if (*n <= 0) {
return 0;
}
if (*incx == 1 && *incy == 1) {
goto L20;
}
/* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL */
/* TO 1 */
ix = 1;
iy = 1;
if (*incx < 0) {
ix = (-(*n) + 1) * *incx + 1;
}
if (*incy < 0) {
iy = (-(*n) + 1) * *incy + 1;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
dtemp = dx[ix];
dx[ix] = dy[iy];
dy[iy] = dtemp;
ix += *incx;
iy += *incy;
/* L10: */
}
return 0;
/* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
/* CLEAN-UP LOOP */
L20:
m = *n % 3;
if (m == 0) {
goto L40;
}
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__) {
dtemp = dx[i__];
dx[i__] = dy[i__];
dy[i__] = dtemp;
/* L30: */
}
if (*n < 3) {
return 0;
}
L40:
mp1 = m + 1;
i__1 = *n;
for (i__ = mp1; i__ <= i__1; i__ += 3) {
dtemp = dx[i__];
dx[i__] = dy[i__];
dy[i__] = dtemp;
dtemp = dx[i__ + 1];
dx[i__ + 1] = dy[i__ + 1];
dy[i__ + 1] = dtemp;
dtemp = dx[i__ + 2];
dx[i__ + 2] = dy[i__ + 2];
dy[i__ + 2] = dtemp;
/* L50: */
}
return 0;
} /* dswap_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
dgemc(integer *m, integer *n, doublereal *a, integer *lda, doublereal *b, integer *ldb, logical *trans)
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, i1, i2;
/* Local variables */
integer i, j, mm, mmp1;
/* This subroutine copies a double precision real */
/* M by N matrix stored in A to double precision real B. */
/* If TRANS is true, B is assigned A transpose. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = a_dim1 + 1;
a -= a_offset;
b_dim1 = *ldb;
/* Function Body */
if (*trans) {
i1 = *n;
for (j = 1; j <= i1; ++j) {
/* USES UNROLLED LOOPS */
/* from JACK DONGARRA, LINPACK, 3/11/78. */
mm = *m % 7;
if (mm == 0) {
goto L80;
}
i2 = mm;
for (i = 1; i <= i2; ++i) {
b[j-1 + (i-1) * b_dim1] = a[i + j * a_dim1];
/* L70: */
}
if (*m < 7) {
goto L99;
}
L80:
mmp1 = mm + 1;
i2 = *m;
for (i = mmp1; i <= i2; i += 7) {
b[j-1 + (i-1) * b_dim1] = a[i + j * a_dim1];
b[j-1 + (i-1 + 1) * b_dim1] = a[i + 1 + j * a_dim1];
b[j-1 + (i-1 + 2) * b_dim1] = a[i + 2 + j * a_dim1];
b[j-1 + (i-1 + 3) * b_dim1] = a[i + 3 + j * a_dim1];
b[j-1 + (i-1 + 4) * b_dim1] = a[i + 4 + j * a_dim1];
b[j-1 + (i-1 + 5) * b_dim1] = a[i + 5 + j * a_dim1];
b[j-1 + (i-1 + 6) * b_dim1] = a[i + 6 + j * a_dim1];
/* L90: */
}
L99:
/* L100: */
;
}
} else {
i1 = *n;
for (j = 1; j <= i1; ++j) {
/* USES UNROLLED LOOPS */
/* from JACK DONGARRA, LINPACK, 3/11/78. */
mm = *m % 7;
if (mm == 0) {
goto L180;
}
i2 = mm;
for (i = 1; i <= i2; ++i) {
b[(j-1) * b_dim1 + i-1] = a[i + j * a_dim1];
/* L170: */
}
if (*m < 7) {
goto L199;
}
L180:
mmp1 = mm + 1;
i2 = *m;
for (i = mmp1; i <= i2; i += 7) {
b[i-1 + (j-1) * b_dim1] = a[i + j * a_dim1];
b[i-1 + 1 + (j-1) * b_dim1] = a[i + 1 + j * a_dim1];
b[i-1 + 2 + (j-1) * b_dim1] = a[i + 2 + j * a_dim1];
b[i-1 + 3 + (j-1) * b_dim1] = a[i + 3 + j * a_dim1];
b[i-1 + 4 + (j-1) * b_dim1] = a[i + 4 + j * a_dim1];
b[i-1 + 5 + (j-1) * b_dim1] = a[i + 5 + j * a_dim1];
b[i-1 + 6 + (j-1) * b_dim1] = a[i + 6 + j * a_dim1];
/* L190: */
}
L199:
/* L200: */
;
}
}
return 0;
} /* dgemc_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* BLAS-2 routines needed in the computation of Floquet multipliers */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
xerbla(char *srname, integer *info, integer srname_len)
{
/* Format strings */
/* Builtin functions */
/* Fortran I/O blocks */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* XERBLA is an error handler for the Level 2 BLAS routines. */
/* It is called by the Level 2 BLAS routines if an input parameter is */
/* invalid. */
/* Installers should consider modifying the STOP statement in order to */
/* call system-specific exception-handling facilities. */
/* Parameters */
/* ========== */
/* SRNAME - CHARACTER*6. */
/* On entry, SRNAME specifies the name of the routine which */
/* called XERBLA. */
/* INFO - INTEGER. */
/* On entry, INFO specifies the position of the invalid */
/* parameter in the parameter-list of the calling routine. */
/* Auxiliary routine for Level 2 Blas. */
/* Written on 20-July-1986. */
/* .. Executable Statements .. */
fprintf(fp6, "On entry to %c%c%c%c%c%c parameter number %ld had an illegal value\n",
srname[0],srname[1],srname[2],srname[3],srname[4],srname[5],(*info));
EXTERNAL_LIBRARY_FATAL_ERROR("O8AUTO", 0);
return 0;
/* End of XERBLA. */
} /* xerbla_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
logical
lsame(char *ca, char *cb, integer ca_len, integer cb_len)
{
/* System generated locals */
logical ret_val;
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* LSAME tests if CA is the same letter as CB regardless of case. */
/* CB is assumed to be an upper case letter. LSAME returns .TRUE. if */
/* CA is either the same as CB or the equivalent lower case letter. */
/* N.B. This version of the routine is only correct for ASCII code. */
/* Installers must modify the routine for other character-codes. */
/* For EBCDIC systems the constant IOFF must be changed to -64. */
/* For CDC systems using 6-12 bit representations, the system- */
/* specific code in comments must be activated. */
/* Parameters */
/* ========== */
/* CA - CHARACTER*1 */
/* CB - CHARACTER*1 */
/* On entry, CA and CB specify characters to be compared. */
/* Unchanged on exit. */
/* Auxiliary routine for Level 2 Blas. */
/* -- Written on 20-July-1986 */
/* Richard Hanson, Sandia National Labs. */
/* Jeremy Du Croz, Nag Central Office. */
/* .. Parameters .. */
/* .. Intrinsic Functions .. */
/* .. Executable Statements .. */
/* Test if the characters are equal */
ret_val = *(unsigned char *)ca == *(unsigned char *)cb;
/* Now test for equivalence */
if (! ret_val) {
ret_val = *(unsigned char *)ca - 32 == *(unsigned char *)cb;
}
return ret_val;
/* The following comments contain code for CDC systems using 6-12 bit */
/* representations. */
/* .. Parameters .. */
/* INTEGER ICIRFX */
/* .. Scalar Arguments .. */
/* CHARACTER*1 CB */
/* .. Array Arguments .. */
/* CHARACTER*1 CA(*) */
/* .. Local Scalars .. */
/* INTEGER IVAL */
/* .. Intrinsic Functions .. */
/* INTRINSIC ICHAR, CHAR */
/* .. Executable Statements .. */
/* See if the first character in string CA equals string CB. */
/* LSAME = CA(1) .EQ. CB .AND. CA(1) .NE. CHAR(ICIRFX) */
/* IF (LSAME) RETURN */
/* The characters are not identical. Now check them for equivalence.
*/
/* Look for the 'escape' character, circumflex, followed by the */
/* letter. */
/* IVAL = ICHAR(CA(2)) */
/* IF (IVAL.GE.ICHAR('A') .AND. IVAL.LE.ICHAR('Z')) THEN */
/* LSAME = CA(1) .EQ. CHAR(ICIRFX) .AND. CA(2) .EQ. CB */
/* END IF */
/* RETURN */
/* End of LSAME. */
} /* lsame_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* BLAS-3 routines needed in the computation of Floquet multipliers */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
dgemm(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c__, integer *ldc, integer transa_len, integer transb_len)
{
/* System generated locals */
integer a_dim1, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2, i__3;
/* Local variables */
integer info;
logical nota, notb;
doublereal temp;
integer i__, j, l, ncola;
integer nrowa, nrowb;
/* .. Scalar Arguments .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGEMM performs one of the matrix-matrix operations */
/* C := alpha*op( A )*op( B ) + beta*C, */
/* where op( X ) is one of */
/* op( X ) = X or op( X ) = X', */
/* alpha and beta are scalars, and A, B and C are matrices, with op( A )
*/
/* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
*/
/* Parameters */
/* ========== */
/* TRANSA - CHARACTER*1. */
/* On entry, TRANSA specifies the form of op( A ) to be used in
*/
/* the matrix multiplication as follows: */
/* TRANSA = 'N' or 'n', op( A ) = A. */
/* TRANSA = 'T' or 't', op( A ) = A'. */
/* TRANSA = 'C' or 'c', op( A ) = A'. */
/* Unchanged on exit. */
/* TRANSB - CHARACTER*1. */
/* On entry, TRANSB specifies the form of op( B ) to be used in
*/
/* the matrix multiplication as follows: */
/* TRANSB = 'N' or 'n', op( B ) = B. */
/* TRANSB = 'T' or 't', op( B ) = B'. */
/* TRANSB = 'C' or 'c', op( B ) = B'. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix
*/
/* op( A ) and of the matrix C. M must be at least zero.
*/
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix
*/
/* op( B ) and the number of columns of the matrix C. N must be
*/
/* at least zero. */
/* Unchanged on exit. */
/* K - INTEGER. */
/* On entry, K specifies the number of columns of the matrix
*/
/* op( A ) and the number of rows of the matrix op( B ). K must
*/
/* be at least zero. */
/* Unchanged on exit. */
/* ALPHA - DOUBLE PRECISION. */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
*/
/* k when TRANSA = 'N' or 'n', and is m otherwise. */
/* Before entry with TRANSA = 'N' or 'n', the leading m by k
*/
/* part of the array A must contain the matrix A, otherwise
*/
/* the leading k by m part of the array A must contain the
*/
/* matrix A. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared
*/
/* in the calling (sub) program. When TRANSA = 'N' or 'n' then
*/
/* LDA must be at least max( 1, m ), otherwise LDA must be at
*/
/* least max( 1, k ). */
/* Unchanged on exit. */
/* B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
*/
/* n when TRANSB = 'N' or 'n', and is k otherwise. */
/* Before entry with TRANSB = 'N' or 'n', the leading k by n
*/
/* part of the array B must contain the matrix B, otherwise
*/
/* the leading n by k part of the array B must contain the
*/
/* matrix B. */
/* Unchanged on exit. */
/* LDB - INTEGER. */
/* On entry, LDB specifies the first dimension of B as declared
*/
/* in the calling (sub) program. When TRANSB = 'N' or 'n' then
*/
/* LDB must be at least max( 1, k ), otherwise LDB must be at
*/
/* least max( 1, n ). */
/* Unchanged on exit. */
/* BETA - DOUBLE PRECISION. */
/* On entry, BETA specifies the scalar beta. When BETA is
*/
/* supplied as zero then C need not be set on input. */
/* Unchanged on exit. */
/* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). */
/* Before entry, the leading m by n part of the array C must
*/
/* contain the matrix C, except when beta is zero, in which
*/
/* case C need not be set on entry. */
/* On exit, the array C is overwritten by the m by n matrix
*/
/* ( alpha*op( A )*op( B ) + beta*C ). */
/* LDC - INTEGER. */
/* On entry, LDC specifies the first dimension of C as declared
*/
/* in the calling (sub) program. LDC must be at least
*/
/* max( 1, m ). */
/* Unchanged on exit. */
/* Level 3 Blas routine. */
/* -- Written on 8-February-1989. */
/* Jack Dongarra, Argonne National Laboratory. */
/* Iain Duff, AERE Harwell. */
/* Jeremy Du Croz, Numerical Algorithms Group Ltd. */
/* Sven Hammarling, Numerical Algorithms Group Ltd. */
/* .. External Functions .. */
/* .. External Subroutines .. */
/* .. Intrinsic Functions .. */
/* .. Local Scalars .. */
/* .. Parameters .. */
/* .. */
/* .. Executable Statements .. */
/* Set NOTA and NOTB as true if A and B respectively are not
*/
/* transposed and set NROWA, NCOLA and NROWB as the number of rows
*/
/* and columns of A and the number of rows of B respectively.
*/
/* Parameter adjustments */
a_dim1 = *lda;
b_dim1 = *ldb;
b_offset = b_dim1 + 1;
b -= b_offset;
c_dim1 = *ldc;
c_offset = c_dim1 + 1;
c__ -= c_offset;
/* Function Body */
nota = lsame(transa, "N", 1L, 1L);
notb = lsame(transb, "N", 1L, 1L);
if (nota) {
nrowa = *m;
ncola = *k;
} else {
nrowa = *k;
ncola = *m;
}
if (notb) {
nrowb = *k;
} else {
nrowb = *n;
}
/* Test the input parameters. */
info = 0;
if (! nota && ! lsame(transa, "C", 1L, 1L) && ! lsame(transa, "T", 1L,
1L)) {
info = 1;
} else if (! notb && ! lsame(transb, "C", 1L, 1L) && ! lsame(transb,
"T", 1L, 1L)) {
info = 2;
} else if (*m < 0) {
info = 3;
} else if (*n < 0) {
info = 4;
} else if (*k < 0) {
info = 5;
} else if (*lda < max(1,nrowa)) {
info = 8;
} else if (*ldb < max(1,nrowb)) {
info = 10;
} else if (*ldc < max(1,*m)) {
info = 13;
}
if (info != 0) {
xerbla("DGEMM ", &info, 6L);
return 0;
}
/* Quick return if possible. */
if (*m == 0 || *n == 0 || ((*alpha == 0. || *k == 0) && *beta == 1.)) {
return 0;
}
/* And if alpha.eq.zero. */
if (*alpha == 0.) {
if (*beta == 0.) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
c__[i__ + j * c_dim1] = 0.;
/* L10: */
}
/* L20: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
/* L30: */
}
/* L40: */
}
}
return 0;
}
/* Start the operations. */
if (notb) {
if (nota) {
/* Form C := alpha*A*B + beta*C. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (*beta == 0.) {
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
c__[i__ + j * c_dim1] = 0.;
/* L50: */
}
} else if (*beta != 1.) {
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
/* L60: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l) {
if (b[l + j * b_dim1] != 0.) {
temp = *alpha * b[l + j * b_dim1];
i__3 = *m;
for (i__ = 1; i__ <= i__3; ++i__) {
c__[i__ + j * c_dim1] += temp * a[i__-1 + (l-1) * a_dim1];
/* L70: */
}
}
/* L80: */
}
/* L90: */
}
} else {
/* Form C := alpha*A'*B + beta*C */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
temp = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l) {
temp += a[l-1 + (i__-1) * a_dim1] * b[l + j * b_dim1];
/* L100: */
}
if (*beta == 0.) {
c__[i__ + j * c_dim1] = *alpha * temp;
} else {
c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
i__ + j * c_dim1];
}
/* L110: */
}
/* L120: */
}
}
} else {
if (nota) {
/* Form C := alpha*A*B' + beta*C */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (*beta == 0.) {
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
c__[i__ + j * c_dim1] = 0.;
/* L130: */
}
} else if (*beta != 1.) {
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
/* L140: */
}
}
i__2 = *k;
for (l = 1; l <= i__2; ++l) {
if (b[j + l * b_dim1] != 0.) {
temp = *alpha * b[j + l * b_dim1];
i__3 = *m;
for (i__ = 1; i__ <= i__3; ++i__) {
c__[i__ + j * c_dim1] += temp * a[i__-1 + (l-1) * a_dim1];
/* L150: */
}
}
/* L160: */
}
/* L170: */
}
} else {
/* Form C := alpha*A'*B' + beta*C */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
temp = 0.;
i__3 = *k;
for (l = 1; l <= i__3; ++l) {
temp += a[l-1 + (i__-1) * a_dim1] * b[j + l * b_dim1];
/* L180: */
}
if (*beta == 0.) {
c__[i__ + j * c_dim1] = *alpha * temp;
} else {
c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
i__ + j * c_dim1];
}
/* L190: */
}
/* L200: */
}
}
}
return 0;
/* End of DGEMM . */
} /* dgemm_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Demmel-Kahan SVD routines needed for computing the Floquet multipliers */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
ezsvd(doublereal *x, integer *ldx, integer *n, integer *p, doublereal *s, doublereal *e, doublereal *u, integer *ldu, doublereal *v, integer *ldv, doublereal *work, integer *job, integer *info, doublereal *tol)
{
/* System generated locals */
integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset;
/* Local variables */
integer idbg, skip, iidir, ifull;
integer kount, kount1, kount2, limshf;
doublereal maxsin;
integer maxitr;
/* new svd by J. Demmel, W. Kahan */
/* finds singular values of bidiagonal matrices with guaranteed high
*/
/* relative precision */
/* easy to use version of ndsvd ("hard to use" version, below) */
/* with defaults for some ndsvd parameters */
/* all parameters same as linpack dsvdc except for tol: */
/* tol = if positive, desired relative precision in singular values
*/
/* if negative, desired absolute precision in singular values
*/
/* (expressed as abs(tol) * sigma-max) */
/* (in both cases, abs(tol) should be less than 1 and */
/* greater than macheps) */
/* I have tested this software on a SUN 3 in double precision */
/* IEEE arithmetic with macheps about 2.2e-16 and tol=1e-14; */
/* In general I recommend tol 10-100 times larger than macheps. */
/* On the average it appears to be as fast or faster than dsvdc. */
/* I have seen it go 3.5 times faster and 2 times slower at the */
/* extremes. */
/* defaults for ndsvd parameters (see ndsvd for more description of */
/* these parameters) are: */
/* set to no debug output */
/* Parameter adjustments */
x_dim1 = *ldx;
x_offset = x_dim1 + 1;
x -= x_offset;
--s;
--e;
u_dim1 = *ldu;
u_offset = u_dim1 + 1;
u -= u_offset;
v_dim1 = *ldv;
v_offset = v_dim1 + 1;
v -= v_offset;
--work;
/* Function Body */
idbg = 0;
/* use zero-shift normally */
ifull = 0;
/* use normal bidiagonalization code */
skip = 0;
/* choose chase direction normally */
iidir = 0;
/* maximum 30 QR sweeps per singular value */
maxitr = 30;
ndsvd(&x[x_offset], ldx, n, p, &s[1], &e[1], &u[u_offset], ldu, &v[
v_offset], ldv, &work[1], job, info, &maxitr, tol, &idbg, &ifull,
&kount, &kount1, &kount2, &skip, &limshf, &maxsin, &iidir);
return 0;
} /* ezsvd_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
ndrotg(doublereal *f, doublereal *g, doublereal *cs, doublereal *sn)
{
/* Local variables */
doublereal t, tt;
/* faster version of drotg, except g unchanged on return */
/* cs, sn returned so that -sn*f+cs*g = 0 */
/* and returned f = cs*f + sn*g */
/* if g=0, then cs=1 and sn=0 (in case svd adds extra zero row */
/* to bidiagonal, this makes sure last row rotation is trivial) */
/* if f=0 and g.ne.0, then cs=0 and sn=1 without floating point work
*/
/* (in case s(i)=0 in svd so that bidiagonal deflates, this */
/* computes rotation without any floating point operations) */
if (*f == 0.) {
if (*g == 0.) {
/* this case needed in case extra zero row added in svd, s
o */
/* bottom rotation always trivial */
*cs = 1.;
*sn = 0.;
} else {
/* this case needed for s(i)=0 in svd to compute rotation
*/
/* cheaply */
*cs = 0.;
*sn = 1.;
*f = *g;
}
} else {
if (fabs(*f) > fabs(*g)) {
t = *g / *f;
tt = sqrt(t * t + 1.);
*cs = 1. / tt;
*sn = t * *cs;
*f *= tt;
} else {
t = *f / *g;
tt = sqrt(t * t + 1.);
*sn = 1. / tt;
*cs = t * *sn;
*f = *g * tt;
}
}
return 0;
} /* ndrotg_ */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int
ndsvd(doublereal *x, integer *ldx, integer *n, integer *p, doublereal *s, doublereal *e, doublereal *u, integer *ldu, doublereal *v, integer *ldv, doublereal *work, integer *job, integer *info, integer *maxitr, doublereal *tol, integer *idbg, integer *ifull, integer *kount, integer *kount1, integer *kount2, integer *skip, integer *limshf, doublereal *maxsin, integer *iidir)
{
/* System generated locals */
integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2,
i__3;
doublereal d__1, d__2, d__3, d__4;
/* Local variables */
doublereal abse;
integer idir;
doublereal abss;
integer oldm, jobu;
doublereal cosl;
integer iter;
doublereal temp, smin, smax, cosr, sinl, sinr;
doublereal test;
integer nctp1, nrtp1;
doublereal f, g;
integer i__, j, k, l, m;
doublereal t;
doublereal oldcs;
integer oldll, iisub;
doublereal shift, oldsn, sigmn;
integer minnp, maxit;
doublereal sminl;
doublereal sigmx;
logical wantu, wantv;
doublereal gg, lambda;
integer oldacc;
doublereal cs;
integer ll, mm;
doublereal sm;
integer lu;
doublereal sn, mu;
doublereal thresh;
integer lm1, lp1, lll, nct, ncu;
doublereal sll;
integer nrt;
doublereal emm1, smm1;
/* Fortran I/O blocks */
/* LINPACK SVD modified by: */
/* James Demmel W. Kahan */
/* Courant Institute Computer Science Dept. */
/* demmel@acf8.nyu.edu U.C. Berkeley */
/* modified version designed to guarantee relative accuracy of */
/* all singular values of intermediate bidiagonal form */
/* extra input/output parameters in addition to those from LINPACK SVD:
*/
/* extra input paramters: */
/* tol = if positive, desired relative precision in singular values
*/
/* if negative, desired absolute precision in singular values
*/
/* (expressed as abs(tol) * sigma-max) */
/* (abs(tol) should be less than 1 and greater than macheps) */
/* idbg = 0 for no debug output (normal setting) */
/* = 1 convergence, shift decisions (written to standard output)
*/
/* = 2 for above plus before, after qr */
/* ifull= 0 if decision to use zero-shift set normally (normal setting)
*/
/* = 1 if always set to nonzero-shift */
/* = 2 if always set to zero-shift */
/* skip =-1 means standard code but do all work of bidiagonalization
*/
/* (even if input bidiagonal) */
/* 0 means standard code (normal setting) */
/* 1 means assume x is bidiagonal, and skip bidiagonalization
*/
/* entirely */
/* (skip used for timing tests) */
/* iidir = 0 if idir (chase direction) chosen normally */
/* 1 if idir=1 (chase top to bottom) always */
/* 2 if idir=2 (chase bottom to top) always */
/* extra output parameters: */
/* kount =number of qr sweeps taken */
/* kount1=number of passes through inner loop of full qr */
/* kount2=number of passes through inner loop of zero-shift qr */
/* limshf = number of times the shift was greater than its threshold
*/
/* (nct*smin) and had to be decreased */
/* maxsin = maximum sin in inner loop of zero-shift */
/* new version designed to be robust with respect to over/underflow */
/* have fast inner loop when shift is zero, */
/* guarantee relative accuracy of all singular values */
/* dsvdc is a subroutine to reduce a double precision nxp matrix x */
/* by orthogonal transformations u and v to diagonal form. the */
/* diagonal elements s(i) are the singular values of x. the */
/* columns of u are the corresponding left singular vectors, */
/* and the columns of v the right singular vectors. */
/* on entry */
/* x double precision(ldx,p), where ldx.ge.n. */
/* x contains the matrix whose singular value */
/* decomposition is to be computed. x is */
/* destroyed by dsvdc. */
/* ldx integer. */
/* ldx is the leading dimension of the array x. */
/* n integer. */
/* n is the number of rows of the matrix x. */
/* p integer. */
/* p is the number of columns of the matrix x. */
/* ldu integer. */
/* ldu is the leading dimension of the array u. */
/* (see below). */
/* ldv integer. */
/* ldv is the leading dimension of the array v. */
/* (see below). */
/* work double precision(n). */
/* work is a scratch array. */
/* job integer. */
/* job controls the computation of the singular */
/* vectors. it has the decimal expansion ab */
/* with the following meaning */
/* a.eq.0 do not compute the left singular */
/* vectors. */
/* a.eq.1 return the n left singular vectors */
/* in u. */
/* a.ge.2 return the first min(n,p) singular */
/* vectors in u. */
/* blas daxpy,ddot,dscal,dswap,dnrm2,drotg */
/* fortran dabs,dmax1,max0,min0,mod,dsqrt,dsign */
/* prse,ndrotg,sig22,sndrtg,sigmin */
/* internal variables */
/* new variables */
/* double precision sg1,sg2 */
/* set the maximum number of iterations. */
/* maxit = 30 */
/* Parameter adjustments */
x_dim1 = *ldx;
x_offset = x_dim1 + 1;
x -= x_offset;
--s;
--e;
u_dim1 = *ldu;
u_offset = u_dim1 + 1;
u -= u_offset;
v_dim1 = *ldv;
v_offset = v_dim1 + 1;
v -= v_offset;
--work;
/* Function Body */
*kount = 0;
*kount1 = 0;
*kount2 = 0;
*limshf = 0;
*maxsin = (double)0.;
/* determine what is to be computed. */
wantu = FALSE_;
wantv = FALSE_;
jobu = *job % 100 / 10;
ncu = *n;
if (jobu > 1) {
ncu = min(*n,*p);
}
if (jobu != 0) {
wantu = TRUE_;
}
if (*job % 10 != 0) {
wantv = TRUE_;
}
/* reduce x to bidiagonal form, storing the diagonal elements */
/* in s and the super-diagonal elements in e. */
*info = 0;
/* Computing MIN */
i__1 = *n - 1;
nct = min(i__1,*p);
/* Computing MAX */
/* Computing MIN */
i__3 = *p - 2;
i__1 = 0, i__2 = min(i__3,*n);
nrt = max(i__1,i__2);
lu = max(nct,nrt);
if (*skip <= 0) {
if (lu < 1) {
goto L170;
}
i__1 = lu;
for (l = 1; l <= i__1; ++l) {
integer c__1 = 1;
lp1 = l + 1;
if (l > nct) {
goto L20;
}
/* compute the transformation for the l-th column and */
/* place the l-th diagonal in s(l). */
i__2 = *n - l + 1;
s[l] = dnrm2(&i__2, &x[l + l * x_dim1], &c__1);
if (s[l] == 0. && *skip == 0) {
goto L10;
}
if (x[l + l * x_dim1] != 0.) {
s[l] = d_sign(s[l], x[l + l * x_dim1]);
}
i__2 = *n - l + 1;
d__1 = 1. / s[l];
dscal(&i__2, &d__1, &x[l + l * x_dim1], &c__1);
x[l + l * x_dim1] += 1.;
L10:
s[l] = -s[l];
L20:
if (*p < lp1) {
goto L50;
}
i__2 = *p;
for (j = lp1; j <= i__2; ++j) {
integer c__1 = 1;
if (l > nct) {
goto L30;
}
if (s[l] == 0. && *skip == 0) {
goto L30;
}
/* apply the transformation. */
i__3 = *n - l + 1;
t = -ddot(&i__3, &x[l + l * x_dim1], &c__1, &x[l + j *
x_dim1], &c__1) / x[l + l * x_dim1];
i__3 = *n - l + 1;
daxpy(&i__3, &t, &x[l + l * x_dim1], &c__1, &x[l + j *
x_dim1], &c__1);
L30:
/* place the l-th row of x into e for the */
/* subsequent calculation of the row transformation.
*/
e[j] = x[l + j * x_dim1];
/* L40: */
}
L50:
if (! wantu || l > nct) {
goto L70;
}
/* place the transformation in u for subsequent back */
/* multiplication. */
i__2 = *n;
for (i__ = l; i__ <= i__2; ++i__) {
u[i__ + l * u_dim1] = x[i__ + l * x_dim1];
/* L60: */
}
L70:
if (l > nrt) {
goto L150;
}
/* compute the l-th row transformation and place the */
/* l-th super-diagonal in e(l). */
i__2 = *p - l;
e[l] = dnrm2(&i__2, &e[lp1], &c__1);
if (e[l] == 0. && *skip == 0) {
goto L80;
}
if (e[lp1] != 0.) {
e[l] = d_sign(e[l], e[lp1]);
}
i__2 = *p - l;
d__1 = 1. / e[l];
dscal(&i__2, &d__1, &e[lp1], &c__1);
e[lp1] += 1.;
L80:
e[l] = -e[l];
if (lp1 > *n || (e[l] == 0. && *skip == 0)) {
goto L120;
}
/* apply the transformation. */
i__2 = *n;
for (i__ = lp1; i__ <= i__2; ++i__) {
work[i__] = 0.;
/* L90: */
}
i__2 = *p;
for (j = lp1; j <= i__2; ++j) {
i__3 = *n - l;
daxpy(&i__3, &e[j], &x[lp1 + j * x_dim1], &c__1, &work[lp1],
&c__1);
/* L100: */
}
i__2 = *p;
for (j = lp1; j <= i__2; ++j) {
i__3 = *n - l;
d__1 = -e[j] / e[lp1];
daxpy(&i__3, &d__1, &work[lp1], &c__1, &x[lp1 + j * x_dim1],
&c__1);
/* L110: */
}
L120:
if (! wantv) {
goto L140;
}
/* place the transformation in v for subsequent */
/* back multiplication. */
i__2 = *p;
for