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>
        /// <param name="nMaxIt">叠加次数</param>
        /// <param name="eps">控制精度</param>
        /// <return>bool 型，方程组求解是否成功</return>
        public static bool GetRootsetMorbid(Matrix mtxLECoef, Matrix mtxLEConst, Matrix mtxResult, int nMaxIt = 60, double eps = 1.0E-7)
        {
            // 方程的阶数
            int n = mtxLECoef.GetNumColumns();

            // 设定迭代次数, 缺省为60
            int i = nMaxIt;

            // 用全选主元高斯消元法求解
            //LEquations leqs = new LEquations(mtxLECoef, mtxLEConst);
            if (!LEquations.GetRootsetGauss(mtxLECoef, mtxLEConst, mtxResult))
            {
                return false;
            }
            double[] x = mtxResult.GetData();

            double q = 1.0 + eps;
            while (q >= eps)
            {
                // 迭代次数已达最大值，仍为求得结果，求解失败
                if (i == 0)
                {
                    return false;
                }

                // 迭代次数减1
                i = i - 1;

                // 矩阵运算
                Matrix mtxE = mtxLECoef.Multiply(mtxResult);
                Matrix mtxR = mtxLEConst.Subtract(mtxE);

                // 用全选主元高斯消元法求解
                //leqs = new LEquations(mtxLECoef, mtxR);
                Matrix mtxRR = new Matrix();
                if (!LEquations.GetRootsetGauss(mtxLECoef, mtxR, mtxRR))
                {
                    return false;
                }
                double[] r = mtxRR.GetData();

                double qq = 0.0;
                for (int k = 0; k <= n - 1; k++)
                {
                    qq = Math.Abs(r[k]) / (1.0 + Math.Abs(x[k] + r[k]));
                    if (qq > q)
                    {
                        q = qq;
                    }
                }

                for (int k = 0; k <= n - 1; k++)
                {
                    x[k] = x[k] + r[k];
                }
            }

            // 求解成功
            return true;
        }
    }
}