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 GetRootsetTriDiagonal(Matrix mtxLECoef, Matrix mtxLEConst, Matrix mtxResult)
        {
            int k, j;
            double s;

            // 将常数矩阵赋给解矩阵
            mtxResult.SetValue(mtxLEConst);
            double[] pDataConst = mtxResult.GetData();

            int n = mtxLECoef.GetNumColumns();
            if (mtxLECoef.GetNumRows() != n)
            {
                return false;
            }

            // 为系数矩阵对角线数组分配内存
            double[] pDiagData = new double[3 * n - 2];

            // 构造系数矩阵对角线元素数组
            k = j = 0;
            if (k == 0)
            {
                pDiagData[j++] = mtxLECoef.GetElement(k, k);
                pDiagData[j++] = mtxLECoef.GetElement(k, k + 1);
            }
            for (k = 1; k < n - 1; ++k)
            {
                pDiagData[j++] = mtxLECoef.GetElement(k, k - 1);
                pDiagData[j++] = mtxLECoef.GetElement(k, k);
                pDiagData[j++] = mtxLECoef.GetElement(k, k + 1);
            }
            if (k == n - 1)
            {
                pDiagData[j++] = mtxLECoef.GetElement(k, k - 1);
                pDiagData[j++] = mtxLECoef.GetElement(k, k);
            }

            // 求解
            for (k = 0; k <= n - 2; k++)
            {
                j = 3 * k;
                s = pDiagData[j];

                // 求解失败
                if (Math.Abs(s) < float.Epsilon) //  Math.Abs(s) + 1.0 == 1.0)
                {
                    return false;
                }

                pDiagData[j + 1] = pDiagData[j + 1] / s;
                pDataConst[k] = pDataConst[k] / s;
                pDiagData[j + 3] = pDiagData[j + 3] - pDiagData[j + 2] * pDiagData[j + 1];
                pDataConst[k + 1] = pDataConst[k + 1] - pDiagData[j + 2] * pDataConst[k];
            }

            s = pDiagData[3 * n - 3];
            if (Math.Abs(s) < float.Epsilon)//  s == 0.0)
            {
                return false;
            }

            // 调整
            pDataConst[n - 1] = pDataConst[n - 1] / s;
            for (k = n - 2; k >= 0; k--)
            {
                pDataConst[k] = pDataConst[k] - pDiagData[3 * k + 1] * pDataConst[k + 1];
            }
            return true;
        }
 

    }
}