using System;

namespace Zhou.CSharp.Algorithm
{
    /// <summary>
    /// 求解线性方程组的类 LEquations
    /// 原作 周长发
    /// 改编 深度混淆
    /// </summary>
    public static class LEquations
    {
        /// <summary>
        /// 全选主元高斯消去法
        /// </summary>
        /// <param name="mtxLECoef">指定的系数矩阵</param>
        /// <param name="mtxLEConst">指定的常数矩阵</param>
        /// <param name="mtxResult">Matrix对象，返回方程组的解</param>
        /// <return>bool 型，方程组求解是否成功</return>
        public static bool GetRootsetGauss(Matrix mtxLECoef, Matrix mtxLEConst, Matrix mtxResult)
        {
            int u, k, i, j, nIs = 0, p, q;
            double d, t;

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

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

            // 消元
            u = 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 (Math.Abs(d) < float.Epsilon)//  d== 0.0)
                {
                    u = 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 (u == 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 (Math.Abs(d) < float.Epsilon)// 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;
        }

    }
 }