rg.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,457 行 · 第 1/3 页
C
1,457 行
#include "f2c.h"
#include "netlib.h"
extern double sqrt(double); /* #include <math.h> */
static void balanc_(integer *nm, integer *n, doublereal *a, integer *low, integer *igh, doublereal *scale);
static void balbak_(integer *nm, integer *n, integer *low, integer *igh, doublereal *scale, integer *m, doublereal *z);
static void cdiv_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci);
static void elmhes_(integer *nm, integer *n, integer *low, integer *igh, doublereal *a, integer *int_);
static void eltran_(integer *nm, integer *n, integer *low, integer *igh, doublereal *a, integer *int_, doublereal *z);
static void hqr_(integer *nm, integer *n, integer *low, integer *igh, doublereal *h,
doublereal *wr, doublereal *wi, integer *ierr);
static void hqr2_(integer *nm, integer *n, integer *low, integer *igh, doublereal *h,
doublereal *wr, doublereal *wi, doublereal *z, integer *ierr);
/* Modified by Peter Vanroose, Sept 2001: manual optimisation and clean-up */
/* Table of constant values */
static doublereal c_b130 = 0.;
/* ====================================================================== */
/* NIST Guide to Available Math Software. */
/* Fullsource for module RG from package EISPACK. */
/* Retrieved from NETLIB on Thu Jan 23 06:12:53 1997. */
/* ====================================================================== */
/* Subroutine */ void rg_(nm, n, a, wr, wi, matz, z, iv1, fv1, ierr)
integer *nm, *n;
doublereal *a, *wr, *wi;
integer *matz;
doublereal *z;
integer *iv1;
doublereal *fv1;
integer *ierr;
{
/* Local variables */
static 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 */
/* array parameters as declared in the calling program */
/* 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. */
/* ------------------------------------------------------------------ */
if (*n <= *nm) {
goto L10;
}
*ierr = *n * 10;
goto L50;
L10:
balanc_(nm, n, a, &is1, &is2, fv1);
elmhes_(nm, n, &is1, &is2, a, iv1);
if (*matz != 0) {
goto L20;
}
/* .......... find eigenvalues only .......... */
hqr_(nm, n, &is1, &is2, a, wr, wi, ierr);
goto L50;
/* .......... find both eigenvalues and eigenvectors .......... */
L20:
eltran_(nm, n, &is1, &is2, a, iv1, z);
hqr2_(nm, n, &is1, &is2, a, wr, wi, z, ierr);
if (*ierr != 0) {
goto L50;
}
balbak_(nm, n, &is1, &is2, fv1, n, z);
L50:
return;
} /* rg_ */
/* Subroutine */ void balanc_(nm, n, a, low, igh, scale)
integer *nm, *n;
doublereal *a;
integer *low, *igh;
doublereal *scale;
{
/* Local variables */
static integer iexc;
static doublereal c, f, g;
static integer i, j, k, l, m;
static doublereal r, s, radix, b2;
static 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 */
/* array parameters as declared in the calling program */
/* 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. */
/* ------------------------------------------------------------------ */
radix = 16.;
b2 = radix * radix;
k = 0;
l = *n;
goto L100;
/* .......... in-line procedure for row and column exchange .......... */
L20:
scale[m] = (doublereal) j+1;
if (j == m) {
goto L50;
}
for (i = 0; i < l; ++i) {
f = a[i + j * *nm];
a[i + j * *nm] = a[i + m * *nm];
a[i + m * *nm] = f;
}
for (i = k; i < *n; ++i) {
f = a[j + i * *nm];
a[j + i * *nm] = a[m + i * *nm];
a[m + i * *nm] = f;
}
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:
for (j = l-1; j >= 0; --j) {
for (i = 0; i < l; ++i) {
if (i != j && a[j + i * *nm] != 0.) {
goto L120; /* continue outer loop */
}
}
m = l-1;
iexc = 1;
goto L20;
L120:
;
}
goto L140;
/* .......... search for columns isolating an eigenvalue and push them left .......... */
L130:
++k;
L140:
for (j = k; j < l; ++j) {
for (i = k; i < l; ++i) {
if (i != j && a[i + j * *nm] != 0.) {
goto L170; /* continue outer loop */
}
}
m = k;
iexc = 2;
goto L20;
L170:
;
}
/* .......... now balance the submatrix in rows k to l .......... */
for (i = k; i < l; ++i) {
scale[i] = 1.;
}
/* .......... iterative loop for norm reduction .......... */
L190:
noconv = FALSE_;
for (i = k; i < l; ++i) {
c = 0.;
r = 0.;
for (j = k; j < l; ++j) {
if (j != i) {
c += abs(a[j + i * *nm]);
r += abs(a[i + j * *nm]);
}
}
/* .......... guard against zero c or r due to underflow .......... */
if (c == 0. || r == 0.) {
continue;
}
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) {
continue;
}
g = 1. / f;
scale[i] *= f;
noconv = TRUE_;
for (j = k; j < *n; ++j) {
a[i + j * *nm] *= g;
}
for (j = 0; j < l; ++j) {
a[j + i * *nm] *= f;
}
}
if (noconv) {
goto L190;
}
L280:
*low = k+1;
*igh = l;
return;
} /* balanc_ */
/* Subroutine */ void balbak_(nm, n, low, igh, scale, m, z)
integer *nm, *n, *low, *igh;
doublereal *scale;
integer *m;
doublereal *z;
{
/* Local variables */
static integer i, j, k;
static doublereal s;
static 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 */
/* array parameters as declared in the calling program */
/* 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. */
/* ------------------------------------------------------------------ */
if (*m == 0) {
return; /* exit from balbak_ */
}
if (*igh == *low) {
goto L120;
}
for (i = *low-1; i < *igh; ++i) {
s = scale[i];
/* .......... left hand eigenvectors are back transformed */
/* if the foregoing statement is replaced by s=1.0d0/scale(i). .......... */
for (j = 0; j < *m; ++j) {
z[i + j * *nm] *= s;
}
}
/* ......... for i=low-1 step -1 until 1, */
/* igh+1 step 1 until n do -- .......... */
L120:
for (ii = 0; ii < *n; ++ii) {
i = ii;
if (i+1 >= *low && i < *igh) {
continue;
}
if (i+1 < *low) {
i = *low - ii - 2;
}
k = (integer) scale[i] - 1;
if (k != i)
for (j = 0; j < *m; ++j) {
s = z[i + j * *nm];
z[i + j * *nm] = z[k + j * *nm];
z[k + j * *nm] = s;
}
}
} /* balbak_ */
/* Subroutine */ void cdiv_(ar, ai, br, bi, cr, ci)
doublereal *ar, *ai, *br, *bi, *cr, *ci;
{
/* Local variables */
static doublereal s, ais, bis, ars, brs;
/* complex division, (cr,ci) = (ar,ai)/(br,bi) */
s = abs(*br) + abs(*bi);
ars = *ar / s;
ais = *ai / s;
brs = *br / s;
bis = *bi / s;
s = brs * brs + bis * bis;
*cr = (ars * brs + ais * bis) / s;
*ci = (ais * brs - ars * bis) / s;
} /* cdiv_ */
/* Subroutine */ void elmhes_(nm, n, low, igh, a, int_)
integer *nm, *n, *low, *igh;
doublereal *a;
integer *int_;
{
/* Local variables */
static integer i, j, m;
static doublereal x, y;
static integer la, mm1;
/* 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 */
/* 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 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. */
/* ------------------------------------------------------------------ */
la = *igh - 1;
if (la < *low + 1) {
return; /* exit from elmhes_ */
}
for (m = *low; m < la; ++m) {
mm1 = m-1;
x = 0.;
i = m;
for (j = m; j < *igh; ++j) {
if (abs(a[j + mm1 * *nm]) > abs(x)) {
x = a[j + mm1 * *nm];
i = j;
}
}
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?