C#，数值计算，计算实非对称矩阵的所有特征值和特征向量，简化为Hes-senberg形式，然后进行QR迭代
1 文本格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// Computes all eigenvalues and eigenvectors of a real nonsymmetric matrix by
/// reduction to Hes- senberg form followed by QR iteration. - PTC
/// </summary>
public class Unsymmeig
{
private int n { get; set; }
private double[,] a { get; set; }
private double[,] zz { get; set; }
public Complex[] wri { get; set; }
private double[] scale { get; set; }
private int[] perm { get; set; }
private bool yesvecs { get; set; }
private bool hessen { get; set; }
/// <summary>
/// Computes all eigenvalues and(optionally) eigenvectors of a real
/// nonsymmetric matrix a[0..n - 1][0..n - 1] by reduction to Hessenberg form
/// followed by QR iteration.If yesvecs is input as true (the default), then
/// the eigenvectors are computed.Otherwise, only the eigenvalues are
/// computed.If hessen is input as false (the default), the matrix is first
/// reduced to Hessenberg form.Otherwise it is assumed that the matrix is
/// already in Hessen- berg from.On output, wri[0..n - 1] contains the
/// eigenvalues of a sorted into descending order, while zz[0..n - 1][0..n - 1] is
/// a matrix whose columns contain the corresponding eigenvectors. For a
/// complex eigenvalue, only the eigenvector corresponding to the eigen- value
/// with a positive imaginary part is stored, with the real part in
/// zz[0..n-1][i] and the imaginary part in h.zz[0..n - 1][i + 1]. The eigenvectors
/// are not normalized.
/// </summary>
/// <param name="aa"></param>
/// <param name="yesvec"></param>
/// <param name="hessenb"></param>
public Unsymmeig(double[,] aa, bool yesvec = true, bool hessenb = false)
{
this.n = aa.GetLength(0);
this.a = aa;
this.zz = new double[n, n];
this.wri = new Complex[n];
this.scale = new double[n];
this.perm = new int[n];
this.yesvecs = yesvec;
this.hessen = hessenb;
for (int i = 0; i < n; i++)
{
scale[i] = 1.0;
}
balance();
if (!hessen)
{
elmhes();
}
if (yesvecs)
{
for (int i = 0; i < n; i++)
{
zz[i, i] = 1.0;
}
if (!hessen)
{
eltran();
}
hqr2();
balbak();
sortvecs();
}
else
{
hqr();
sort();
}
}
/// <summary>
/// Given a matrix a[0..n - 1][0..n-1], this routine replaces it by a balanced
/// matrix with identical eigenvalues.A symmetric matrix is already balanced
/// and is unaffected by this procedure.
/// </summary>
public void balance()
{
// https://en.cppreference.com/w/cpp/types/climits
double RADIX = 2.0;// numeric_limits<double>.radix;
bool done = false;
double sqrdx = RADIX * RADIX;
while (!done)
{
done = true;
for (int i = 0; i < n; i++)
{
double r = 0.0;
double c = 0.0;
for (int j = 0; j < n; j++)
{
if (j != i)
{
c += Math.Abs(a[j, i]);
r += Math.Abs(a[i, j]);
}
}
if (c != 0.0 && r != 0.0)
{
double g = r / RADIX;
double f = 1.0;
double s = c + r;
while (c < g)
{
f *= RADIX;
c *= sqrdx;
}
g = r * RADIX;
while (c > g)
{
f /= RADIX;
c /= sqrdx;
}
if ((c + r) / f < 0.95 * s)
{
done = false;
g = 1.0 / f;
scale[i] *= f;
for (int j = 0; j < n; j++)
{
a[i, j] *= g;
}
for (int j = 0; j < n; j++)
{
a[j, i] *= f;
}
}
}
}
}
}
/// <summary>
/// Forms the eigenvectors of a real nonsymmetric matrix by back transforming
/// those of the corre- sponding balanced matrix determined by balance.
/// </summary>
public void balbak()
{
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
zz[i, j] *= scale[i];
}
}
}
/// <summary>
/// Reduction to Hessenberg form by the elimination method.Replaces the real
/// nonsymmetric matrix a[0..n - 1][0..n - 1] by an upper Hessenberg matrix with
/// identical eigenvalues.Rec- ommended, but not required, is that this
/// routine be preceded by balance. On output, the Hessenberg matrix is in
/// elements a[i][j] with i <= j+1. Elements with i > j+1 are to be thought of
/// as zero, but are returned with random values.
/// </summary>
public void elmhes()
{
for (int m = 1; m < n - 1; m++)
{
double x = 0.0;
int i = m;
for (int j = m; j < n; j++)
{
if (Math.Abs(a[j, m - 1]) > Math.Abs(x))
{
x = a[j, m - 1];
i = j;
}
}
perm[m] = i;
if (i != m)
{
for (int j = m - 1; j < n; j++)
{
Globals.SWAP(ref a[i, j], ref a[m, j]);
}
for (int j = 0; j < n; j++)
{
Globals.SWAP(ref a[j, i], ref a[j, m]);
}
}
if (x != 0.0)
{
for (i = m + 1; i < n; i++)
{
double y = a[i, m - 1];
if (y != 0.0)
{
y /= x;
a[i, m - 1] = y;
for (int j = m; j < n; j++)
{
a[i, j] -= y * a[m, j];
}
for (int j = 0; j < n; j++)
{
a[j, m] += y * a[j, i];
}
}
}
}
}
}
/// <summary>
/// This routine accumulates the stabilized elementary similarity
/// transformations used in the reduction to upper Hessenberg form by elmhes.
/// The multipliers that were used in the reduction are obtained from the lower
/// triangle (below the subdiagonal) of a.The transformations are permuted
/// according to the permutations stored in perm by elmhes.
/// </summary>
public void eltran()
{
for (int mp = n - 2; mp > 0; mp--)
{
for (int k = mp + 1; k < n; k++)
{
zz[k, mp] = a[k, mp - 1];
}
int i = perm[mp];
if (i != mp)
{
for (int j = mp; j < n; j++)
{
zz[mp, j] = zz[i, j];
zz[i, j] = 0.0;
}
zz[i, mp] = 1.0;
}
}
}
public void hqr()
{
double r = 0.0;
double q = 0.0;
double p = 0.0;
double anorm = 0.0;
const double EPS = float.Epsilon;
for (int i = 0; i < n; i++)
{
for (int j = Math.Max(i - 1, 0); j < n; j++)
{
anorm += Math.Abs(a[i, j]);
}
}
int nn = n - 1;
double t = 0.0;
while (nn >= 0)
{
int its = 0;
int l;
do
{
for (l = nn; l > 0; l--)
{
double s = Math.Abs(a[l - 1, l - 1]) + Math.Abs(a[l, l]);
//if (s == 0.0)
if (s <= float.Epsilon)
{
s = anorm;
}
if (Math.Abs(a[l, l - 1]) <= EPS * s)
{
a[l, l - 1] = 0.0;
break;
}
}
double x = a[nn, nn];
if (l == nn)
{
wri[nn--] = new Complex(x + t);
}
else
{
double y = a[nn - 1, nn - 1];
double w = a[nn, nn - 1] * a[nn - 1, nn];
if (l == nn - 1)
{
p = 0.5 * (y - x);
q = p * p + w;
double z = Math.Sqrt(Math.Abs(q));
x += t;
if (q >= 0.0)
{
z = p + Globals.SIGN(z, p);
wri[nn - 1] = wri[nn] = new Complex(x + z);
if (z != 0.0)
{
wri[nn] = new Complex(x - w / z);
}
}
else
{
wri[nn] = new Complex(x + p, -z);
wri[nn - 1] = (wri[nn].conj());
}
nn -= 2;
}
else
{
if (its == 30)
{
throw new Exception("Too many iterations in hqr");
}
if (its == 10 || its == 20)
{
t += x;
for (int i = 0; i < nn + 1; i++)
{
a[i, i] -= x;
}
double s = Math.Abs(a[nn, nn - 1]) + Math.Abs(a[nn - 1, nn - 2]);
y = x = 0.75 * s;
w = -0.4375 * s * s;
}
++its;
int m = nn - 2;
for (; m >= l; m--)
{
double z = a[m, m];
r = x - z;
double s = y - z;
p = (r * s - w) / a[m + 1, m] + a[m, m + 1];
q = a[m + 1, m + 1] - z - r - s;
r = a[m + 2, m + 1];
s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);
p /= s;
q /= s;
r /= s;
if (m == l)
{
break;
}
double u = Math.Abs(a[m, m - 1]) * (Math.Abs(q) + Math.Abs(r));
double v = Math.Abs(p) * (Math.Abs(a[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(a[m + 1, m + 1]));
if (u <= EPS * v)
{
break;
}
}
for (int i = m; i < nn - 1; i++)
{
a[i + 2, i] = 0.0;
if (i != m) a[i + 2, i - 1] = 0.0;
}
for (int k = m; k < nn; k++)
{
if (k != m)
{
p = a[k, k - 1];
q = a[k + 1, k - 1];
r = 0.0;
if (k + 1 != nn) r = a[k + 2, k - 1];
if ((x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r)) != 0.0)
{
p /= x;
q /= x;
r /= x;
}
}
double s = Globals.SIGN(Math.Sqrt(p * p + q * q + r * r), p);
if ((s) != 0.0)
{
if (k == m)
{
if (l != m)
{
a[k, k - 1] = -a[k, k - 1];
}
}
else
{
a[k, k - 1] = -s * x;
}
p += s;
x = p / s;
y = q / s;
double z = r / s;
q /= p;
r /= p;
for (int j = k; j < nn + 1; j++)
{
p = a[k, j] + q * a[k + 1, j];
if (k + 1 != nn)
{
p += r * a[k + 2, j];
a[k + 2, j] -= p * z;
}
a[k + 1, j] -= p * y;
a[k, j] -= p * x;
}
int mmin = nn < k + 3 ? nn : k + 3;
for (int i = l; i < mmin + 1; i++)
{
p = x * a[i, k] + y * a[i, k + 1];
if (k + 1 != nn)
{
p += z * a[i, k + 2];
a[i, k + 2] -= p * r;
}
a[i, k + 1] -= p * q;
a[i, k] -= p;
}
}
}
}
}
} while (l + 1 < nn);
}
}
public void hqr2()
{
double z = 0.0;
double r = 0.0;
double q = 0.0;
double p = 0.0;
double anorm = 0.0;
double EPS = float.Epsilon;
for (int i = 0; i < n; i++)
{
for (int j = Math.Max(i - 1, 0); j < n; j++)
{
anorm += Math.Abs(a[i, j]);
}
}
int nn = n - 1;
double t = 0.0;
while (nn >= 0)
{
int its = 0;
int l;
do
{
for (l = nn; l > 0; l--)
{
double s = Math.Abs(a[l - 1, l - 1]) + Math.Abs(a[l, l]);
//if (s == 0.0)
if (s <= float.Epsilon)
{
s = anorm;
}
if (Math.Abs(a[l, l - 1]) <= EPS * s)
{
a[l, l - 1] = 0.0;
break;
}
}
double x = a[nn, nn];
if (l == nn)
{
a[nn, nn] = (x + t);
wri[nn] = new Complex(a[nn, nn]);
nn--;
}
else
{
double y = a[nn - 1, nn - 1];
double w = a[nn, nn - 1] * a[nn - 1, nn];
if (l == nn - 1)
{
p = 0.5 * (y - x);
q = p * p + w;
z = Math.Sqrt(Math.Abs(q));
x += t;
a[nn, nn] = x;
a[nn - 1, nn - 1] = y + t;
if (q >= 0.0)
{
z = p + Globals.SIGN(z, p);
wri[nn - 1] = wri[nn] = new Complex(x + z);
if (z != 0.0) wri[nn] = new Complex(x - w / z);
x = a[nn, nn - 1];
double s = Math.Abs(x) + Math.Abs(z);
p = x / s;
q = z / s;
r = Math.Sqrt(p * p + q * q);
p /= r;
q /= r;
for (int j = nn - 1; j < n; j++)
{
z = a[nn - 1, j];
a[nn - 1, j] = q * z + p * a[nn, j];
a[nn, j] = q * a[nn, j] - p * z;
}
for (int i = 0; i <= nn; i++)
{
z = a[i, nn - 1];
a[i, nn - 1] = q * z + p * a[i, nn];
a[i, nn] = q * a[i, nn] - p * z;
}
for (int i = 0; i < n; i++)
{
z = zz[i, nn - 1];
zz[i, nn - 1] = q * z + p * zz[i, nn];
zz[i, nn] = q * zz[i, nn] - p * z;
}
}
else
{
wri[nn] = new Complex(x + p, -z);
wri[nn - 1] = wri[nn].conj();
}
nn -= 2;
}
else
{
if (its == 30) throw new Exception("Too many iterations in hqr");
if (its == 10 || its == 20)
{
t += x;
for (int i = 0; i < nn + 1; i++)
{
a[i, i] -= x;
}
double s = Math.Abs(a[nn, nn - 1]) + Math.Abs(a[nn - 1, nn - 2]);
y = x = 0.75 * s;
w = -0.4375 * s * s;
}
++its;
int m = nn - 2;
for (; m >= l; m--)
{
z = a[m, m];
r = x - z;
double s = y - z;
p = (r * s - w) / a[m + 1, m] + a[m, m + 1];
q = a[m + 1, m + 1] - z - r - s;
r = a[m + 2, m + 1];
s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);
p /= s;
q /= s;
r /= s;
if (m == l)
{
break;
}
double u = Math.Abs(a[m, m - 1]) * (Math.Abs(q) + Math.Abs(r));
double v = Math.Abs(p) * (Math.Abs(a[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(a[m + 1, m + 1]));
if (u <= EPS * v)
{
break;
}
}
for (int i = m; i < nn - 1; i++)
{
a[i + 2, i] = 0.0;
if (i != m)
{
a[i + 2, i - 1] = 0.0;
}
}
for (int k = m; k < nn; k++)
{
if (k != m)
{
p = a[k, k - 1];
q = a[k + 1, k - 1];
r = 0.0;
if (k + 1 != nn) r = a[k + 2, k - 1];
if ((x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r)) != 0.0)
{
p /= x;
q /= x;
r /= x;
}
}
double s = Globals.SIGN(Math.Sqrt(p * p + q * q + r * r), p);
if ((s) != 0.0)
{
if (k == m)
{
if (l != m)
{
a[k, k - 1] = -a[k, k - 1];
}
}
else
{
a[k, k - 1] = -s * x;
}
p += s;
x = p / s;
y = q / s;
z = r / s;
q /= p;
r /= p;
for (int j = k; j < n; j++)
{
p = a[k, j] + q * a[k + 1, j];
if (k + 1 != nn)
{
p += r * a[k + 2, j];
a[k + 2, j] -= p * z;
}
a[k + 1, j] -= p * y;
a[k, j] -= p * x;
}
int mmin = nn < k + 3 ? nn : k + 3;
for (int i = 0; i < mmin + 1; i++)
{
p = x * a[i, k] + y * a[i, k + 1];
if (k + 1 != nn)
{
p += z * a[i, k + 2];
a[i, k + 2] -= p * r;
}
a[i, k + 1] -= p * q;
a[i, k] -= p;
}
for (int i = 0; i < n; i++)
{
p = x * zz[i, k] + y * zz[i, k + 1];
if (k + 1 != nn)
{
p += z * zz[i, k + 2];
zz[i, k + 2] -= p * r;
}
zz[i, k + 1] -= p * q;
zz[i, k] -= p;
}
}
}
}
}
} while (l + 1 < nn);
}
if (anorm != 0.0)
{
for (nn = n - 1; nn >= 0; nn--)
{
p = (wri[nn].re);
q = (wri[nn].im);
int na = nn - 1;
//if (q == 0.0)
if (Math.Abs(q) <= float.Epsilon)
{
int m = nn;
a[nn, nn] = 1.0;
for (int i = nn - 1; i >= 0; i--)
{
double w = a[i, i] - p;
r = 0.0;
for (int j = m; j <= nn; j++)
{
r += a[i, j] * a[j, nn];
}
double s = 0.0;
if ((wri[i].im) < 0.0)
{
z = w;
s = r;
}
else
{
m = i;
//if ((wri[i].im) == 0.0)
if (Math.Abs(wri[i].im) <= float.Epsilon)
{
t = w;
//if (t == 0.0)
if (Math.Abs(t) <= float.Epsilon)
{
t = EPS * anorm;
}
a[i, nn] = -r / t;
}
else
{
double x = a[i, i + 1];
double y = a[i + 1, i];
q = Globals.SQR((wri[i].re) - p) + Globals.SQR((wri[i].im));
t = (x * s - z * r) / q;
a[i, nn] = t;
if (Math.Abs(x) > Math.Abs(z))
{
a[i + 1, nn] = (-r - w * t) / x;
}
else
{
a[i + 1, nn] = (-s - y * t) / z;
}
}
t = Math.Abs(a[i, nn]);
if (EPS * t * t > 1)
{
for (int j = i; j <= nn; j++)
{
a[j, nn] /= t;
}
}
}
}
}
else if (q < 0.0)
{
int m = na;
if (Math.Abs(a[nn, na]) > Math.Abs(a[na, nn]))
{
a[na, na] = q / a[nn, na];
a[na, nn] = -(a[nn, nn] - p) / a[nn, na];
}
else
{
Complex temp = new Complex(0.0, -a[na, nn]) / new Complex(a[na, na] - p, q);
a[na, na] = (temp.re);
a[na, nn] = (temp.im);
}
a[nn, na] = 0.0;
a[nn, nn] = 1.0;
for (int i = nn - 2; i >= 0; i--)
{
double w = a[i, i] - p;
double ra = 0.0;
double sa = 0.0;
double s = 0.0;
for (int j = m; j <= nn; j++)
{
ra += a[i, j] * a[j, na];
sa += a[i, j] * a[j, nn];
}
if ((wri[i].im) < 0.0)
{
z = w;
r = ra;
s = sa;
}
else
{
m = i;
//if ((wri[i].im) == 0.0)
if (Math.Abs(wri[i].im) <= float.Epsilon)
{
Complex temp = new Complex(-ra, -sa) / new Complex(w, q);
a[i, na] = (temp.re);
a[i, nn] = (temp.im);
}
else
{
double x = a[i, i + 1];
double y = a[i + 1, i];
double vr = Globals.SQR((wri[i].re) - p) + Globals.SQR((wri[i].im)) - q * q;
double vi = 2.0 * q * ((wri[i].re) - p);
//if (vr == 0.0 && vi == 0.0)
if (Math.Abs(vr) <= float.Epsilon && Math.Abs(vi) <= float.Epsilon)
{
vr = EPS * anorm * (Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z));
}
Complex temp = new Complex(x * r - z * ra + q * sa, x * s - z * sa - q * ra) / new Complex(vr, vi);
a[i, na] = (temp.re);
a[i, nn] = (temp.im);
if (Math.Abs(x) > Math.Abs(z) + Math.Abs(q))
{
a[i + 1, na] = (-ra - w * a[i, na] + q * a[i, nn]) / x;
a[i + 1, nn] = (-sa - w * a[i, nn] - q * a[i, na]) / x;
}
else
{
Complex tempx = new Complex(-r - y * a[i, na], -s - y * a[i, nn]) / new Complex(z, q);
a[i + 1, na] = (tempx.re);
a[i + 1, nn] = (tempx.im);
}
}
}
t = Math.Max(Math.Abs(a[i, na]), Math.Abs(a[i, nn]));
if (EPS * t * t > 1)
{
for (int j = i; j <= nn; j++)
{
a[j, na] /= t;
a[j, nn] /= t;
}
}
}
}
}
for (int j = n - 1; j >= 0; j--)
{
for (int i = 0; i < n; i++)
{
z = 0.0;
for (int k = 0; k <= j; k++)
{
z += zz[i, k] * a[k, j];
}
zz[i, j] = z;
}
}
}
}
public void sort()
{
for (int j = 1; j < n; j++)
{
Complex x = wri[j];
int i = j - 1;
for (; i >= 0; i--)
{
if ((wri[i].re) >= (x.re))
{
break;
}
wri[i + 1] = wri[i];
}
wri[i + 1] = x;
}
}
public void sortvecs()
{
double[] temp = new double[n];
for (int j = 1; j < n; j++)
{
Complex x = wri[j];
for (int k = 0; k < n; k++)
{
temp[k] = zz[k, j];
}
int i = j - 1;
for (; i >= 0; i--)
{
if ((wri[i].re) >= (x.re))
{
break;
}
wri[i + 1] = wri[i];
for (int k = 0; k < n; k++)
{
zz[k, i + 1] = zz[k, i];
}
}
wri[i + 1] = x;
for (int k = 0; k < n; k++)
{
zz[k, i + 1] = temp[k];
}
}
}
}
}
2 代码格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// Computes all eigenvalues and eigenvectors of a real nonsymmetric matrix by
/// reduction to Hes- senberg form followed by QR iteration. - PTC
/// </summary>
public class Unsymmeig
{
private int n { get; set; }
private double[,] a { get; set; }
private double[,] zz { get; set; }
public Complex[] wri { get; set; }
private double[] scale { get; set; }
private int[] perm { get; set; }
private bool yesvecs { get; set; }
private bool hessen { get; set; }
/// <summary>
/// Computes all eigenvalues and(optionally) eigenvectors of a real
/// nonsymmetric matrix a[0..n - 1][0..n - 1] by reduction to Hessenberg form
/// followed by QR iteration.If yesvecs is input as true (the default), then
/// the eigenvectors are computed.Otherwise, only the eigenvalues are
/// computed.If hessen is input as false (the default), the matrix is first
/// reduced to Hessenberg form.Otherwise it is assumed that the matrix is
/// already in Hessen- berg from.On output, wri[0..n - 1] contains the
/// eigenvalues of a sorted into descending order, while zz[0..n - 1][0..n - 1] is
/// a matrix whose columns contain the corresponding eigenvectors. For a
/// complex eigenvalue, only the eigenvector corresponding to the eigen- value
/// with a positive imaginary part is stored, with the real part in
/// zz[0..n-1][i] and the imaginary part in h.zz[0..n - 1][i + 1]. The eigenvectors
/// are not normalized.
/// </summary>
/// <param name="aa"></param>
/// <param name="yesvec"></param>
/// <param name="hessenb"></param>
public Unsymmeig(double[,] aa, bool yesvec = true, bool hessenb = false)
{
this.n = aa.GetLength(0);
this.a = aa;
this.zz = new double[n, n];
this.wri = new Complex[n];
this.scale = new double[n];
this.perm = new int[n];
this.yesvecs = yesvec;
this.hessen = hessenb;
for (int i = 0; i < n; i++)
{
scale[i] = 1.0;
}
balance();
if (!hessen)
{
elmhes();
}
if (yesvecs)
{
for (int i = 0; i < n; i++)
{
zz[i, i] = 1.0;
}
if (!hessen)
{
eltran();
}
hqr2();
balbak();
sortvecs();
}
else
{
hqr();
sort();
}
}
/// <summary>
/// Given a matrix a[0..n - 1][0..n-1], this routine replaces it by a balanced
/// matrix with identical eigenvalues.A symmetric matrix is already balanced
/// and is unaffected by this procedure.
/// </summary>
public void balance()
{
// https://en.cppreference.com/w/cpp/types/climits
double RADIX = 2.0;// numeric_limits<double>.radix;
bool done = false;
double sqrdx = RADIX * RADIX;
while (!done)
{
done = true;
for (int i = 0; i < n; i++)
{
double r = 0.0;
double c = 0.0;
for (int j = 0; j < n; j++)
{
if (j != i)
{
c += Math.Abs(a[j, i]);
r += Math.Abs(a[i, j]);
}
}
if (c != 0.0 && r != 0.0)
{
double g = r / RADIX;
double f = 1.0;
double s = c + r;
while (c < g)
{
f *= RADIX;
c *= sqrdx;
}
g = r * RADIX;
while (c > g)
{
f /= RADIX;
c /= sqrdx;
}
if ((c + r) / f < 0.95 * s)
{
done = false;
g = 1.0 / f;
scale[i] *= f;
for (int j = 0; j < n; j++)
{
a[i, j] *= g;
}
for (int j = 0; j < n; j++)
{
a[j, i] *= f;
}
}
}
}
}
}
/// <summary>
/// Forms the eigenvectors of a real nonsymmetric matrix by back transforming
/// those of the corre- sponding balanced matrix determined by balance.
/// </summary>
public void balbak()
{
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
zz[i, j] *= scale[i];
}
}
}
/// <summary>
/// Reduction to Hessenberg form by the elimination method.Replaces the real
/// nonsymmetric matrix a[0..n - 1][0..n - 1] by an upper Hessenberg matrix with
/// identical eigenvalues.Rec- ommended, but not required, is that this
/// routine be preceded by balance. On output, the Hessenberg matrix is in
/// elements a[i][j] with i <= j+1. Elements with i > j+1 are to be thought of
/// as zero, but are returned with random values.
/// </summary>
public void elmhes()
{
for (int m = 1; m < n - 1; m++)
{
double x = 0.0;
int i = m;
for (int j = m; j < n; j++)
{
if (Math.Abs(a[j, m - 1]) > Math.Abs(x))
{
x = a[j, m - 1];
i = j;
}
}
perm[m] = i;
if (i != m)
{
for (int j = m - 1; j < n; j++)
{
Globals.SWAP(ref a[i, j], ref a[m, j]);
}
for (int j = 0; j < n; j++)
{
Globals.SWAP(ref a[j, i], ref a[j, m]);
}
}
if (x != 0.0)
{
for (i = m + 1; i < n; i++)
{
double y = a[i, m - 1];
if (y != 0.0)
{
y /= x;
a[i, m - 1] = y;
for (int j = m; j < n; j++)
{
a[i, j] -= y * a[m, j];
}
for (int j = 0; j < n; j++)
{
a[j, m] += y * a[j, i];
}
}
}
}
}
}
/// <summary>
/// This routine accumulates the stabilized elementary similarity
/// transformations used in the reduction to upper Hessenberg form by elmhes.
/// The multipliers that were used in the reduction are obtained from the lower
/// triangle (below the subdiagonal) of a.The transformations are permuted
/// according to the permutations stored in perm by elmhes.
/// </summary>
public void eltran()
{
for (int mp = n - 2; mp > 0; mp--)
{
for (int k = mp + 1; k < n; k++)
{
zz[k, mp] = a[k, mp - 1];
}
int i = perm[mp];
if (i != mp)
{
for (int j = mp; j < n; j++)
{
zz[mp, j] = zz[i, j];
zz[i, j] = 0.0;
}
zz[i, mp] = 1.0;
}
}
}
public void hqr()
{
double r = 0.0;
double q = 0.0;
double p = 0.0;
double anorm = 0.0;
const double EPS = float.Epsilon;
for (int i = 0; i < n; i++)
{
for (int j = Math.Max(i - 1, 0); j < n; j++)
{
anorm += Math.Abs(a[i, j]);
}
}
int nn = n - 1;
double t = 0.0;
while (nn >= 0)
{
int its = 0;
int l;
do
{
for (l = nn; l > 0; l--)
{
double s = Math.Abs(a[l - 1, l - 1]) + Math.Abs(a[l, l]);
//if (s == 0.0)
if (s <= float.Epsilon)
{
s = anorm;
}
if (Math.Abs(a[l, l - 1]) <= EPS * s)
{
a[l, l - 1] = 0.0;
break;
}
}
double x = a[nn, nn];
if (l == nn)
{
wri[nn--] = new Complex(x + t);
}
else
{
double y = a[nn - 1, nn - 1];
double w = a[nn, nn - 1] * a[nn - 1, nn];
if (l == nn - 1)
{
p = 0.5 * (y - x);
q = p * p + w;
double z = Math.Sqrt(Math.Abs(q));
x += t;
if (q >= 0.0)
{
z = p + Globals.SIGN(z, p);
wri[nn - 1] = wri[nn] = new Complex(x + z);
if (z != 0.0)
{
wri[nn] = new Complex(x - w / z);
}
}
else
{
wri[nn] = new Complex(x + p, -z);
wri[nn - 1] = (wri[nn].conj());
}
nn -= 2;
}
else
{
if (its == 30)
{
throw new Exception("Too many iterations in hqr");
}
if (its == 10 || its == 20)
{
t += x;
for (int i = 0; i < nn + 1; i++)
{
a[i, i] -= x;
}
double s = Math.Abs(a[nn, nn - 1]) + Math.Abs(a[nn - 1, nn - 2]);
y = x = 0.75 * s;
w = -0.4375 * s * s;
}
++its;
int m = nn - 2;
for (; m >= l; m--)
{
double z = a[m, m];
r = x - z;
double s = y - z;
p = (r * s - w) / a[m + 1, m] + a[m, m + 1];
q = a[m + 1, m + 1] - z - r - s;
r = a[m + 2, m + 1];
s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);
p /= s;
q /= s;
r /= s;
if (m == l)
{
break;
}
double u = Math.Abs(a[m, m - 1]) * (Math.Abs(q) + Math.Abs(r));
double v = Math.Abs(p) * (Math.Abs(a[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(a[m + 1, m + 1]));
if (u <= EPS * v)
{
break;
}
}
for (int i = m; i < nn - 1; i++)
{
a[i + 2, i] = 0.0;
if (i != m) a[i + 2, i - 1] = 0.0;
}
for (int k = m; k < nn; k++)
{
if (k != m)
{
p = a[k, k - 1];
q = a[k + 1, k - 1];
r = 0.0;
if (k + 1 != nn) r = a[k + 2, k - 1];
if ((x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r)) != 0.0)
{
p /= x;
q /= x;
r /= x;
}
}
double s = Globals.SIGN(Math.Sqrt(p * p + q * q + r * r), p);
if ((s) != 0.0)
{
if (k == m)
{
if (l != m)
{
a[k, k - 1] = -a[k, k - 1];
}
}
else
{
a[k, k - 1] = -s * x;
}
p += s;
x = p / s;
y = q / s;
double z = r / s;
q /= p;
r /= p;
for (int j = k; j < nn + 1; j++)
{
p = a[k, j] + q * a[k + 1, j];
if (k + 1 != nn)
{
p += r * a[k + 2, j];
a[k + 2, j] -= p * z;
}
a[k + 1, j] -= p * y;
a[k, j] -= p * x;
}
int mmin = nn < k + 3 ? nn : k + 3;
for (int i = l; i < mmin + 1; i++)
{
p = x * a[i, k] + y * a[i, k + 1];
if (k + 1 != nn)
{
p += z * a[i, k + 2];
a[i, k + 2] -= p * r;
}
a[i, k + 1] -= p * q;
a[i, k] -= p;
}
}
}
}
}
} while (l + 1 < nn);
}
}
public void hqr2()
{
double z = 0.0;
double r = 0.0;
double q = 0.0;
double p = 0.0;
double anorm = 0.0;
double EPS = float.Epsilon;
for (int i = 0; i < n; i++)
{
for (int j = Math.Max(i - 1, 0); j < n; j++)
{
anorm += Math.Abs(a[i, j]);
}
}
int nn = n - 1;
double t = 0.0;
while (nn >= 0)
{
int its = 0;
int l;
do
{
for (l = nn; l > 0; l--)
{
double s = Math.Abs(a[l - 1, l - 1]) + Math.Abs(a[l, l]);
//if (s == 0.0)
if (s <= float.Epsilon)
{
s = anorm;
}
if (Math.Abs(a[l, l - 1]) <= EPS * s)
{
a[l, l - 1] = 0.0;
break;
}
}
double x = a[nn, nn];
if (l == nn)
{
a[nn, nn] = (x + t);
wri[nn] = new Complex(a[nn, nn]);
nn--;
}
else
{
double y = a[nn - 1, nn - 1];
double w = a[nn, nn - 1] * a[nn - 1, nn];
if (l == nn - 1)
{
p = 0.5 * (y - x);
q = p * p + w;
z = Math.Sqrt(Math.Abs(q));
x += t;
a[nn, nn] = x;
a[nn - 1, nn - 1] = y + t;
if (q >= 0.0)
{
z = p + Globals.SIGN(z, p);
wri[nn - 1] = wri[nn] = new Complex(x + z);
if (z != 0.0) wri[nn] = new Complex(x - w / z);
x = a[nn, nn - 1];
double s = Math.Abs(x) + Math.Abs(z);
p = x / s;
q = z / s;
r = Math.Sqrt(p * p + q * q);
p /= r;
q /= r;
for (int j = nn - 1; j < n; j++)
{
z = a[nn - 1, j];
a[nn - 1, j] = q * z + p * a[nn, j];
a[nn, j] = q * a[nn, j] - p * z;
}
for (int i = 0; i <= nn; i++)
{
z = a[i, nn - 1];
a[i, nn - 1] = q * z + p * a[i, nn];
a[i, nn] = q * a[i, nn] - p * z;
}
for (int i = 0; i < n; i++)
{
z = zz[i, nn - 1];
zz[i, nn - 1] = q * z + p * zz[i, nn];
zz[i, nn] = q * zz[i, nn] - p * z;
}
}
else
{
wri[nn] = new Complex(x + p, -z);
wri[nn - 1] = wri[nn].conj();
}
nn -= 2;
}
else
{
if (its == 30) throw new Exception("Too many iterations in hqr");
if (its == 10 || its == 20)
{
t += x;
for (int i = 0; i < nn + 1; i++)
{
a[i, i] -= x;
}
double s = Math.Abs(a[nn, nn - 1]) + Math.Abs(a[nn - 1, nn - 2]);
y = x = 0.75 * s;
w = -0.4375 * s * s;
}
++its;
int m = nn - 2;
for (; m >= l; m--)
{
z = a[m, m];
r = x - z;
double s = y - z;
p = (r * s - w) / a[m + 1, m] + a[m, m + 1];
q = a[m + 1, m + 1] - z - r - s;
r = a[m + 2, m + 1];
s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);
p /= s;
q /= s;
r /= s;
if (m == l)
{
break;
}
double u = Math.Abs(a[m, m - 1]) * (Math.Abs(q) + Math.Abs(r));
double v = Math.Abs(p) * (Math.Abs(a[m - 1, m - 1]) + Math.Abs(z) + Math.Abs(a[m + 1, m + 1]));
if (u <= EPS * v)
{
break;
}
}
for (int i = m; i < nn - 1; i++)
{
a[i + 2, i] = 0.0;
if (i != m)
{
a[i + 2, i - 1] = 0.0;
}
}
for (int k = m; k < nn; k++)
{
if (k != m)
{
p = a[k, k - 1];
q = a[k + 1, k - 1];
r = 0.0;
if (k + 1 != nn) r = a[k + 2, k - 1];
if ((x = Math.Abs(p) + Math.Abs(q) + Math.Abs(r)) != 0.0)
{
p /= x;
q /= x;
r /= x;
}
}
double s = Globals.SIGN(Math.Sqrt(p * p + q * q + r * r), p);
if ((s) != 0.0)
{
if (k == m)
{
if (l != m)
{
a[k, k - 1] = -a[k, k - 1];
}
}
else
{
a[k, k - 1] = -s * x;
}
p += s;
x = p / s;
y = q / s;
z = r / s;
q /= p;
r /= p;
for (int j = k; j < n; j++)
{
p = a[k, j] + q * a[k + 1, j];
if (k + 1 != nn)
{
p += r * a[k + 2, j];
a[k + 2, j] -= p * z;
}
a[k + 1, j] -= p * y;
a[k, j] -= p * x;
}
int mmin = nn < k + 3 ? nn : k + 3;
for (int i = 0; i < mmin + 1; i++)
{
p = x * a[i, k] + y * a[i, k + 1];
if (k + 1 != nn)
{
p += z * a[i, k + 2];
a[i, k + 2] -= p * r;
}
a[i, k + 1] -= p * q;
a[i, k] -= p;
}
for (int i = 0; i < n; i++)
{
p = x * zz[i, k] + y * zz[i, k + 1];
if (k + 1 != nn)
{
p += z * zz[i, k + 2];
zz[i, k + 2] -= p * r;
}
zz[i, k + 1] -= p * q;
zz[i, k] -= p;
}
}
}
}
}
} while (l + 1 < nn);
}
if (anorm != 0.0)
{
for (nn = n - 1; nn >= 0; nn--)
{
p = (wri[nn].re);
q = (wri[nn].im);
int na = nn - 1;
//if (q == 0.0)
if (Math.Abs(q) <= float.Epsilon)
{
int m = nn;
a[nn, nn] = 1.0;
for (int i = nn - 1; i >= 0; i--)
{
double w = a[i, i] - p;
r = 0.0;
for (int j = m; j <= nn; j++)
{
r += a[i, j] * a[j, nn];
}
double s = 0.0;
if ((wri[i].im) < 0.0)
{
z = w;
s = r;
}
else
{
m = i;
//if ((wri[i].im) == 0.0)
if (Math.Abs(wri[i].im) <= float.Epsilon)
{
t = w;
//if (t == 0.0)
if (Math.Abs(t) <= float.Epsilon)
{
t = EPS * anorm;
}
a[i, nn] = -r / t;
}
else
{
double x = a[i, i + 1];
double y = a[i + 1, i];
q = Globals.SQR((wri[i].re) - p) + Globals.SQR((wri[i].im));
t = (x * s - z * r) / q;
a[i, nn] = t;
if (Math.Abs(x) > Math.Abs(z))
{
a[i + 1, nn] = (-r - w * t) / x;
}
else
{
a[i + 1, nn] = (-s - y * t) / z;
}
}
t = Math.Abs(a[i, nn]);
if (EPS * t * t > 1)
{
for (int j = i; j <= nn; j++)
{
a[j, nn] /= t;
}
}
}
}
}
else if (q < 0.0)
{
int m = na;
if (Math.Abs(a[nn, na]) > Math.Abs(a[na, nn]))
{
a[na, na] = q / a[nn, na];
a[na, nn] = -(a[nn, nn] - p) / a[nn, na];
}
else
{
Complex temp = new Complex(0.0, -a[na, nn]) / new Complex(a[na, na] - p, q);
a[na, na] = (temp.re);
a[na, nn] = (temp.im);
}
a[nn, na] = 0.0;
a[nn, nn] = 1.0;
for (int i = nn - 2; i >= 0; i--)
{
double w = a[i, i] - p;
double ra = 0.0;
double sa = 0.0;
double s = 0.0;
for (int j = m; j <= nn; j++)
{
ra += a[i, j] * a[j, na];
sa += a[i, j] * a[j, nn];
}
if ((wri[i].im) < 0.0)
{
z = w;
r = ra;
s = sa;
}
else
{
m = i;
//if ((wri[i].im) == 0.0)
if (Math.Abs(wri[i].im) <= float.Epsilon)
{
Complex temp = new Complex(-ra, -sa) / new Complex(w, q);
a[i, na] = (temp.re);
a[i, nn] = (temp.im);
}
else
{
double x = a[i, i + 1];
double y = a[i + 1, i];
double vr = Globals.SQR((wri[i].re) - p) + Globals.SQR((wri[i].im)) - q * q;
double vi = 2.0 * q * ((wri[i].re) - p);
//if (vr == 0.0 && vi == 0.0)
if (Math.Abs(vr) <= float.Epsilon && Math.Abs(vi) <= float.Epsilon)
{
vr = EPS * anorm * (Math.Abs(w) + Math.Abs(q) + Math.Abs(x) + Math.Abs(y) + Math.Abs(z));
}
Complex temp = new Complex(x * r - z * ra + q * sa, x * s - z * sa - q * ra) / new Complex(vr, vi);
a[i, na] = (temp.re);
a[i, nn] = (temp.im);
if (Math.Abs(x) > Math.Abs(z) + Math.Abs(q))
{
a[i + 1, na] = (-ra - w * a[i, na] + q * a[i, nn]) / x;
a[i + 1, nn] = (-sa - w * a[i, nn] - q * a[i, na]) / x;
}
else
{
Complex tempx = new Complex(-r - y * a[i, na], -s - y * a[i, nn]) / new Complex(z, q);
a[i + 1, na] = (tempx.re);
a[i + 1, nn] = (tempx.im);
}
}
}
t = Math.Max(Math.Abs(a[i, na]), Math.Abs(a[i, nn]));
if (EPS * t * t > 1)
{
for (int j = i; j <= nn; j++)
{
a[j, na] /= t;
a[j, nn] /= t;
}
}
}
}
}
for (int j = n - 1; j >= 0; j--)
{
for (int i = 0; i < n; i++)
{
z = 0.0;
for (int k = 0; k <= j; k++)
{
z += zz[i, k] * a[k, j];
}
zz[i, j] = z;
}
}
}
}
public void sort()
{
for (int j = 1; j < n; j++)
{
Complex x = wri[j];
int i = j - 1;
for (; i >= 0; i--)
{
if ((wri[i].re) >= (x.re))
{
break;
}
wri[i + 1] = wri[i];
}
wri[i + 1] = x;
}
}
public void sortvecs()
{
double[] temp = new double[n];
for (int j = 1; j < n; j++)
{
Complex x = wri[j];
for (int k = 0; k < n; k++)
{
temp[k] = zz[k, j];
}
int i = j - 1;
for (; i >= 0; i--)
{
if ((wri[i].re) >= (x.re))
{
break;
}
wri[i + 1] = wri[i];
for (int k = 0; k < n; k++)
{
zz[k, i + 1] = zz[k, i];
}
}
wri[i + 1] = x;
for (int k = 0; k < n; k++)
{
zz[k, i + 1] = temp[k];
}
}
}
}
}