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 GetRootsetDjn(Matrix mtxLECoef, Matrix mtxLEConst, Matrix mtxResult)
        {
            int i, j, r, k, u, v, w, k1, k2, k3;
            double p;

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

            // 非对称系数矩阵，不能用本方法求解
            if (Math.Abs(pDataCoef[0]) < float.Epsilon)//  pDataCoef[0] == 0.0)
            {
                return false;
            }
            for (i = 1; i <= n - 1; i++)
            {
                u = i * n;
                pDataCoef[u] = pDataCoef[u] / pDataCoef[0];
            }

            for (i = 1; i <= n - 2; i++)
            {
                u = i * n + i;
                for (j = 1; j <= i; j++)
                {
                    v = i * n + j - 1;
                    r = (j - 1) * n + j - 1;
                    pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[v] * pDataCoef[r];
                }

                p = pDataCoef[u];
                if (Math.Abs(p) < float.Epsilon)
                {
                    return false;
                }
                for (k = i + 1; k <= n - 1; k++)
                {
                    u = k * n + i;
                    for (j = 1; j <= i; j++)
                    {
                        v = k * n + j - 1;
                        r = i * n + j - 1;
                        w = (j - 1) * n + j - 1;
                        pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[r] * pDataCoef[w];
                    }

                    pDataCoef[u] = pDataCoef[u] / p;
                }
            }

            u = n * n - 1;
            for (j = 1; j <= n - 1; j++)
            {
                v = (n - 1) * n + j - 1;
                w = (j - 1) * n + j - 1;
                pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[v] * pDataCoef[w];
            }

            p = pDataCoef[u];
            if (Math.Abs(p) < float.Epsilon)
            {
                return false;
            }
            for (j = 0; j <= m - 1; j++)
            {
                for (i = 1; i <= n - 1; i++)
                {
                    u = i * m + j;
                    for (k = 1; k <= i; k++)
                    {
                        v = i * n + k - 1;
                        w = (k - 1) * m + j;
                        pDataConst[u] = pDataConst[u] - pDataCoef[v] * pDataConst[w];
                    }
                }
            }

            for (i = 1; i <= n - 1; i++)
            {
                u = (i - 1) * n + i - 1;
                for (j = i; j <= n - 1; j++)
                {
                    v = (i - 1) * n + j;
                    w = j * n + i - 1;
                    pDataCoef[v] = pDataCoef[u] * pDataCoef[w];
                }
            }

            for (j = 0; j <= m - 1; j++)
            {
                u = (n - 1) * m + j;
                pDataConst[u] = pDataConst[u] / p;
                for (k = 1; k <= n - 1; k++)
                {
                    k1 = n - k;
                    k3 = k1 - 1;
                    u = k3 * m + j;
                    for (k2 = k1; k2 <= n - 1; k2++)
                    {
                        v = k3 * n + k2;
                        w = k2 * m + j;
                        pDataConst[u] = pDataConst[u] - pDataCoef[v] * pDataConst[w];
                    }

                    pDataConst[u] = pDataConst[u] / pDataCoef[k3 * n + k3];
                }
            }

            return true;
        }
 

    }
}