在大型稀疏方程组中，最常见的是
带状方程组
，其系数矩阵是带状矩阵，非零元素仅集中在对角线附近的带状区域内。
特别的，当上下带宽相等时我们A称方程组为等带宽方程组。总带宽为3的等带宽方程组我们又叫三对角方程组。
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=”nBandWidth”>带宽</param>
/// <param name=”mtxResult”>Matrix对象，返回方程组解矩阵</param>
/// <return>bool 型，方程组求解是否成功</return>
public static bool GetRootsetBand(Matrix mtxLECoef, Matrix mtxLEConst, int nBandWidth, Matrix mtxResult)
{
int r, k, i, j, nis = 0, u, v;
double p, t;
// 带宽必须为奇数
if ((nBandWidth – 1) % 2 != 0)
{
return false;
}
// 将常数矩阵赋给解矩阵
mtxResult.SetValue(mtxLEConst);
double[] pDataConst = mtxResult.GetData();
// 方程组属性
int n = mtxLECoef.GetNumColumns();
int m = mtxLEConst.GetNumColumns();
if (mtxLECoef.GetNumRows() != n)
{
return false;
}
// 带宽数组：带型矩阵的有效数据
double[] pBandData = new double[nBandWidth * n];
// 半带宽
r = (nBandWidth – 1) / 2;
// 构造带宽数组
for (i = 0; i < n; ++i)
{
j = 0;
for (k = Math.Max(0, i – r); k < Math.Max(0, i – r) + nBandWidth; ++k)
{
if (k < n)
{
pBandData[i * nBandWidth + j++] = mtxLECoef.GetElement(i, k);
}
else
{
pBandData[i * nBandWidth + j++] = 0;
}
}
}
// 求解
for (k = 0; k <= n – 2; k++)
{
p = 0.0;
for (i = k; i <= r; i++)
{
t = Math.Abs(pBandData[i * nBandWidth]);
if (t > p)
{
p = t;
nis = i;
}
}
if (Math.Abs(p) < float.Epsilon)// p == 0.0)
{
return false;
}
for (j = 0; j <= m – 1; j++)
{
u = k * m + j;
v = nis * m + j;
t = pDataConst[u];
pDataConst[u] = pDataConst[v];
pDataConst[v] = t;
}
for (j = 0; j <= nBandWidth – 1; j++)
{
u = k * nBandWidth + j;
v = nis * nBandWidth + j;
t = pBandData[u];
pBandData[u] = pBandData[v];
pBandData[v] = t;
}
for (j = 0; j <= m – 1; j++)
{
u = k * m + j;
pDataConst[u] = pDataConst[u] / pBandData[k * nBandWidth];
}
for (j = 1; j <= nBandWidth – 1; j++)
{
u = k * nBandWidth + j;
pBandData[u] = pBandData[u] / pBandData[k * nBandWidth];
}
for (i = k + 1; i <= r; i++)
{
t = pBandData[i * nBandWidth];
for (j = 0; j <= m – 1; j++)
{
u = i * m + j;
v = k * m + j;
pDataConst[u] = pDataConst[u] – t * pDataConst[v];
}
for (j = 1; j <= nBandWidth – 1; j++)
{
u = i * nBandWidth + j;
v = k * nBandWidth + j;
pBandData[u – 1] = pBandData[u] – t * pBandData[v];
}
u = i * nBandWidth + nBandWidth – 1; pBandData[u] = 0.0;
}
if (r != (n – 1))
{
r = r + 1;
}
}
p = pBandData[(n – 1) * nBandWidth];
if (Math.Abs(p) < float.Epsilon)// p == 0.0)
{
return false;
}
for (j = 0; j <= m – 1; j++)
{
u = (n – 1) * m + j;
pDataConst[u] = pDataConst[u] / p;
}
r = 1;
for (i = n – 2; i >= 0; i–)
{
for (k = 0; k <= m – 1; k++)
{
u = i * m + k;
for (j = 1; j <= r; j++)
{
v = i * nBandWidth + j;
nis = (i + j) * m + k;
pDataConst[u] = pDataConst[u] – pBandData[v] * pDataConst[nis];
}
}
if (r != (nBandWidth – 1))
{
r = r + 1;
}
}
return true;
}
}
}