using System;
namespace Zhou.CSharp.Algorithm
{
/// <summary>
/// 矩阵类
/// 作者：周长发
/// 改进：深度混淆
/// https://blog.csdn.net/beijinghorn
/// </summary>
public partial class Matrix
{
/// <summary>
/// 实对称三对角阵的全部特征值与特征向量的计算
/// </summary>
/// <param name="src">源矩阵</param>
/// <param name="dblB">一维数组，长度为矩阵的阶数，传入对称三对角阵的主对角线元素；返回时存放全部特征值。</param>
/// <param name="dblC">一维数组，长度为矩阵的阶数，前n-1个元素传入对称三对角阵的次对角线元素</param>
/// <param name="mtxQ">如果传入单位矩阵，则返回实对称三对角阵的特征值向量矩阵；如果传入MakeSymTri函数求得的矩阵A的豪斯荷尔德变换的乘积矩阵Q，则返回矩阵A的特征值向量矩阵。其中第i列为与数组dblB中第j个特征值对应的特征向量。</param>
/// <param name="nMaxIt">迭代次数</param>
/// <param name="eps">计算精度</param>
/// <returns>求解是否成功</returns>
public static bool ComputeEvSymTri(Matrix src, double[] dblB, double[] dblC, Matrix mtxQ, int nMaxIt, double eps)
{
int i, j, k, m, it, u, v;
double d, f, h, g, p, r, e, s;
// 初值
int n = mtxQ.GetNumColumns();
dblC[n - 1] = 0.0;
d = 0.0;
f = 0.0;
// 迭代计算
for (j = 0; j <= n - 1; j++)
{
it = 0;
h = eps * (Math.Abs(dblB[j]) + Math.Abs(dblC[j]));
if (h > d)
{
d = h;
}
m = j;
while ((m <= n - 1) && (Math.Abs(dblC[m]) > d))
{
m = m + 1;
}
if (m != j)
{
do
{
if (it == nMaxIt)
{
return false;
}
it = it + 1;
g = dblB[j];
p = (dblB[j + 1] - g) / (2.0 * dblC[j]);
r = Math.Sqrt(p * p + 1.0);
if (p >= 0.0)
{
dblB[j] = dblC[j] / (p + r);
}
else
{
dblB[j] = dblC[j] / (p - r);
}
h = g - dblB[j];
for (i = j + 1; i <= n - 1; i++)
{
dblB[i] = dblB[i] - h;
}
f = f + h;
p = dblB[m];
e = 1.0;
s = 0.0;
for (i = m - 1; i >= j; i--)
{
g = e * dblC[i];
h = e * p;
if (Math.Abs(p) >= Math.Abs(dblC[i]))
{
e = dblC[i] / p;
r = Math.Sqrt(e * e + 1.0);
dblC[i + 1] = s * p * r;
s = e / r;
e = 1.0 / r;
}
else
{
e = p / dblC[i];
r = Math.Sqrt(e * e + 1.0);
dblC[i + 1] = s * dblC[i] * r;
s = 1.0 / r;
e = e / r;
}
p = e * dblB[i] - s * g;
dblB[i + 1] = h + s * (e * g + s * dblB[i]);
for (k = 0; k <= n - 1; k++)
{
u = k * n + i + 1;
v = u - 1;
h = mtxQ[u];
mtxQ[u] = s * mtxQ[v] + e * h;
mtxQ[v] = e * mtxQ[v] - s * h;
}
}
dblC[j] = s * p;
dblB[j] = e * p;
} while (Math.Abs(dblC[j]) > d);
}
dblB[j] = dblB[j] + f;
}
for (i = 0; i <= n - 1; i++)
{
k = i;
p = dblB[i];
if (i + 1 <= n - 1)
{
j = i + 1;
while ((j <= n - 1) && (dblB[j] <= p))
{
k = j;
p = dblB[j];
j = j + 1;
}
}
if (k != i)
{
dblB[k] = dblB[i];
dblB[i] = p;
for (j = 0; j <= n - 1; j++)
{
u = j * n + i;
v = j * n + k;
p = mtxQ[u];
mtxQ[u] = mtxQ[v];
mtxQ[v] = p;
}
}
}
return true;
}
}
}
