/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:19 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "sevun.h"
#include <stdlib.h>
void /*FUNCTION*/ sevun(
float *a,
long lda,
long n,
float evalr[],
float evali[],
long iflag[])
{
#define A(I_,J_) (*(a+(I_)*(lda)+(J_)))
LOGICAL32 notlas;
long int en, enm2, i, igh, itn, its, j, k, l, low, ltype, m, mp2,
na;
float norm, p, q, r, s, t, tst1, tst2, w, x, y, zz;
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const Evali = &evali[0] - 1;
float *const Evalr = &evalr[0] - 1;
long *const Iflag = &iflag[0] - 1;
/* end of OFFSET VECTORS */
/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
* ALL RIGHTS RESERVED.
* Based on Government Sponsored Research NAS7-03001.
*>> 1996-03-30 SEVUN Krogh MIN0 => MIN
*>> 1994-10-20 SEVUN Krogh Changes to use M77CON
*>> 1992-04-24 SEVUN CLL Minor edits.
*>> 1992-04-23 SEVUN Krogh Made DP version compatible with SP version.
*>> 1991-10-25 SEVUN Krogh Initial version, converted from EISPACK.
*
* This subroutine uses slight modifcations of EISPACK routines
* BALANC, ELMHES, and HQR to get the eigenvalues of a general real
* matrix by the QR method. The first two are are encapsulated in
* SEVBH. HQR, which is inline here, is a translation of the ALGOL
* procedure HQR, Num. Math. 14, 219-231(1970) by Martin, Peters, and
* Wilkinson.
* Handbook for Auto. Comp., Vol.ii-Linear Algebra, 359-371(1971).
*
* On input
* A contains the input matrix whose eigenvalues are desired.
* LDA 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.
*
* On output
* A has been destroyed.
* EVALR and EVALI contain real and imaginary parts, respectively, of
* the eigenvalues. The eigenvalues are given in order of
* increasing real parts. When real parts are equal they are
* given in order of increasing absolute complex part. Complex
* conjugate pairs of values appear consecutively with
* the eigenvalue having the positive imaginary part first. If
* an error exit is made, the eigenvalues should be correct
* (but unordered) for indices IFLAG(1)+1,...,N.
* IFLAG(1) is set to
* 1 If all eigenvalues are real.
* 2 If some eigenvalues are complex.
* 3 If N < 1 on the initial entry.
* 4 If the limit of 30*N iterations is exhausted.
* ------------------------------------------------------------------
*--S replaces "?": ?EVUN, ?EVBH
* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------
* */
sevbh( a, lda, n, &low, &igh, iflag, evalr );
ltype = 1;
norm = 0.0e0;
k = 1;
/* .......... store roots isolated by SEVBH
* and compute matrix norm .......... */
for (i = 1; i <= n; i++)
{
for (j = k; j <= n; j++)
{
norm += fabsf( A(j - 1,i - 1) );
}
k = i;
if ((i < low) || (i > igh))
{
Evalr[i] = A(i - 1,i - 1);
Evali[i] = 0.0e0;
}
}
en = igh;
t = 0.0e0;
itn = 30*n;
/* .......... search for next eigenvalues .......... */
L_60:
if (en < low)
goto L_300;
its = 0;
na = en - 1;
enm2 = na - 1;
/* .......... look for single small sub-diagonal element */
L_70:
for (l = en; l >= (low + 1); l--)
{
s = fabsf( A(l - 2,l - 2) ) + fabsf( A(l - 1,l - 1) );
if (s == 0.0e0)
s = norm;
tst1 = s;
tst2 = tst1 + fabsf( A(l - 2,l - 1) );
if (tst2 == tst1)
goto L_100;
}
/* .......... form shift .......... */
L_100:
x = A(en - 1,en - 1);
if (l == en)
goto L_270;
y = A(na - 1,na - 1);
w = A(na - 1,en - 1)*A(en - 1,na - 1);
if (l == na)
goto L_280;
if (itn <= 0)
{
/* .......... set error -- all eigenvalues have not
* converged after 30*N iterations .......... */
ermsg( "SEVUN", en, 0, "ERROR NO. is index of eigenvalue causing convergence failure."
, '.' );
Iflag[1] = 4;
if (en <= 0)
Iflag[1] = 3;
return;
}
if (its != 10 && its != 20)
goto L_130;
/* .......... form exceptional shift .......... */
t += x;
for (i = low; i <= en; i++)
{
A(i - 1,i - 1) -= x;
}
s = fabsf( A(na - 1,en - 1) ) + fabsf( A(enm2 - 1,na - 1) );
x = 0.75e0*s;
y = x;
w = -0.4375e0*s*s;
L_130:
its += 1;
itn -= 1;
/* .......... look for two consecutive small sub-diagonal elements. */
for (m = enm2; m >= l; m--)
{
zz = A(m - 1,m - 1);
r = x - zz;
s = y - zz;
p = (r*s - w)/A(m - 1,m) + A(m,m - 1);
q = A(m,m) - zz - r - s;
r = A(m,m + 1);
s = fabsf( p ) + fabsf( q ) + fabsf( r );
p /= s;
q /= s;
r /= s;
if (m == l)
goto L_150;
tst1 = fabsf( p )*(fabsf( A(m - 2,m - 2) ) + fabsf( zz ) + fabsf( A(m,m) ));
tst2 = tst1 + fabsf( A(m - 2,m - 1) )*(fabsf( q ) + fabsf( r ));
if (tst2 == tst1)
goto L_150;
}
L_150:
mp2 = m + 2;
for (i = mp2; i <= en; i++)
{
A(i - 3,i - 1) = 0.0e0;
if (i != mp2)
A(i - 4,i - 1) = 0.0e0;
}
/* .......... double QR step involving rows L to EN and
* columns M to EN .......... */
for (k = m; k <= na; k++)
{
notlas = k != na;
if (k != m)
{
p = A(k - 2,k - 1);
q = A(k - 2,k);
r = 0.0e0;
if (notlas)
r = A(k - 2,k + 1);
x = fabsf( p ) + fabsf( q ) + fabsf( r );
if (x == 0.0e0)
goto L_260;
p /= x;
q /= x;
r /= x;
}
s = signf( sqrtf( p*p + q*q + r*r ), p );
if (k != m)
{
A(k - 2,k - 1) = -s*x;
}
else
{
if (l != m)
A(k - 2,k - 1) = -A(k - 2,k - 1);
}
p += s;
x = p/s;
y = q/s;
zz = r/s;
q /= p;
r /= p;
if (notlas)
{
/* .......... row modification .......... */
for (j = k; j <= n; j++)
{
p = A(j - 1,k - 1) + q*A(j - 1,k) + r*A(j - 1,k + 1);
A(j - 1,k - 1) += -p*x;
A(j - 1,k) += -p*y;
A(j - 1,k + 1) += -p*zz;
}
j = min( en, k + 3 );
/* .......... column modification .......... */
for (i = 1; i <= j; i++)
{
p = x*A(k - 1,i - 1) + y*A(k,i - 1) + zz*A(k + 1,i - 1);
A(k - 1,i - 1) -= p;
A(k,i - 1) += -p*q;
A(k + 1,i - 1) += -p*r;
}
}
else
{
/* .......... row modification .......... */
for (j = k; j <= n; j++)
{
p = A(j - 1,k - 1) + q*A(j - 1,k);
A(j - 1,k - 1) += -p*x;
A(j - 1,k) += -p*y;
}
j = min( en, k + 3 );
/* .......... column modification .......... */
for (i = 1; i <= j; i++)
{
p = x*A(k - 1,i - 1) + y*A(k,i - 1);
A(k - 1,i - 1) -= p;
A(k,i - 1) += -p*q;
}
}
L_260:
;
}
goto L_70;
/* .......... one root found .......... */
L_270:
Evalr[en] = x + t;
Evali[en] = 0.0e0;
en = na;
goto L_60;
/* .......... two roots found .......... */
L_280:
p = (y - x)/2.0e0;
q = p*p + w;
zz = sqrtf( fabsf( q ) );
x += t;
if (q >= 0.0e0)
{
/* .......... real pair .......... */
zz = p + signf( zz, p );
Evalr[na] = x + zz;
Evalr[en] = Evalr[na];
if (zz != 0.0e0)
Evalr[en] = x - w/zz;
Evali[na] = 0.0e0;
Evali[en] = 0.0e0;
}
else
{
/* .......... complex pair .......... */
ltype = 2;
Evalr[na] = x + p;
Evalr[en] = x + p;
Evali[na] = zz;
Evali[en] = -zz;
}
en = enm2;
goto L_60;
L_300:
;
/*-- Begin mask code changes
* Set up for Shell sort
* Sort so real parts are algebraically increasing
* For = real parts, so abs(imag. parts) are increasing
* For both =, sort on index -- preserves complex pair order
*-- End mask code changes */
for (i = 1; i <= n; i++)
{
Iflag[i] = i;
}
l = 1;
for (k = 1; k <= n; k++)
{
l = 3*l + 1;
if (l >= n)
goto L_2020;
}
L_2020:
l = max( 1, (l - 4)/9 );
L_2030:
for (j = l + 1; j <= n; j++)
{
k = Iflag[j];
p = Evalr[k];
i = j - l;
L_2040:
switch (SARITHIF(p - Evalr[Iflag[i]]))
{
case -1: goto L_2070;
case 0: goto L_2050;
case 1: goto L_2080;
}
L_2050:
switch (SARITHIF(fabsf( Evali[k] ) - fabsf( Evali[Iflag[i]] )))
{
case -1: goto L_2070;
case 0: goto L_2060;
case 1: goto L_2080;
}
L_2060:
if (k > Iflag[i])
goto L_2080;
L_2070:
Iflag[i + l] = Iflag[i];
i -= l;
if (i > 0)
goto L_2040;
L_2080:
Iflag[i + l] = k;
}
l = (l - 1)/3;
if (l != 0)
goto L_2030;
/* Indices in IFLAG now give the desired order --
* Move entries to get this order. */
L_2110:
for (i = l + 1; i <= n; i++)
{
if (Iflag[i] != i)
{
l = i;
m = i;
p = Evalr[i];
q = Evali[i];
L_2120:
k = Iflag[m];
Iflag[m] = m;
if (k != l)
{
Evalr[m] = Evalr[k];
Evali[m] = Evali[k];
m = k;
goto L_2120;
}
else
{
Evalr[m] = p;
Evali[m] = q;
goto L_2110;
}
}
}
Iflag[1] = ltype;
return;
#undef A
} /* end of function */