/* cbal.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 cbal_(integer *nm, integer *n, doublereal *ar,
doublereal *ai, integer *low, integer *igh, doublereal *scale)
{
/* System generated locals */
integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2;
doublereal d__1, d__2;
/* Local variables */
static integer iexc;
static doublereal c__, f, g;
static integer i__, j, k, l, m;
static doublereal r__, s, radix, b2;
static integer jj;
static logical noconv;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE */
/* CBALANCE, WHICH IS A COMPLEX VERSION OF 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 COMPLEX MATRIX AND ISOLATES */
/* EIGENVALUES WHENEVER POSSIBLE. */
/* 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. */
/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
/* RESPECTIVELY, OF THE COMPLEX MATRIX TO BE BALANCED. */
/* ON OUTPUT */
/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
/* RESPECTIVELY, OF THE BALANCED MATRIX. */
/* LOW AND IGH ARE TWO INTEGERS SUCH THAT AR(I,J) AND AI(I,J) */
/* ARE 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 CBALANCE APPEARS IN */
/* CBAL IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS */
/* K,L HAVE BEEN REVERSED.) */
/* ARITHMETIC IS REAL THROUGHOUT. */
/* 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;
ai_dim1 = *nm;
ai_offset = ai_dim1 + 1;
ai -= ai_offset;
ar_dim1 = *nm;
ar_offset = ar_dim1 + 1;
ar -= ar_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 = ar[i__ + j * ar_dim1];
ar[i__ + j * ar_dim1] = ar[i__ + m * ar_dim1];
ar[i__ + m * ar_dim1] = f;
f = ai[i__ + j * ai_dim1];
ai[i__ + j * ai_dim1] = ai[i__ + m * ai_dim1];
ai[i__ + m * ai_dim1] = f;
/* L30: */
}
i__1 = *n;
for (i__ = k; i__ <= i__1; ++i__) {
f = ar[j + i__ * ar_dim1];
ar[j + i__ * ar_dim1] = ar[m + i__ * ar_dim1];
ar[m + i__ * ar_dim1] = f;
f = ai[j + i__ * ai_dim1];
ai[j + i__ * ai_dim1] = ai[m + i__ * ai_dim1];
ai[m + i__ * ai_dim1] = f;
/* L40: */
}
L50:
switch (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 (ar[j + i__ * ar_dim1] != 0. || ai[j + i__ * ai_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 (ar[i__ + j * ar_dim1] != 0. || ai[i__ + j * ai_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__ = c__ + (d__1 = ar[j + i__ * ar_dim1], abs(d__1)) + (d__2 =
ai[j + i__ * ai_dim1], abs(d__2));
r__ = r__ + (d__1 = ar[i__ + j * ar_dim1], abs(d__1)) + (d__2 =
ai[i__ + j * ai_dim1], abs(d__2));
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) {
ar[i__ + j * ar_dim1] *= g;
ai[i__ + j * ai_dim1] *= g;
/* L250: */
}
i__2 = l;
for (j = 1; j <= i__2; ++j) {
ar[j + i__ * ar_dim1] *= f;
ai[j + i__ * ai_dim1] *= f;
/* L260: */
}
L270:
;
}
if (noconv) {
goto L190;
}
L280:
*low = k;
*igh = l;
return 0;
} /* cbal_ */
syntax highlighted by Code2HTML, v. 0.9.1