/* comhes.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 comhes_(integer *nm, integer *n, integer *low, integer *
igh, doublereal *ar, doublereal *ai, integer *int__)
{
/* System generated locals */
integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3;
doublereal d__1, d__2;
/* Local variables */
extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
, doublereal *, doublereal *, doublereal *);
static integer i__, j, m, la;
static doublereal xi, yi, xr, yr;
static integer mm1, kp1, mp1;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMHES, */
/* NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
/* GIVEN A COMPLEX 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 */
/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRIX. */
/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
/* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */
/* SET LOW=1, IGH=N. */
/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
/* RESPECTIVELY, OF THE COMPLEX INPUT MATRIX. */
/* ON OUTPUT */
/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
/* RESPECTIVELY, OF THE HESSENBERG MATRIX. THE */
/* MULTIPLIERS WHICH WERE USED IN THE REDUCTION */
/* ARE STORED IN THE REMAINING TRIANGLES UNDER THE */
/* HESSENBERG MATRIX. */
/* INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
/* INTERCHANGED IN THE REDUCTION. */
/* ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
/* 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 */
ai_dim1 = *nm;
ai_offset = ai_dim1 + 1;
ai -= ai_offset;
ar_dim1 = *nm;
ar_offset = ar_dim1 + 1;
ar -= ar_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;
xr = 0.;
xi = 0.;
i__ = m;
i__2 = *igh;
for (j = m; j <= i__2; ++j) {
if ((d__1 = ar[j + mm1 * ar_dim1], abs(d__1)) + (d__2 = ai[j +
mm1 * ai_dim1], abs(d__2)) <= abs(xr) + abs(xi)) {
goto L100;
}
xr = ar[j + mm1 * ar_dim1];
xi = ai[j + mm1 * ai_dim1];
i__ = j;
L100:
;
}
int__[m] = i__;
if (i__ == m) {
goto L130;
}
/* .......... INTERCHANGE ROWS AND COLUMNS OF AR AND AI ..........
*/
i__2 = *n;
for (j = mm1; j <= i__2; ++j) {
yr = ar[i__ + j * ar_dim1];
ar[i__ + j * ar_dim1] = ar[m + j * ar_dim1];
ar[m + j * ar_dim1] = yr;
yi = ai[i__ + j * ai_dim1];
ai[i__ + j * ai_dim1] = ai[m + j * ai_dim1];
ai[m + j * ai_dim1] = yi;
/* L110: */
}
i__2 = *igh;
for (j = 1; j <= i__2; ++j) {
yr = ar[j + i__ * ar_dim1];
ar[j + i__ * ar_dim1] = ar[j + m * ar_dim1];
ar[j + m * ar_dim1] = yr;
yi = ai[j + i__ * ai_dim1];
ai[j + i__ * ai_dim1] = ai[j + m * ai_dim1];
ai[j + m * ai_dim1] = yi;
/* L120: */
}
/* .......... END INTERCHANGE .......... */
L130:
if (xr == 0. && xi == 0.) {
goto L180;
}
mp1 = m + 1;
i__2 = *igh;
for (i__ = mp1; i__ <= i__2; ++i__) {
yr = ar[i__ + mm1 * ar_dim1];
yi = ai[i__ + mm1 * ai_dim1];
if (yr == 0. && yi == 0.) {
goto L160;
}
cdiv_(&yr, &yi, &xr, &xi, &yr, &yi);
ar[i__ + mm1 * ar_dim1] = yr;
ai[i__ + mm1 * ai_dim1] = yi;
i__3 = *n;
for (j = m; j <= i__3; ++j) {
ar[i__ + j * ar_dim1] = ar[i__ + j * ar_dim1] - yr * ar[m + j
* ar_dim1] + yi * ai[m + j * ai_dim1];
ai[i__ + j * ai_dim1] = ai[i__ + j * ai_dim1] - yr * ai[m + j
* ai_dim1] - yi * ar[m + j * ar_dim1];
/* L140: */
}
i__3 = *igh;
for (j = 1; j <= i__3; ++j) {
ar[j + m * ar_dim1] = ar[j + m * ar_dim1] + yr * ar[j + i__ *
ar_dim1] - yi * ai[j + i__ * ai_dim1];
ai[j + m * ai_dim1] = ai[j + m * ai_dim1] + yr * ai[j + i__ *
ai_dim1] + yi * ar[j + i__ * ar_dim1];
/* L150: */
}
L160:
;
}
L180:
;
}
L200:
return 0;
} /* comhes_ */
syntax highlighted by Code2HTML, v. 0.9.1