1 文本格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// Computes all eigenvalues and eigenvectors of
/// a real symmetric matrix by Jacobi's method.
/// </summary>
public class Jacobi
{
private int n { get; set; }
private double[,] a;
private double[,] v;
private double[] d;
private int nrot { get; set; }
private double EPS { get; set; }
/// <summary>
/// Computes all eigenvalues and eigenvectors of a real symmetric matrix
/// a[0..n - 1][0..n - 1]. On output, d[0..n - 1] contains the eigenvalues of a
/// sorted into descending order, while v[0..n - 1][0..n - 1] is a matrix whose
/// columns contain the corresponding normalized eigenvectors.nrot contains
/// the number of Jacobi rotations that were required.Only the upper triangle
/// of a is accessed.
/// </summary>
/// <param name="aa"></param>
/// <exception cref="Exception"></exception>
public Jacobi(double[,] aa)
{
this.n = aa.GetLength(0);
this.a = aa;
this.v = new double[n, n];
this.d = new double[n];
this.nrot = 0;
this.EPS = float.Epsilon;
double[] b = new double[n];
double[] z = new double[n];
for (int ip = 0; ip < n; ip++)
{
for (int iq = 0; iq < n; iq++)
{
v[ip, iq] = 0.0;
}
v[ip, ip] = 1.0;
}
for (int ip = 0; ip < n; ip++)
{
b[ip] = d[ip] = a[ip, ip];
z[ip] = 0.0;
}
for (int i = 1; i <= 50; i++)
{
double sm = 0.0;
for (int ip = 0; ip < n - 1; ip++)
{
for (int iq = ip + 1; iq < n; iq++)
{
sm += Math.Abs(a[ip, iq]);
}
}
//if (sm == 0.0)
if (Math.Abs(sm) <= float.Epsilon)
{
eigsrt( d, v);
return;
}
double tresh;
if (i < 4)
{
tresh = 0.2 * sm / (n * n);
}
else
{
tresh = 0.0;
}
for (int ip = 0; ip < n - 1; ip++)
{
for (int iq = ip + 1; iq < n; iq++)
{
double g = 100.0 * Math.Abs(a[ip, iq]);
if (i > 4 && g <= EPS * Math.Abs(d[ip]) && g <= EPS * Math.Abs(d[iq]))
{
a[ip, iq] = 0.0;
}
else if (Math.Abs(a[ip, iq]) > tresh)
{
double h = d[iq] - d[ip];
double t;
if (g <= EPS * Math.Abs(h))
{
t = (a[ip, iq]) / h;
}
else
{
double theta = 0.5 * h / (a[ip, iq]);
t = 1.0 / (Math.Abs(theta) + Math.Sqrt(1.0 + theta * theta));
if (theta < 0.0)
{
t = -t;
}
}
double c = 1.0 / Math.Sqrt(1 + t * t);
double s = t * c;
double tau = s / (1.0 + c);
h = t * a[ip, iq];
z[ip] -= h;
z[iq] += h;
d[ip] -= h;
d[iq] += h;
a[ip, iq] = 0.0;
for (int j = 0; j < ip; j++)
{
rot( a, s, tau, j, ip, j, iq);
}
for (int j = ip + 1; j < iq; j++)
{
rot( a, s, tau, ip, j, j, iq);
}
for (int j = iq + 1; j < n; j++)
{
rot( a, s, tau, ip, j, iq, j);
}
for (int j = 0; j < n; j++)
{
rot( v, s, tau, j, ip, j, iq);
}
++nrot;
}
}
}
for (int ip = 0; ip < n; ip++)
{
b[ip] += z[ip];
d[ip] = b[ip];
z[ip] = 0.0;
}
}
throw new Exception("Too many iterations in routine jacobi");
}
public void rot(double[,] a, double s, double tau, int i, int j, int k, int l)
{
double g = a[i, j];
double h = a[k, l];
a[i, j] = g - s * (h + g * tau);
a[k, l] = h + s * (g - h * tau);
}
/// <summary>
/// Given the eigenvalues d[0..n - 1] and(optionally) the eigenvectors
/// v[0..n - 1][0..n - 1] as determined by Jacobi or tqli, this routine sorts the
/// eigenvalues into descending order and rearranges the columns of v
/// correspondingly.The method is straight insertion.
/// </summary>
/// <param name="d"></param>
/// <param name="v"></param>
public static void eigsrt(double[] d, double[,] v)
{
int k;
int n = d.Length;
for (int i = 0; i < n - 1; i++)
{
double p = d[k = i];
for (int j = i; j < n; j++)
{
if (d[j] >= p)
{
p = d[k = j];
}
}
if (k != i)
{
d[k] = d[i];
d[i] = p;
if (v != null)
{
for (int j = 0; j < n; j++)
{
p = v[j, i];
v[j, i] = v[j, k];
v[j, k] = p;
}
}
}
}
}
}
}
2 代码格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// Computes all eigenvalues and eigenvectors of
/// a real symmetric matrix by Jacobi's method.
/// </summary>
public class Jacobi
{
private int n { get; set; }
private double[,] a;
private double[,] v;
private double[] d;
private int nrot { get; set; }
private double EPS { get; set; }
/// <summary>
/// Computes all eigenvalues and eigenvectors of a real symmetric matrix
/// a[0..n - 1][0..n - 1]. On output, d[0..n - 1] contains the eigenvalues of a
/// sorted into descending order, while v[0..n - 1][0..n - 1] is a matrix whose
/// columns contain the corresponding normalized eigenvectors.nrot contains
/// the number of Jacobi rotations that were required.Only the upper triangle
/// of a is accessed.
/// </summary>
/// <param name="aa"></param>
/// <exception cref="Exception"></exception>
public Jacobi(double[,] aa)
{
this.n = aa.GetLength(0);
this.a = aa;
this.v = new double[n, n];
this.d = new double[n];
this.nrot = 0;
this.EPS = float.Epsilon;
double[] b = new double[n];
double[] z = new double[n];
for (int ip = 0; ip < n; ip++)
{
for (int iq = 0; iq < n; iq++)
{
v[ip, iq] = 0.0;
}
v[ip, ip] = 1.0;
}
for (int ip = 0; ip < n; ip++)
{
b[ip] = d[ip] = a[ip, ip];
z[ip] = 0.0;
}
for (int i = 1; i <= 50; i++)
{
double sm = 0.0;
for (int ip = 0; ip < n - 1; ip++)
{
for (int iq = ip + 1; iq < n; iq++)
{
sm += Math.Abs(a[ip, iq]);
}
}
//if (sm == 0.0)
if (Math.Abs(sm) <= float.Epsilon)
{
eigsrt( d, v);
return;
}
double tresh;
if (i < 4)
{
tresh = 0.2 * sm / (n * n);
}
else
{
tresh = 0.0;
}
for (int ip = 0; ip < n - 1; ip++)
{
for (int iq = ip + 1; iq < n; iq++)
{
double g = 100.0 * Math.Abs(a[ip, iq]);
if (i > 4 && g <= EPS * Math.Abs(d[ip]) && g <= EPS * Math.Abs(d[iq]))
{
a[ip, iq] = 0.0;
}
else if (Math.Abs(a[ip, iq]) > tresh)
{
double h = d[iq] - d[ip];
double t;
if (g <= EPS * Math.Abs(h))
{
t = (a[ip, iq]) / h;
}
else
{
double theta = 0.5 * h / (a[ip, iq]);
t = 1.0 / (Math.Abs(theta) + Math.Sqrt(1.0 + theta * theta));
if (theta < 0.0)
{
t = -t;
}
}
double c = 1.0 / Math.Sqrt(1 + t * t);
double s = t * c;
double tau = s / (1.0 + c);
h = t * a[ip, iq];
z[ip] -= h;
z[iq] += h;
d[ip] -= h;
d[iq] += h;
a[ip, iq] = 0.0;
for (int j = 0; j < ip; j++)
{
rot( a, s, tau, j, ip, j, iq);
}
for (int j = ip + 1; j < iq; j++)
{
rot( a, s, tau, ip, j, j, iq);
}
for (int j = iq + 1; j < n; j++)
{
rot( a, s, tau, ip, j, iq, j);
}
for (int j = 0; j < n; j++)
{
rot( v, s, tau, j, ip, j, iq);
}
++nrot;
}
}
}
for (int ip = 0; ip < n; ip++)
{
b[ip] += z[ip];
d[ip] = b[ip];
z[ip] = 0.0;
}
}
throw new Exception("Too many iterations in routine jacobi");
}
public void rot(double[,] a, double s, double tau, int i, int j, int k, int l)
{
double g = a[i, j];
double h = a[k, l];
a[i, j] = g - s * (h + g * tau);
a[k, l] = h + s * (g - h * tau);
}
/// <summary>
/// Given the eigenvalues d[0..n - 1] and(optionally) the eigenvectors
/// v[0..n - 1][0..n - 1] as determined by Jacobi or tqli, this routine sorts the
/// eigenvalues into descending order and rearranges the columns of v
/// correspondingly.The method is straight insertion.
/// </summary>
/// <param name="d"></param>
/// <param name="v"></param>
public static void eigsrt(double[] d, double[,] v)
{
int k;
int n = d.Length;
for (int i = 0; i < n - 1; i++)
{
double p = d[k = i];
for (int j = i; j < n; j++)
{
if (d[j] >= p)
{
p = d[k = j];
}
}
if (k != i)
{
d[k] = d[i];
d[i] = p;
if (v != null)
{
for (int j = 0; j < n; j++)
{
p = v[j, i];
v[j, i] = v[j, k];
v[j, k] = p;
}
}
}
}
}
}
}
