C# 高斯列主元消去法

Question

public bool GetRootsetGauss(Matrix mtxResult)
{ 
int l,k,i,j,p,q;
            int J=0;
double d,t; // 方程组的属性，将常数矩阵赋给解矩阵
mtxResult.SetValue(mtxLEConst);
double[] pDataCoef = mtxLECoef.GetData();
double[] pDataConst = mtxResult.GetData();
int n = GetNumUnknowns(); // 消元
l=1;
for (k=0;k<=n-2;k++)
{
                d = 0;
                for (j = k; j < n - 1; j++)
                {
                    t = Math.Abs(pDataCoef[n*j+k]);
                    if (t > d)
                    { d = t; J = j; }
                }                if (d == 0)
                    l = 0;
                else
                {
                    for (i = 0; i < n; i++)
                    {
                        double temp1 = pDataCoef[n *J + i];
                        pDataCoef[n * J + i] = pDataCoef[n * k + i];
                        pDataCoef[n * k + i] = temp1; 
                    }
                    double temp2 = pDataCoef[k];
                    pDataCoef[k]=pDataCoef[J];
                    pDataCoef[J] = temp2;
                }
        
// 求解失败
if (l==0)
{ 
return false;
}
        
d=pDataCoef[k*n+k];
for (j=k+1;j<=n-1;j++)
{ 
p=k*n+j; 
pDataCoef[p]=pDataCoef[p]/d;
}
        
pDataConst[k]=pDataConst[k]/d;
for (i=k+1;i<=n-1;i++)
{ 
for (j=k+1;j<=n-1;j++)
{ 
p=i*n+j;
pDataCoef[p]=pDataCoef[p]-pDataCoef[i*n+k]*pDataCoef[k*n+j];
}
            
pDataConst[i]=pDataConst[i]-pDataCoef[i*n+k]*pDataConst[k];
}
}
    
// 求解失败
d=pDataCoef[(n-1)*n+n-1];
if (d == 0.0)
{ 
return false;
} // 求解
pDataConst[n-1]=pDataConst[n-1]/d;
for (i=n-2;i>=0;i--)
{ 
t=0.0;
for (j=i+1;j<=n-1;j++)
t=t+pDataCoef[i*n+j]*pDataConst[j];
                pDataConst[i] = (pDataConst[i] - t);
}  
return true;
}
Martix为 矩阵类
有两个构造函数
但结果不正确
为何？

Answer 1

明白了 应该是这样的        public bool GetRootsetGauss2(Matrix mtxResult)
        {
            int l, k, i, j, nIs = 0, p, q;
            double d, t;            // 方程组的属性，将常数矩阵赋给解矩阵
            mtxResult.SetValue(mtxLEConst);
            double[] pDataCoef = mtxLECoef.GetData();
            double[] pDataConst = mtxResult.GetData();
            int n = GetNumUnknowns();            // 临时缓冲区，存放列数
            int[] pnJs = new int[n];            // 消元
            l = 1;
            for (k = 0; k <= n - 2; 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 (d == 0.0)
                    l = 0;
                else
                {
                    if (pnJs[k] != k)
                    {
                        for (i = 0; i <= n - 1; i++)
                        {
                            p = i * n + k;
                            q = i * n + pnJs[k];
                            t = pDataCoef[p];
                            pDataCoef[p] = pDataCoef[q];
                            pDataCoef[q] = t;
                        }
                    }                    if (nIs != k)
                    {
                        for (j = k; j <= n - 1; j++)
                        {
                            p = k * n + j;
                            q = nIs * n + j;
                            t = pDataCoef[p];
                            pDataCoef[p] = pDataCoef[q];
                            pDataCoef[q] = t;
                        }                        t = pDataConst[k];
                        pDataConst[k] = pDataConst[nIs];
                        pDataConst[nIs] = t;
                    }
                }                // 求解失败
                if (l == 0)
                {
                    return false;
                }                d = pDataCoef[k * n + k];
                for (j = k + 1; j <= n - 1; j++)
                {
                    p = k * n + j;
                    pDataCoef[p] = pDataCoef[p] / d;
                }                pDataConst[k] = pDataConst[k] / d;
                for (i = k + 1; i <= n - 1; i++)
                {
                    for (j = k + 1; j <= n - 1; j++)
                    {
                        p = i * n + j;
                        pDataCoef[p] = pDataCoef[p] - pDataCoef[i * n + k] * pDataCoef[k * n + j];
                    }                    pDataConst[i] = pDataConst[i] - pDataCoef[i * n + k] * pDataConst[k];
                }
            }            // 求解失败
            d = pDataCoef[(n - 1) * n + n - 1];
            if (d == 0.0)
            {
                return false;
            }            // 求解
            pDataConst[n - 1] = pDataConst[n - 1] / d;
            for (i = n - 2; i >= 0; i--)
            {
                t = 0.0;
                for (j = i + 1; j <= n - 1; j++)
                    t = t + pDataCoef[i * n + j] * pDataConst[j];
                pDataConst[i] = pDataConst[i] - t;
            }            // 调整解的位置
            pnJs[n - 1] = n - 1;
            for (k = n - 1; k >= 0; k--)
            {
                if (pnJs[k] != k)
                {
                    t = pDataConst[k];
                    pDataConst[k] = pDataConst[pnJs[k]];
                    pDataConst[pnJs[k]] = t;
                }
            }            return true;
        }