C#，码海拾贝（29）——求解“大型稀疏方程组”的“全选主元高斯

大型稀疏矩阵线性化方程组的数值求解问题广泛存在于工程实践尤其是计算机仿真领域如水力管网计算，电力系统的大型导纳矩阵计算，高阶偏微分方程的数值求解，以及铸件充型过程与凝固过程的数值模拟等。

经常出现在科学和工程计算中, 因此寻找稀疏线性方程组的高效计算方法及其并行算法是提高科学与工程应用问题计算效率的有效途径。

这方面的研究始于 20 世纪 60 年代然而至今仍未见任何公诸于世的高效的数值解法计算准确~ 快速又节省内存。

所谓稀疏矩阵，是指矩阵中大多数元素为零元素。而大型稀疏矩阵，顾名思义是指矩阵的阶数很大的稀疏矩阵比如阶数大于 10 万甚至 100 万或更大。鉴于稀疏矩阵的特点在存储时不必存储矩阵中的零元素只存储矩阵中的非零元素这种存储方式称之为稀疏矩阵的压缩存储. 通常稀疏矩阵的存储方式有两种: 三元组法与十字链表法。

大型稀疏矩阵线性化方程组的数值求解问题广泛存在于工程实践尤其是计算机仿真领域如水力管网计算，电力系统的大型导纳矩阵计算，高阶偏微分方程的数值求解，以及铸件充型过程与凝固过程的数值模拟等。
经常出现在科学和工程计算中, 因此寻找稀疏线性方程组的高效计算方法及其并行算法是提高科学与工程应用问题计算效率的有效途径。
这方面的研究始于 20 世纪 60 年代然而至今仍未见任何公诸于世的高效的数值解法计算准确~ 快速又节省内存。
所谓稀疏矩阵，是指矩阵中大多数元素为零元素。而大型稀疏矩阵，顾名思义是指矩阵的阶数很大的稀疏矩阵比如阶数大于 10 万甚至 100 万或更大。鉴于稀疏矩阵的特点在存储时不必存储矩阵中的零元素只存储矩阵中的非零元素这种存储方式称之为稀疏矩阵的压缩存储. 通常稀疏矩阵的存储方式有两种: 三元组法与十字链表法。

using System;

namespace Zhou.CSharp.Algorithm
{
/// <summary>
/// 求解线性方程组的类 LEquations
/// 原作周长发
/// 改编深度混淆
/// </summary>
public static partial class LEquations
{

/// <summary>
/// 求解大型稀疏方程组的全选主元高斯－约去消去法
/// </summary>
/// <param name="mtxLECoef">指定的系数矩阵</param>
/// <param name="mtxLEConst">指定的常数矩阵</param>
/// <param name="mtxResult">Matrix引用对象，返回方程组解矩阵</param>
/// <return>bool 型，方程组求解是否成功</return>
public static bool GetRootsetGgje(Matrix mtxLECoef, Matrix mtxLEConst, Matrix mtxResult)
{
int i, j, k, nIs = 0, u, v;
double d, t;

// 方程组属性，将常数矩阵赋给解矩阵
Matrix mtxCoef = new Matrix(mtxLECoef);
mtxResult.SetValue(mtxLEConst);
int n = mtxCoef.GetNumColumns();
double[] pDataCoef = mtxCoef.GetData();
double[] pDataConst = mtxResult.GetData();

// 临时缓冲区，存放变换的列数
int[] pnJs = new int[n];

// 消元
for (k = 0; k <= n - 1; k++)
{
d = 0.0;
for (i = k; i <= n - 1; i++)
{
for (j = k; j <= n - 1; j++)
{
t = Math.Abs(pDataCoef[i * n + j]);
if (t > d)
{
d = t;
pnJs[k] = j;
nIs = i;
}
}
}

if (Math.Abs(d) < float.Epsilon)
{
return false;
}

if (nIs != k)
{
for (j = k; j <= n - 1; j++)
{
u = k * n + j;
v = nIs * n + j;
t = pDataCoef[u];
pDataCoef[u] = pDataCoef[v];
pDataCoef[v] = t;
}

t = pDataConst[k];
pDataConst[k] = pDataConst[nIs];
pDataConst[nIs] = t;
}

if (pnJs[k] != k)
{
for (i = 0; i <= n - 1; i++)
{
u = i * n + k;
v = i * n + pnJs[k];
t = pDataCoef[u];
pDataCoef[u] = pDataCoef[v];
pDataCoef[v] = t;
}
}

t = pDataCoef[k * n + k];
for (j = k + 1; j <= n - 1; j++)
{
u = k * n + j;
if (Math.Abs(pDataCoef[u]) > float.Epsilon)
{
pDataCoef[u] = pDataCoef[u] / t;
}
}

pDataConst[k] = pDataConst[k] / t;
for (j = k + 1; j <= n - 1; j++)
{
u = k * n + j;
if (Math.Abs(pDataCoef[u]) > float.Epsilon)
{
for (i = 0; i <= n - 1; i++)
{
v = i * n + k;
if ((i != k) && Math.Abs(pDataCoef[v]) > float.Epsilon)
{
nIs = i * n + j;
pDataCoef[nIs] = pDataCoef[nIs] - pDataCoef[v] * pDataCoef[u];
}
}
}
}

for (i = 0; i <= n - 1; i++)
{
u = i * n + k;
if ((i != k) && Math.Abs(pDataCoef[u]) > float.Epsilon)
{
pDataConst[i] = pDataConst[i] - pDataCoef[u] * pDataConst[k];
}
}
}

// 调整
for (k = n - 1; k >= 0; k--)
{
if (k != pnJs[k])
{
t = pDataConst[k];
pDataConst[k] = pDataConst[pnJs[k]];
pDataConst[pnJs[k]] = t;
}
}

return true;
}

}
}