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="mtxQ">Matrix对象，返回豪斯荷尔德变换的Q矩阵</param>
        /// <param name="mtxR">Matrix对象，返回豪斯荷尔德变换的R矩阵</param>
        /// <return>bool 型，方程组求解是否成功</return>
        public static bool GetRootsetMqr(Matrix mtxLECoef, Matrix mtxLEConst, Matrix mtxResult, Matrix mtxQ, Matrix mtxR)
        {
            // 方程组的方程数和未知数个数
            int m = mtxLECoef.GetNumRows();
            int n = mtxLECoef.GetNumColumns();
            // 奇异方程组
            if (m < n)
            {
                return false;
            }

            // 将解向量初始化为常数向量
            mtxResult.SetValue(mtxLEConst);
            double[] pDataConst = mtxResult.GetData();

            // 构造临时矩阵，用于QR分解
            mtxR.SetValue(mtxLECoef);
            double[] pDataCoef = mtxR.GetData();

            // QR分解
            if (!Matrix.SplitQR(mtxR, mtxQ))
            {
                return false;
            }

            // 临时缓冲区
            double[] c = new double[n];
            double[] q = mtxQ.GetData();

            // 求解
            for (int i = 0; i <= n - 1; i++)
            {
                double d = 0.0;
                for (int j = 0; j <= m - 1; j++)
                {
                    d = d + q[j * m + i] * pDataConst[j];
                }
                c[i] = d;
            }

            pDataConst[n - 1] = c[n - 1] / pDataCoef[n * n - 1];
            for (int i = n - 2; i >= 0; i--)
            {
                double d = 0.0;
                for (int j = i + 1; j <= n - 1; j++)
                {
                    d = d + pDataCoef[i * n + j] * pDataConst[j];
                }
                pDataConst[i] = (c[i] - d) / pDataCoef[i * n + i];
            }

            return true;
        }
}

}