#include <stdlib.h>
#include "f2c.h"
#include <math.h>
#define MAX(a,b) ((a)>(b)?(a):(b))
#define COMPZERO 1e-13
extern int NODE;
extern double NEWT_ERR;
int (*rhs)();
typedef struct {
int nunstab,nstab,n; /* dimension of unstable, stable, phasespace */
int eleft,eright; /* addresses of variables holding left and right
equilibria */
int u0; /* address of first phase variable */
double cprev[1600],a[400];
int iflag[4];
double fb[400]; /* values of the boundary projections
first nstab are proj to unstable mfld
at left ane then the proj to the stabl
mfld on the right */
} HOMOCLIN;
HOMOCLIN my_hom;
static integer c__1 = 1;
static integer c__2 = 2;
static integer c__5 = 5;
static integer c_n1 = -1;
static integer c__9 = 9;
static integer c__0 = 0;
static integer c__20 = 20;
double hom_bcs(double x)
{
int i=(int)x;
if(i>=my_hom.n)return 0.0;
return my_hom.fb[i];
}
do_projection(double *x0,double t0,double *x1,double t1)
{
double y[400],f[400],fnew[400];
double bound[400];
double eps=NEWT_ERR,del,dsy,yold;
double *fb;
int i,j,n=my_hom.n,nunstab=my_hom.nunstab,nstab=my_hom.nstab;
int imfld=-1,itrans=2,is=1; /* stored in transpose form
so dont transpose */
/* first we take care of the left */
/* printf("n=%d nu=%d ns=%d no=%d \n",
n,nunstab,nstab,NODE); */
for(i=0;i<n;i++){
y[i]=x0[i+my_hom.eleft];
}
for(i=n;i<NODE;i++)
y[i]=x0[i];
/* Jacobian around the left equilibrium */
rhs(t0,y,f,NODE);
for(i=0;i<n;i++){
del=eps*MAX(eps,fabs(y[i]));
dsy=1/del;
yold=y[i];
y[i]=y[i]+del;
rhs(t0,y,fnew,NODE);
for(j=0;j<n;j++){
my_hom.a[j*n+i]=dsy*(fnew[j]-f[j]);
/* printf("%d %d a=%g \n",i,j,my_hom.a[j*n+i]); */
}
y[i]=yold;
}
projection_(bound, &imfld, &is, &itrans);
for(i=0;i<nstab;i++){
my_hom.fb[i]=0.0;
for(j=0;j<n;j++){
my_hom.fb[i]+=((x0[j]-y[j])*bound[i+j*20]);
}
}
/* okay - thats all the BCs associated with the left
end. I.e. The projection onto the UNSTABLE MFLD
*/
/* Now the right side !! */
for(i=0;i<n;i++)
y[i]=x0[i+my_hom.eright];
for(i=n;i<NODE;i++)
y[i]=x0[i];
/* Jacobian around the right equilibrium */
rhs(t0,y,f,NODE);
for(i=0;i<n;i++){
del=eps*MAX(eps,fabs(y[i]));
dsy=1/del;
yold=y[i];
y[i]=y[i]+del;
rhs(t0,y,fnew,NODE);
for(j=0;j<n;j++)
my_hom.a[j*n+i]=dsy*(fnew[j]-f[j]);
y[i]=yold;
}
imfld=1;
is=2;
projection_(bound, &imfld, &is, &itrans);
for(i=0;i<nunstab;i++){
my_hom.fb[i+nstab]=0.0;
for(j=0;j<n;j++){
my_hom.fb[i+nstab]+=((x1[j]-y[j])*bound[i+nstab+j*20]);
}
}
}
int pdfdu_(double *a,int n)
{
int i,j;
for(j=0;j<n;j++)
for(i=0;i<n;i++){
a[i+j*20]=my_hom.a[i+j*n];
}
}
/* ----------------------------------------------------- */
/* Subroutine */ int projection_(bound, imfd, is, itrans)
doublereal *bound;
integer *imfd, *is, *itrans;
{
/* System generated locals */
integer i__1, i__2, i__3;
/* Local variables */
static doublereal cnow[400] /* was [20][20] */;
static integer type__[20];
static doublereal a[400] /* was [20][20] */, d__[400] /* was [20][
20] */;
extern /* Subroutine */ int f04aef_();
static integer i__, j, k, n, ifail;
static doublereal v[400] /* was [20][20] */;
static integer mcond;
extern /* Subroutine */ int pdfdu_();
static integer k1, k2, m0;
static doublereal aa[400] /* was [20][20] */, bb[400] /* was [20][
20] */, ei[20], er[20];
extern /* Subroutine */ int orthes_(), ortran_(), hqr3loc_();
static doublereal eps, ort[20], wkspace[20], dum1[400] /* was [20][
20] */, dum2[400] /* was [20][20] */;
/* ----------------------------------------------------- */
/* Compute NUNSTAB (or NSTAB) projection boundary condition functions */
/*onto to the UNSTABLE (or STABLE) manifold of the appropriate equilibrium
*/
/* IMFD = -1 stable eigenspace */
/* = 1 unstable eigenspace */
/* ITRANS = 1 use transpose of A */
/* = 2 otherwise */
/* IS = I (1 or 2) implies use the ith equilibrium in XEQUIB */
/* using the normalization in Beyn 1990 (4.4) to ensure */
/* continuity w.r.t parameters (using NAG routine F04AEF). */
/* For the purposes of this routine the "previous point on the */
/* branch" is at the values of par at which the routine was last */
/* called with the same values of IS and ITRANS. */
/* Parameter adjustments */
bound -= 21;
/* Function Body */
eps = COMPZERO; /* Machine epsilon */
n = my_hom.n;
pdfdu_(a,n);
/* compute transpose of A if ITRANS=1 */
if (*itrans == 1) {
i__1 = my_hom.n;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = my_hom.n;
for (j = 1; j <= i__2; ++j) {
bb[i__ + j * 20 - 21] = a[j + i__ * 20 - 21];
}
}
i__1 = my_hom.n;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = my_hom.n;
for (j = 1; j <= i__2; ++j) {
a[i__ + j * 20 - 21] = bb[i__ + j * 20 - 21];
}
}
}
/* Compute basis V to put A in upper Hessenberg form */
orthes_(&c__20, &n, &c__1, &n, a, ort);
ortran_(&c__20, &n, &c__1, &n, a, ort, v);
/* force A to be upper Hessenberg */
if (n > 2) {
i__1 = n;
for (i__ = 3; i__ <= i__1; ++i__) {
i__2 = i__ - 2;
for (j = 1; j <= i__2; ++j) {
a[i__ + j * 20 - 21] = 0.;
}
}
}
/* Computes basis to put A in "Quasi Upper-Triangular form" */
/* with the positive (negative) eigenvalues first if IMFD =-1 (=1) */
hqr3loc_(a, v, &c__20, &c__1, &n, &eps, er, ei, type__, &c__20, &c__20,
imfd);
/* put the basis in the appropriate part of the matrix CNOW */
if (*imfd == 1) {
k1 = n - my_hom.nunstab + 1;
k2 = n;
} else {
k1 = 1;
k2 = my_hom.nstab;
}
mcond = k2 - k1 + 1;
m0 = k1 - 1;
i__1 = k2;
for (i__ = k1; i__ <= i__1; ++i__) {
i__2 = n;
for (j = 1; j <= i__2; ++j) {
cnow[i__ + j * 20 - 21] = v[j + (i__ - k1 + 1) * 20 - 21];
}
}
/* set previous matrix to be the present one if this is the first call */
if (my_hom.iflag[*is + (*itrans << 1) - 3] == 0) {
i__1 = k2;
for (i__ = k1; i__ <= i__1; ++i__) {
i__2 = my_hom.n;
for (j = 1; j <= i__2; ++j) {
my_hom.cprev[i__ + (j + (*is + (*itrans << 1)) * 20) * 20 -
1221] = cnow[i__ + j * 20 - 21];
bound[i__ + j * 20] = cnow[i__ + j * 20 - 21];
/* printf(" i=%d j=%d b=%g \n",i__,j,bound[i__ + j * 20]); */
}
}
my_hom.iflag[*is + (*itrans << 1) - 3] = 1;
goto L400;
}
/* calculate the (transpose of the) BEYN matrix D and hence BOUND */
i__1 = mcond;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = mcond;
for (j = 1; j <= i__2; ++j) {
dum1[i__ + j * 20 - 21] = 0.;
dum2[i__ + j * 20 - 21] = 0.;
i__3 = my_hom.n;
for (k = 1; k <= i__3; ++k) {
dum1[i__ + j * 20 - 21] += my_hom.cprev[i__ + m0 + (k + (*is
+ (*itrans << 1)) * 20) * 20 - 1221] * cnow[j + m0 +
k * 20 - 21];
dum2[i__ + j * 20 - 21] += my_hom.cprev[i__ + m0 + (k + (*is
+ (*itrans << 1)) * 20) * 20 - 1221] * my_hom.cprev[j
+ m0 + (k + (*is + (*itrans << 1)) * 20) * 20 - 1221];
}
}
}
if (mcond != 1) {
ifail = 0;
f04aef_(dum1, &c__20, dum2, &c__20, &mcond, &mcond, d__, &c__20,
wkspace, aa, &c__20, bb, &c__20, &ifail);
} else {
d__[0] = dum2[0] / dum1[0];
}
i__1 = mcond;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = my_hom.n;
for (j = 1; j <= i__2; ++j) {
bound[i__ + m0 + j * 20] = (float)0.;
i__3 = mcond;
for (k = 1; k <= i__3; ++k) {
bound[i__ + m0 + j * 20] += d__[k + i__ * 20 - 21] * cnow[k +
m0 + j * 20 - 21];
}
}
}
i__1 = k2;
for (i__ = k1; i__ <= i__1; ++i__) {
i__2 = my_hom.n;
for (j = 1; j <= i__2; ++j) {
my_hom.cprev[i__ + (j + (*is + (*itrans << 1)) * 20) * 20 - 1221]
= bound[i__ + j * 20];
}
}
L400:
return 0;
} /* projection_ */
int hqr3loc_(a, v, n, nlow, nup, eps, er, ei, type__, na, nv,
imfd)
doublereal *a, *v;
integer *n, *nlow, *nup;
doublereal *eps, *er, *ei;
integer *type__, *na, *nv, *imfd;
{
/* System generated locals */
integer a_dim1, a_offset, v_dim1, v_offset, i__1;
doublereal d__1, d__2, d__3;
/* Local variables */
static logical fail;
static integer i__, l;
static doublereal p, q, r__, s, t, w, x, y, z__, e1, e2;
extern /* Subroutine */ int split_();
static integer nl, it, mu, nu;
extern /* Subroutine */ int exchng_(), qrstep_();
/* Parameter adjustments */
--type__;
--ei;
--er;
a_dim1 = *na;
a_offset = a_dim1 + 1;
a -= a_offset;
v_dim1 = *nv;
v_offset = v_dim1 + 1;
v -= v_offset;
/* Function Body */
i__1 = *nup;
for (i__ = *nlow; i__ <= i__1; ++i__) {
type__[i__] = -1;
/* L10: */
}
t = 0.;
/* MAIN LOOP. FIND AND ORDER EIGENVALUES. */
nu = *nup;
L20:
if (nu < *nlow) {
goto L240;
}
it = 0;
/* QR LOOP. FIND NEGLIGIBLE ELEMENTS AND PERFORM */
/* QR STEPS. */
L30:
/* SEARCH BACK FOR NEGLIGIBLE ELEMENTS. */
l = nu;
L40:
if (l == *nlow) {
goto L50;
}
if ((d__1 = a[l + (l - 1) * a_dim1], abs(d__1)) <= *eps * ((d__2 = a[l -
1 + (l - 1) * a_dim1], abs(d__2)) + (d__3 = a[l + l * a_dim1],
abs(d__3)))) {
goto L50;
}
--l;
goto L40;
L50:
/* TEST TO SEE IF AN EIGENVALUE OR A 2X2 BLOCK */
/* HAS BEEN FOUND. */
x = a[nu + nu * a_dim1];
if (l == nu) {
goto L160;
}
y = a[nu - 1 + (nu - 1) * a_dim1];
w = a[nu + (nu - 1) * a_dim1] * a[nu - 1 + nu * a_dim1];
if (l == nu - 1) {
goto L100;
}
/* TEST ITERATION COUNT. IF IT IS 30 QUIT. IF */
/* IT IS 10 OR 20 SET UP AN AD-HOC SHIFT. */
if (it == 30) {
goto L240;
}
if (it != 10 && it != 20) {
goto L70;
}
/* AD-HOC SHIFT. */
t += x;
i__1 = nu;
for (i__ = *nlow; i__ <= i__1; ++i__) {
a[i__ + i__ * a_dim1] -= x;
/* L60: */
}
s = (d__1 = a[nu + (nu - 1) * a_dim1], abs(d__1)) + (d__2 = a[nu - 1 + (
nu - 2) * a_dim1], abs(d__2));
x = s * .75;
y = x;
/* Computing 2nd power */
d__1 = s;
w = d__1 * d__1 * -.4375;
L70:
++it;
/* LOOK FOR TWO CONSECUTIVE SMALL SUB-DIAGONAL */
/* ELEMENTS. */
nl = nu - 2;
L80:
z__ = a[nl + nl * a_dim1];
r__ = x - z__;
s = y - z__;
p = (r__ * s - w) / a[nl + 1 + nl * a_dim1] + a[nl + (nl + 1) * a_dim1];
q = a[nl + 1 + (nl + 1) * a_dim1] - z__ - r__ - s;
r__ = a[nl + 2 + (nl + 1) * a_dim1];
s = abs(p) + abs(q) + abs(r__);
p /= s;
q /= s;
r__ /= s;
if (nl == l) {
goto L90;
}
if ((d__1 = a[nl + (nl - 1) * a_dim1], abs(d__1)) * (abs(q) + abs(r__)) <=
*eps * abs(p) * ((d__2 = a[nl - 1 + (nl - 1) * a_dim1], abs(d__2)
) + abs(z__) + (d__3 = a[nl + 1 + (nl + 1) * a_dim1], abs(d__3))))
{
goto L90;
}
--nl;
goto L80;
L90:
/* PERFORM A QR STEP BETWEEN NL AND NU. */
qrstep_(&a[a_offset], &v[v_offset], &p, &q, &r__, &nl, &nu, n, na, nv);
goto L30;
/* 2X2 BLOCK FOUND. */
L100:
if (nu != *nlow + 1) {
a[nu - 1 + (nu - 2) * a_dim1] = 0.;
}
a[nu + nu * a_dim1] += t;
a[nu - 1 + (nu - 1) * a_dim1] += t;
type__[nu] = 0;
type__[nu - 1] = 0;
mu = nu;
/* LOOP TO POSITION 2X2 BLOCK. */
L110:
nl = mu - 1;
/* ATTEMPT TO SPLIT THE BLOCK INTO TWO REAL */
/* EIGENVALUES. */
split_(&a[a_offset], &v[v_offset], n, &nl, &e1, &e2, na, nv);
/* IF THE SPLIT WAS SUCCESSFUL, GO AND ORDER THE */
/* REAL EIGENVALUES. */
if (a[mu + (mu - 1) * a_dim1] == 0.) {
goto L170;
}
/* TEST TO SEE IF THE BLOCK IS PROPERLY POSITIONED, */
/* AND IF NOT EXCHANGE IT */
if (mu == *nup) {
goto L230;
}
if (mu == *nup - 1) {
goto L130;
}
if (a[mu + 2 + (mu + 1) * a_dim1] == 0.) {
goto L130;
}
/* THE NEXT BLOCK IS 2X2. */
/* IF (A(MU-1,MU-1)*A(MU,MU)-A(MU-1,MU)*A(MU,MU-1).GE.A(MU+1, */
/* * MU+1)*A(MU+2,MU+2)-A(MU+1,MU+2)*A(MU+2,MU+1)) GO TO 230 */
if (*imfd == 1) {
if (a[mu - 1 + (mu - 1) * a_dim1] + a[mu + mu * a_dim1] >= a[mu + 1 +
(mu + 1) * a_dim1] + a[mu + 2 + (mu + 2) * a_dim1]) {
goto L230;
}
} else {
if (a[mu - 1 + (mu - 1) * a_dim1] + a[mu + mu * a_dim1] <= a[mu + 1 +
(mu + 1) * a_dim1] + a[mu + 2 + (mu + 2) * a_dim1]) {
goto L230;
}
}
exchng_(&a[a_offset], &v[v_offset], n, &nl, &c__2, &c__2, eps, &fail, na,
nv);
if (! fail) {
goto L120;
}
type__[nl] = -1;
type__[nl + 1] = -1;
type__[nl + 2] = -1;
type__[nl + 3] = -1;
goto L240;
L120:
mu += 2;
goto L150;
L130:
/* THE NEXT BLOCK IS 1X1. */
/* IF (A(MU-1,MU-1)*A(MU,MU)-A(MU-1,MU)*A(MU,MU-1).GE.A(MU+1, */
/* * MU+1)**2) GO TO 230 */
if (*imfd == 1) {
if (a[mu - 1 + (mu - 1) * a_dim1] + a[mu + mu * a_dim1] >= a[mu + 1 +
(mu + 1) * a_dim1] * 2.) {
goto L230;
}
} else {
if (a[mu - 1 + (mu - 1) * a_dim1] + a[mu + mu * a_dim1] <= a[mu + 1 +
(mu + 1) * a_dim1] * 2.) {
goto L230;
}
}
exchng_(&a[a_offset], &v[v_offset], n, &nl, &c__2, &c__1, eps, &fail, na,
nv);
if (! fail) {
goto L140;
}
type__[nl] = -1;
type__[nl + 1] = -1;
type__[nl + 2] = -1;
goto L240;
L140:
++mu;
L150:
goto L110;
/* SINGLE EIGENVALUE FOUND. */
L160:
nl = 0;
a[nu + nu * a_dim1] += t;
if (nu != *nlow) {
a[nu + (nu - 1) * a_dim1] = 0.;
}
type__[nu] = 0;
mu = nu;
/* LOOP TO POSITION ONE OR TWO REAL EIGENVALUES. */
L170:
/* POSITION THE EIGENVALUE LOCATED AT A(NL,NL). */
L180:
if (mu == *nup) {
goto L220;
}
if (mu == *nup - 1) {
goto L200;
}
if (a[mu + 2 + (mu + 1) * a_dim1] == 0.) {
goto L200;
}
/* THE NEXT BLOCK IS 2X2. */
/* IF (A(MU,MU)**2.GE.A(MU+1,MU+1)*A(MU+2,MU+2)-A(MU+1,MU+2)* */
/* * A(MU+2,MU+1)) GO TO 220 */
if (*imfd == 1) {
if (a[mu + mu * a_dim1] * 2. >= a[mu + 1 + (mu + 1) * a_dim1] + a[mu
+ 2 + (mu + 2) * a_dim1]) {
goto L220;
}
} else {
if (a[mu + mu * a_dim1] * 2. <= a[mu + 1 + (mu + 1) * a_dim1] + a[mu
+ 2 + (mu + 2) * a_dim1]) {
goto L220;
}
}
exchng_(&a[a_offset], &v[v_offset], n, &mu, &c__1, &c__2, eps, &fail, na,
nv);
if (! fail) {
goto L190;
}
type__[mu] = -1;
type__[mu + 1] = -1;
type__[mu + 2] = -1;
goto L240;
L190:
mu += 2;
goto L210;
L200:
/* THE NEXT BLOCK IS 1X1. */
/* IF (ABS(A(MU,MU)).GE.ABS(A(MU+1,MU+1))) GO TO 220 */
if (*imfd == 1) {
if (a[mu + mu * a_dim1] >= a[mu + 1 + (mu + 1) * a_dim1]) {
goto L220;
}
} else {
if (a[mu + mu * a_dim1] <= a[mu + 1 + (mu + 1) * a_dim1]) {
goto L220;
}
}
exchng_(&a[a_offset], &v[v_offset], n, &mu, &c__1, &c__1, eps, &fail, na,
nv);
++mu;
L210:
goto L180;
L220:
mu = nl;
nl = 0;
if (mu != 0) {
goto L170;
}
/* GO BACK AND GET THE NEXT EIGENVALUE. */
L230:
nu = l - 1;
goto L20;
/* ALL THE EIGENVALUES HAVE BEEN FOUND AND ORDERED. */
/* COMPUTE THEIR VALUES AND TYPE. */
L240:
if (nu < *nlow) {
goto L260;
}
i__1 = nu;
for (i__ = *nlow; i__ <= i__1; ++i__) {
a[i__ + i__ * a_dim1] += t;
/* L250: */
}
L260:
nu = *nup;
L270:
if (type__[nu] != -1) {
goto L280;
}
--nu;
goto L310;
L280:
if (nu == *nlow) {
goto L290;
}
if (a[nu + (nu - 1) * a_dim1] == 0.) {
goto L290;
}
/* 2X2 BLOCK. */
i__1 = nu - 1;
split_(&a[a_offset], &v[v_offset], n, &i__1, &e1, &e2, na, nv);
if (a[nu + (nu - 1) * a_dim1] == 0.) {
goto L290;
}
er[nu] = e1;
ei[nu - 1] = e2;
er[nu - 1] = er[nu];
ei[nu] = -ei[nu - 1];
type__[nu - 1] = 1;
type__[nu] = 2;
nu += -2;
goto L300;
L290:
/* SINGLE ROOT. */
er[nu] = a[nu + nu * a_dim1];
ei[nu] = 0.;
--nu;
L300:
L310:
if (nu >= *nlow) {
goto L270;
}
return 0;
} /* hqr3loc_ */
/* Subroutine */ int split_(a, v, n, l, e1, e2, na, nv)
doublereal *a, *v;
integer *n, *l;
doublereal *e1, *e2;
integer *na, *nv;
{
/* System generated locals */
integer a_dim1, a_offset, v_dim1, v_offset, i__1;
doublereal d__1, d__2;
/* Builtin functions */
double sqrt();
/* Local variables */
static integer i__, j;
static doublereal p, q, r__, t, u, w, x, y, z__;
static integer l1;
/* Parameter adjustments */
a_dim1 = *na;
a_offset = a_dim1 + 1;
a -= a_offset;
v_dim1 = *nv;
v_offset = v_dim1 + 1;
v -= v_offset;
/* Function Body */
x = a[*l + 1 + (*l + 1) * a_dim1];
y = a[*l + *l * a_dim1];
w = a[*l + (*l + 1) * a_dim1] * a[*l + 1 + *l * a_dim1];
p = (y - x) / 2.;
/* Computing 2nd power */
d__1 = p;
q = d__1 * d__1 + w;
if (q >= 0.) {
goto L10;
}
/* COMPLEX EIGENVALUE. */
*e1 = p + x;
*e2 = sqrt(-q);
return 0;
L10:
/* TWO REAL EIGENVALUES. SET UP TRANSFORMATION. */
z__ = sqrt(q);
if (p < 0.) {
goto L20;
}
z__ = p + z__;
goto L30;
L20:
z__ = p - z__;
L30:
if (z__ == 0.) {
goto L40;
}
r__ = -w / z__;
goto L50;
L40:
r__ = 0.;
L50:
if ((d__1 = x + z__, abs(d__1)) >= (d__2 = x + r__, abs(d__2))) {
z__ = r__;
}
y = y - x - z__;
x = -z__;
t = a[*l + (*l + 1) * a_dim1];
u = a[*l + 1 + *l * a_dim1];
if (abs(y) + abs(u) <= abs(t) + abs(x)) {
goto L60;
}
q = u;
p = y;
goto L70;
L60:
q = x;
p = t;
L70:
/* Computing 2nd power */
d__1 = p;
/* Computing 2nd power */
d__2 = q;
r__ = sqrt(d__1 * d__1 + d__2 * d__2);
if (r__ > 0.) {
goto L80;
}
*e1 = a[*l + *l * a_dim1];
*e2 = a[*l + 1 + (*l + 1) * a_dim1];
a[*l + 1 + *l * a_dim1] = 0.;
return 0;
L80:
p /= r__;
q /= r__;
/* PREMULTIPLY. */
i__1 = *n;
for (j = *l; j <= i__1; ++j) {
z__ = a[*l + j * a_dim1];
a[*l + j * a_dim1] = p * z__ + q * a[*l + 1 + j * a_dim1];
a[*l + 1 + j * a_dim1] = p * a[*l + 1 + j * a_dim1] - q * z__;
/* L90: */
}
/* POSTMULTIPLY. */
l1 = *l + 1;
i__1 = l1;
for (i__ = 1; i__ <= i__1; ++i__) {
z__ = a[i__ + *l * a_dim1];
a[i__ + *l * a_dim1] = p * z__ + q * a[i__ + (*l + 1) * a_dim1];
a[i__ + (*l + 1) * a_dim1] = p * a[i__ + (*l + 1) * a_dim1] - q * z__;
/* L100: */
}
/* ACCUMULATE THE TRANSFORMATION IN V. */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
z__ = v[i__ + *l * v_dim1];
v[i__ + *l * v_dim1] = p * z__ + q * v[i__ + (*l + 1) * v_dim1];
v[i__ + (*l + 1) * v_dim1] = p * v[i__ + (*l + 1) * v_dim1] - q * z__;
/* L110: */
}
a[*l + 1 + *l * a_dim1] = 0.;
*e1 = a[*l + *l * a_dim1];
*e2 = a[*l + 1 + (*l + 1) * a_dim1];
return 0;
} /* split_ */
/* Subroutine */ int exchng_(a, v, n, l, b1, b2, eps, fail, na, nv)
doublereal *a, *v;
integer *n, *l, *b1, *b2;
doublereal *eps;
logical *fail;
integer *na, *nv;
{
/* System generated locals */
integer a_dim1, a_offset, v_dim1, v_offset, i__1;
doublereal d__1, d__2, d__3;
/* Builtin functions */
double sqrt();
/* Local variables */
static integer i__, j, m;
static doublereal p, q, r__, s, w, x, y, z__;
static integer l1, it;
extern /* Subroutine */ int qrstep_();
/* Parameter adjustments */
a_dim1 = *na;
a_offset = a_dim1 + 1;
a -= a_offset;
v_dim1 = *nv;
v_offset = v_dim1 + 1;
v -= v_offset;
/* Function Body */
*fail = FALSE_;
if (*b1 == 2) {
goto L70;
}
if (*b2 == 2) {
goto L40;
}
/* INTERCHANGE 1X1 AND 1X1 BLOCKS. */
l1 = *l + 1;
q = a[*l + 1 + (*l + 1) * a_dim1] - a[*l + *l * a_dim1];
p = a[*l + (*l + 1) * a_dim1];
/* Computing MAX */
d__1 = abs(p), d__2 = abs(q);
r__ = max(d__1,d__2);
if (r__ == 0.) {
return 0;
}
p /= r__;
q /= r__;
/* Computing 2nd power */
d__1 = p;
/* Computing 2nd power */
d__2 = q;
r__ = sqrt(d__1 * d__1 + d__2 * d__2);
p /= r__;
q /= r__;
i__1 = *n;
for (j = *l; j <= i__1; ++j) {
s = p * a[*l + j * a_dim1] + q * a[*l + 1 + j * a_dim1];
a[*l + 1 + j * a_dim1] = p * a[*l + 1 + j * a_dim1] - q * a[*l + j *
a_dim1];
a[*l + j * a_dim1] = s;
/* L10: */
}
i__1 = l1;
for (i__ = 1; i__ <= i__1; ++i__) {
s = p * a[i__ + *l * a_dim1] + q * a[i__ + (*l + 1) * a_dim1];
a[i__ + (*l + 1) * a_dim1] = p * a[i__ + (*l + 1) * a_dim1] - q * a[
i__ + *l * a_dim1];
a[i__ + *l * a_dim1] = s;
/* L20: */
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
s = p * v[i__ + *l * v_dim1] + q * v[i__ + (*l + 1) * v_dim1];
v[i__ + (*l + 1) * v_dim1] = p * v[i__ + (*l + 1) * v_dim1] - q * v[
i__ + *l * v_dim1];
v[i__ + *l * v_dim1] = s;
/* L30: */
}
a[*l + 1 + *l * a_dim1] = 0.;
return 0;
L40:
/* INTERCHANGE 1X1 AND 2X2 BLOCKS. */
x = a[*l + *l * a_dim1];
p = 1.;
q = 1.;
r__ = 1.;
i__1 = *l + 2;
qrstep_(&a[a_offset], &v[v_offset], &p, &q, &r__, l, &i__1, n, na, nv);
it = 0;
L50:
++it;
if (it <= 30) {
goto L60;
}
*fail = TRUE_;
return 0;
L60:
p = a[*l + *l * a_dim1] - x;
q = a[*l + 1 + *l * a_dim1];
r__ = 0.;
i__1 = *l + 2;
qrstep_(&a[a_offset], &v[v_offset], &p, &q, &r__, l, &i__1, n, na, nv);
if ((d__1 = a[*l + 2 + (*l + 1) * a_dim1], abs(d__1)) > *eps * ((d__2 = a[
*l + 1 + (*l + 1) * a_dim1], abs(d__2)) + (d__3 = a[*l + 2 + (*l
+ 2) * a_dim1], abs(d__3)))) {
goto L50;
}
a[*l + 2 + (*l + 1) * a_dim1] = 0.;
return 0;
L70:
/* INTERCHANGE 2X2 AND B2XB2 BLOCKS. */
m = *l + 2;
if (*b2 == 2) {
++m;
}
x = a[*l + 1 + (*l + 1) * a_dim1];
y = a[*l + *l * a_dim1];
w = a[*l + 1 + *l * a_dim1] * a[*l + (*l + 1) * a_dim1];
p = 1.;
q = 1.;
r__ = 1.;
qrstep_(&a[a_offset], &v[v_offset], &p, &q, &r__, l, &m, n, na, nv);
it = 0;
L80:
++it;
if (it <= 30) {
goto L90;
}
*fail = TRUE_;
return 0;
L90:
z__ = a[*l + *l * a_dim1];
r__ = x - z__;
s = y - z__;
p = (r__ * s - w) / a[*l + 1 + *l * a_dim1] + a[*l + (*l + 1) * a_dim1];
q = a[*l + 1 + (*l + 1) * a_dim1] - z__ - r__ - s;
r__ = a[*l + 2 + (*l + 1) * a_dim1];
s = abs(p) + abs(q) + abs(r__);
p /= s;
q /= s;
r__ /= s;
qrstep_(&a[a_offset], &v[v_offset], &p, &q, &r__, l, &m, n, na, nv);
if ((d__1 = a[m - 1 + (m - 2) * a_dim1], abs(d__1)) > *eps * ((d__2 = a[m
- 1 + (m - 1) * a_dim1], abs(d__2)) + (d__3 = a[m - 2 + (m - 2) *
a_dim1], abs(d__3)))) {
goto L80;
}
a[m - 1 + (m - 2) * a_dim1] = 0.;
return 0;
} /* exchng_ */
/* Subroutine */ int qrstep_(a, v, p, q, r__, nl, nu, n, na, nv)
doublereal *a, *v, *p, *q, *r__;
integer *nl, *nu, *n, *na, *nv;
{
/* System generated locals */
integer a_dim1, a_offset, v_dim1, v_offset, i__1, i__2;
doublereal d__1, d__2, d__3;
/* Builtin functions */
double sqrt();
/* Local variables */
static logical last;
static integer i__, j, k;
static doublereal s, x, y, z__;
static integer nl2, nl3, num1;
/* INTERNAL VARIABLES. */
/* Parameter adjustments */
a_dim1 = *na;
a_offset = a_dim1 + 1;
a -= a_offset;
v_dim1 = *nv;
v_offset = v_dim1 + 1;
v -= v_offset;
/* Function Body */
nl2 = *nl + 2;
i__1 = *nu;
for (i__ = nl2; i__ <= i__1; ++i__) {
a[i__ + (i__ - 2) * a_dim1] = 0.;
/* L10: */
}
if (nl2 == *nu) {
goto L30;
}
nl3 = *nl + 3;
i__1 = *nu;
for (i__ = nl3; i__ <= i__1; ++i__) {
a[i__ + (i__ - 3) * a_dim1] = 0.;
/* L20: */
}
L30:
num1 = *nu - 1;
i__1 = num1;
for (k = *nl; k <= i__1; ++k) {
/* DETERMINE THE TRANSFORMATION. */
last = k == num1;
if (k == *nl) {
goto L40;
}
*p = a[k + (k - 1) * a_dim1];
*q = a[k + 1 + (k - 1) * a_dim1];
*r__ = 0.;
if (! last) {
*r__ = a[k + 2 + (k - 1) * a_dim1];
}
x = abs(*p) + abs(*q) + abs(*r__);
if (x == 0.) {
goto L130;
}
*p /= x;
*q /= x;
*r__ /= x;
L40:
/* Computing 2nd power */
d__1 = *p;
/* Computing 2nd power */
d__2 = *q;
/* Computing 2nd power */
d__3 = *r__;
s = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
if (*p < 0.) {
s = -s;
}
if (k == *nl) {
goto L50;
}
a[k + (k - 1) * a_dim1] = -s * x;
goto L60;
L50:
if (*nl != 1) {
a[k + (k - 1) * a_dim1] = -a[k + (k - 1) * a_dim1];
}
L60:
*p += s;
x = *p / s;
y = *q / s;
z__ = *r__ / s;
*q /= *p;
*r__ /= *p;
/* PREMULTIPLY. */
i__2 = *n;
for (j = k; j <= i__2; ++j) {
*p = a[k + j * a_dim1] + *q * a[k + 1 + j * a_dim1];
if (last) {
goto L70;
}
*p += *r__ * a[k + 2 + j * a_dim1];
a[k + 2 + j * a_dim1] -= *p * z__;
L70:
a[k + 1 + j * a_dim1] -= *p * y;
a[k + j * a_dim1] -= *p * x;
/* L80: */
}
/* POSTMULTIPLY. */
/* Computing MIN */
i__2 = k + 3;
j = min(i__2,*nu);
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__) {
*p = x * a[i__ + k * a_dim1] + y * a[i__ + (k + 1) * a_dim1];
if (last) {
goto L90;
}
*p += z__ * a[i__ + (k + 2) * a_dim1];
a[i__ + (k + 2) * a_dim1] -= *p * *r__;
L90:
a[i__ + (k + 1) * a_dim1] -= *p * *q;
a[i__ + k * a_dim1] -= *p;
/* L100: */
}
/* ACCUMULATE THE TRANSFORMATION IN V. */
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
*p = x * v[i__ + k * v_dim1] + y * v[i__ + (k + 1) * v_dim1];
if (last) {
goto L110;
}
*p += z__ * v[i__ + (k + 2) * v_dim1];
v[i__ + (k + 2) * v_dim1] -= *p * *r__;
L110:
v[i__ + (k + 1) * v_dim1] -= *p * *q;
v[i__ + k * v_dim1] -= *p;
/* L120: */
}
L130:
;
}
return 0;
} /* qrstep_ */
/* Subroutine */ int orthes_(nm, n, low, igh, a, ort)
integer *nm, *n, *low, *igh;
doublereal *a, *ort;
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3;
doublereal d__1;
/* Builtin functions */
double sqrt(), d_sign();
/* Local variables */
static doublereal f, g, h__;
static integer i__, j, m;
static doublereal scale;
static integer la, ii, jj, mp, kp1;
/* Parameter adjustments */
a_dim1 = *nm;
a_offset = a_dim1 + 1;
a -= a_offset;
--ort;
/* Function Body */
la = *igh - 1;
kp1 = *low + 1;
if (la < kp1) {
goto L200;
}
i__1 = la;
for (m = kp1; m <= i__1; ++m) {
h__ = 0.;
ort[m] = 0.;
scale = 0.;
/* .......... scale column (algol tol then not needed) ..........
*/
i__2 = *igh;
for (i__ = m; i__ <= i__2; ++i__) {
/* L90: */
scale += (d__1 = a[i__ + (m - 1) * a_dim1], abs(d__1));
}
if (scale == 0.) {
goto L180;
}
mp = m + *igh;
/* .......... for i=igh step -1 until m do -- .......... */
i__2 = *igh;
for (ii = m; ii <= i__2; ++ii) {
i__ = mp - ii;
ort[i__] = a[i__ + (m - 1) * a_dim1] / scale;
h__ += ort[i__] * ort[i__];
/* L100: */
}
d__1 = sqrt(h__);
g = -d_sign(&d__1, &ort[m]);
h__ -= ort[m] * g;
ort[m] -= g;
/* .......... form (i-(u*ut)/h) * a .......... */
i__2 = *n;
for (j = m; j <= i__2; ++j) {
f = 0.;
/* .......... for i=igh step -1 until m do -- .......... */
i__3 = *igh;
for (ii = m; ii <= i__3; ++ii) {
i__ = mp - ii;
f += ort[i__] * a[i__ + j * a_dim1];
/* L110: */
}
f /= h__;
i__3 = *igh;
for (i__ = m; i__ <= i__3; ++i__) {
/* L120: */
a[i__ + j * a_dim1] -= f * ort[i__];
}
/* L130: */
}
/* .......... form (i-(u*ut)/h)*a*(i-(u*ut)/h) .......... */
i__2 = *igh;
for (i__ = 1; i__ <= i__2; ++i__) {
f = 0.;
/* .......... for j=igh step -1 until m do -- .......... */
i__3 = *igh;
for (jj = m; jj <= i__3; ++jj) {
j = mp - jj;
f += ort[j] * a[i__ + j * a_dim1];
/* L140: */
}
f /= h__;
i__3 = *igh;
for (j = m; j <= i__3; ++j) {
/* L150: */
a[i__ + j * a_dim1] -= f * ort[j];
}
/* L160: */
}
ort[m] = scale * ort[m];
a[m + (m - 1) * a_dim1] = scale * g;
L180:
;
}
L200:
return 0;
} /* orthes_ */
/* Subroutine */ int ortran_(nm, n, low, igh, a, ort, z__)
integer *nm, *n, *low, *igh;
doublereal *a, *ort, *z__;
{
/* System generated locals */
integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3;
/* Local variables */
static doublereal g;
static integer i__, j, kl, mm, mp, mp1;
/* Parameter adjustments */
z_dim1 = *nm;
z_offset = z_dim1 + 1;
z__ -= z_offset;
--ort;
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;
if (a[mp + (mp - 1) * a_dim1] == 0.) {
goto L140;
}
mp1 = mp + 1;
i__2 = *igh;
for (i__ = mp1; i__ <= i__2; ++i__) {
/* L100: */
ort[i__] = a[i__ + (mp - 1) * a_dim1];
}
i__2 = *igh;
for (j = mp; j <= i__2; ++j) {
g = 0.;
i__3 = *igh;
for (i__ = mp; i__ <= i__3; ++i__) {
/* L110: */
g += ort[i__] * z__[i__ + j * z_dim1];
}
/* .......... divisor below is negative of h formed in orthes.
*/
/* double division avoids possible underflow ......
.... */
g = g / ort[mp] / a[mp + (mp - 1) * a_dim1];
i__3 = *igh;
for (i__ = mp; i__ <= i__3; ++i__) {
/* L120: */
z__[i__ + j * z_dim1] += g * ort[i__];
}
/* L130: */
}
L140:
;
}
L200:
return 0;
} /* ortran_ */
syntax highlighted by Code2HTML, v. 0.9.1