C#，码海拾贝（34）——求“赫申伯格矩阵”全部“特征值”的“QR方法”之C#源代码

news2024/10/6 9:20:13

using System;

namespace Zhou.CSharp.Algorithm
{
/// <summary>
/// 矩阵类
/// 作者：周长发
/// 改进：深度混淆
/// https://blog.csdn.net/beijinghorn
/// </summary>
public partial class Matrix
{

/// <summary>
/// 求赫申伯格矩阵全部特征值的QR方法
/// </summary>
/// <param name="src">源矩阵</param>
/// <param name="dblU">一维数组，长度为矩阵的阶数，返回时存放特征值的实部</param>
/// <param name="dblV">一维数组，长度为矩阵的阶数，返回时存放特征值的虚部</param>
/// <param name="nMaxIt">迭代次数</param>
/// <param name="eps">计算精度</param>
/// <returns>求解是否成功</returns>
public static bool ComputeEvHBerg(Matrix src, out double[] dblU, out double[] dblV, int nMaxIt = 100, double eps = 1.0E-7)
{
int m, it, i, j, k, ii, jj, kk;
double b, c, w, g, xy, p, q, r, x, s, e, f, z, y;

int n = src.Columns;
dblU = new double[n];
dblV = new double[n];

it = 0;
m = n;
while (m != 0)
{
int u = m - 1;
while ((u > 0) && (Math.Abs(src[u * n + u - 1]) >
eps * (Math.Abs(src[(u - 1) * n + u - 1]) + Math.Abs(src[u * n + u]))))
{
u = u - 1;
}
ii = (m - 1) * n + m - 1;
jj = (m - 1) * n + m - 2;
kk = (m - 2) * n + m - 1;
int uu = (m - 2) * n + m - 2;
if (u == m - 1)
{
dblU[m - 1] = src[(m - 1) * n + m - 1];
dblV[m - 1] = 0.0;
m = m - 1;
it = 0;
}
else if (u == m - 2)
{
b = -(src[ii] + src[uu]);
c = src[ii] * src[uu] - src[jj] * src[kk];
w = b * b - 4.0 * c;
if (Math.Abs(w) < float.Epsilon)
{
return false;
}
y = Math.Sqrt(Math.Abs(w));
if (w > 0.0)
{
xy = 1.0;
if (b < 0.0)
{
xy = -1.0;
}
dblU[m - 1] = (-b - xy * y) / 2.0;
dblU[m - 2] = c / dblU[m - 1];
dblV[m - 1] = 0.0; dblV[m - 2] = 0.0;
}
else
{
dblU[m - 1] = -b / 2.0;
dblU[m - 2] = dblU[m - 1];
dblV[m - 1] = y / 2.0;
dblV[m - 2] = -dblV[m - 1];
}

m = m - 2;
it = 0;
}
else
{
if (it >= nMaxIt)
{
return false;
}
it = it + 1;
for (j = u + 2; j <= m - 1; j++)
{
src[j * n + j - 2] = 0.0;
}
for (j = u + 3; j <= m - 1; j++)
{
src[j * n + j - 3] = 0.0;
}
for (k = u; k <= m - 2; k++)
{
if (k != u)
{
p = src[k * n + k - 1];
q = src[(k + 1) * n + k - 1];
r = 0.0;
if (k != m - 2)
{
r = src[(k + 2) * n + k - 1];
}
}
else
{
x = src[ii] + src[uu];
y = src[uu] * src[ii] - src[kk] * src[jj];
ii = u * n + u;
jj = u * n + u + 1;
kk = (u + 1) * n + u;
uu = (u + 1) * n + u + 1;
p = src[ii] * (src[ii] - x) + src[jj] * src[kk] + y;
q = src[kk] * (src[ii] + src[uu] - x);
r = src[kk] * src[(u + 2) * n + u + 1];
}

if ((Math.Abs(p) + Math.Abs(q) + Math.Abs(r)) > float.Epsilon)// != 0.0)
{
xy = 1.0;
if (p < 0.0)
{
xy = -1.0;
}
s = xy * Math.Sqrt(p * p + q * q + r * r);
if (k != u)
{
src[k * n + k - 1] = -s;
}
e = -q / s;
f = -r / s;
x = -p / s;
y = -x - f * r / (p + s);
g = e * r / (p + s);
z = -x - e * q / (p + s);
for (j = k; j <= m - 1; j++)
{
ii = k * n + j;
jj = (k + 1) * n + j;
p = x * src[ii] + e * src[jj];
q = e * src[ii] + y * src[jj];
r = f * src[ii] + g * src[jj];
if (k != m - 2)
{
kk = (k + 2) * n + j;
p = p + f * src[kk];
q = q + g * src[kk];
r = r + z * src[kk];
src[kk] = r;
}

src[jj] = q; src[ii] = p;
}

j = k + 3;
if (j >= m - 1)
{
j = m - 1;
}
for (i = u; i <= j; i++)
{
ii = i * n + k;
jj = i * n + k + 1;
p = x * src[ii] + e * src[jj];
q = e * src[ii] + y * src[jj];
r = f * src[ii] + g * src[jj];
if (k != m - 2)
{
kk = i * n + k + 2;
p = p + f * src[kk];
q = q + g * src[kk];
r = r + z * src[kk];
src[kk] = r;
}