C#，码海拾贝（24）——求解“复系数方程组”的“全选主元高斯－约当消去法”之C#源代码，《C#数值计算算法编程》源代码升级改进版

news2024/11/22 11:59:13

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="mtxCoefImag">系数矩阵的虚部矩阵</param>
/// <param name="mtxConstImag">常数矩阵的虚部矩阵</param>
/// <param name="mtxResult">Matrix对象，返回方程组解矩阵的实部矩阵</param>
/// <param name="mtxResultImag">Matrix对象，返回方程组解矩阵的虚部矩阵</param>
/// <return>bool 型，方程组求解是否成功</return>
public static bool GetRootsetGaussJordan(Matrix mtxLECoef, Matrix mtxLEConst, Matrix mtxCoefImag, Matrix mtxConstImag, Matrix mtxResult, Matrix mtxResultImag)
{
int r, k, i, j, nIs = 0, u, v;
double p, q, s, d;

// 方程组的属性，将常数矩阵赋给解矩阵
mtxResult.SetValue(mtxLEConst);
mtxResultImag.SetValue(mtxConstImag);
double[] pDataCoef = mtxLECoef.GetData();
double[] pDataConst = mtxResult.GetData();
double[] pDataCoefImag = mtxCoefImag.GetData();
double[] pDataConstImag = mtxResultImag.GetData();
int n = mtxLECoef.GetNumColumns();
int m = mtxLEConst.GetNumColumns();

// 临时缓冲区，存放变换的列数
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++)
{
u = i * n + j;
p = pDataCoef[u] * pDataCoef[u] + pDataCoefImag[u] * pDataCoefImag[u];
if (p > d)
{
d = p;
pnJs[k] = j;
nIs = i;
}
}
}

// 求解失败
if (Math.Abs(d) < float.Epsilon)// d == 0.0)
{
return false;
}

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

for (j = 0; j <= m - 1; j++)
{
u = k * m + j;
v = nIs * m + j;
p = pDataConst[u];
pDataConst[u] = pDataConst[v];
pDataConst[v] = p;
p = pDataConstImag[u];
pDataConstImag[u] = pDataConstImag[v];
pDataConstImag[v] = p;
}
}

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

v = k * n + k;
for (j = k + 1; j <= n - 1; j++)
{
u = k * n + j;
p = pDataCoef[u] * pDataCoef[v];
q = -pDataCoefImag[u] * pDataCoefImag[v];
s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataCoef[u] + pDataCoefImag[u]);
pDataCoef[u] = (p - q) / d;
pDataCoefImag[u] = (s - p - q) / d;
}

for (j = 0; j <= m - 1; j++)
{
u = k * m + j;
p = pDataConst[u] * pDataCoef[v];
q = -pDataConstImag[u] * pDataCoefImag[v];
s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataConst[u] + pDataConstImag[u]);
pDataConst[u] = (p - q) / d;
pDataConstImag[u] = (s - p - q) / d;
}

for (i = 0; i <= n - 1; i++)
{
if (i != k)
{
u = i * n + k;
for (j = k + 1; j <= n - 1; j++)
{
v = k * n + j;
u = i * n + j;
p = pDataCoef[u] * pDataCoef[v];
q = pDataCoefImag[u] * pDataCoefImag[v];
s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataCoef[v] + pDataCoefImag[v]);
pDataCoef[u] = pDataCoef[u] - p + q;
pDataCoefImag[u] = pDataCoefImag[u] - s + p + q;
}

for (j = 0; j <= m - 1; j++)
{
u = i * m + j;
v = k * m + j;
p = pDataCoef[u] * pDataConst[v]; q = pDataCoefImag[u] * pDataConstImag[v];
s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataConst[v] + pDataConstImag[v]);
pDataConst[u] = pDataConst[u] - p + q;
pDataConstImag[u] = pDataConstImag[u] - s + p + q;
}
}
}
}

// 求解调整
for (k = n - 1; k >= 0; k--)
{
if (pnJs[k] != k)
{
for (j = 0; j <= m - 1; j++)
{
u = k * m + j;
v = pnJs[k] * m + j;
p = pDataConst[u];
pDataConst[u] = pDataConst[v];
pDataConst[v] = p;
p = pDataConstImag[u];
pDataConstImag[u] = pDataConstImag[v];
pDataConstImag[v] = p;
}
}
}