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 GetRootsetGaussJordan(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();
            int m = mtxLEConst.GetNumColumns();

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

            // 消元
            u = 1;
            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)//  d + 1.0 == 1.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;
                        }

                        for (j = 0; j <= m - 1; j++)
                        {
                            p = k * m + j;
                            q = nIs * m + j;
                            t = pDataConst[p];
                            pDataConst[p] = pDataConst[q];
                            pDataConst[q] = 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;
                }

                for (j = 0; j <= m - 1; j++)
                {
                    p = k * m + j;
                    pDataConst[p] = pDataConst[p] / d;
                }

                for (j = k + 1; j <= n - 1; j++)
                {
                    for (i = 0; i <= n - 1; i++)
                    {
                        p = i * n + j;
                        if (i != k)
                        {
                            pDataCoef[p] = pDataCoef[p] - pDataCoef[i * n + k] * pDataCoef[k * n + j];
                        }
                    }
                }
                for (j = 0; j <= m - 1; j++)
                {
                    for (i = 0; i <= n - 1; i++)
                    {
                        p = i * m + j;
                        if (i != k)
                        {
                            pDataConst[p] = pDataConst[p] - pDataCoef[i * n + k] * pDataConst[k * m + j];
                        }
                    }
                }
            }

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

            return true;
        }
    }
}