/* qzhes.f -- translated by f2c (version 19961017).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Subroutine */ int 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;
/* Builtin functions */
double sqrt(doublereal), d_sign(doublereal *, doublereal *);
/* Local variables */
static integer i__, j, k, l;
static doublereal r__, s, t;
static integer l1;
static doublereal u1, u2, v1, v2;
static integer lb, nk1, nm1, nm2;
static 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 */
/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
/* 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], abs(d__1));
/* L20: */
}
if (s == 0.) {
goto L100;
}
s += (d__1 = b[l + l * b_dim1], abs(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], abs(d__1)) + (d__2 = a[l1 + k *
a_dim1], abs(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], abs(d__1)) + (d__2 = b[l1 + l *
b_dim1], abs(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_ */
syntax highlighted by Code2HTML, v. 0.9.1