当前位置：首页 > news >正文

C#，数值计算——计算实对称矩阵所有特征值和特征向量的雅可比（Jacobi）方法与源程序

news 来源：原创 2024/5/20 11:38:49

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;}}}}}}
}