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="eps - 控制精度</param>
        /// <return>bool 型，方程组求解是否成功</return>
        public static bool GetRootsetGaussSeidel(Matrix mtxLECoef, Matrix mtxLEConst, Matrix mtxResult, double eps)
        {
            int i, j, u, v;
            double p, t, s, q;

            // 未知数个数
            int n = mtxLECoef.GetNumColumns();

            // 初始化解向量
            mtxResult.Init(n, 1);
            double[] x = mtxResult.GetData();

            // 系数与常数
            double[] pDataCoef = mtxLECoef.GetData();
            double[] pDataConst = mtxLEConst.GetData();

            // 求解
            for (i = 0; i <= n - 1; i++)
            {
                u = i * n + i;
                p = 0.0;
                x[i] = 0.0;
                for (j = 0; j <= n - 1; j++)
                {
                    if (i != j)
                    {
                        v = i * n + j;
                        p = p + Math.Abs(pDataCoef[v]);
                    }
                }

                if (p >= Math.Abs(pDataCoef[u]))
                {
                    return false;
                }
            }

            // 精度控制
            p = eps + 1.0;
            while (p >= eps)
            {
                p = 0.0;
                for (i = 0; i <= n - 1; i++)
                {
                    t = x[i];
                    s = 0.0;
                    for (j = 0; j <= n - 1; j++)
                    {
                        if (j != i)
                        {
                            s = s + pDataCoef[i * n + j] * x[j];
                        }
                    }
                    x[i] = (pDataConst[i] - s) / pDataCoef[i * n + i];
                    q = Math.Abs(x[i] - t) / (1.0 + Math.Abs(x[i]));
                    if (q > p)
                    {
                        p = q;
                    }
                }
            }

            return true;
        }
 

}

}