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>
        public static void GetRootsetGrad(Matrix mtxLECoef, Matrix mtxLEConst, Matrix mtxResult, double eps)
        {
            int i, k;
            double alpha, beta, d, e;

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

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

            // 构造临时矩阵
            Matrix mtxP = new Matrix(n, 1);
            double[] p = mtxP.GetData();

            double[] pDataCoef = mtxLECoef.GetData();
            double[] pDataConst = mtxLEConst.GetData();

            double[] r = new double[n];

            for (i = 0; i <= n - 1; i++)
            {
                x[i] = 0.0;
                p[i] = pDataConst[i];
                r[i] = pDataConst[i];
            }

            i = 0;
            while (i <= n - 1)
            {
                Matrix mtxS = mtxLECoef.Multiply(mtxP);
                double[] s = mtxS.GetData();

                d = 0.0;
                e = 0.0;
                for (k = 0; k <= n - 1; k++)
                {
                    d = d + p[k] * pDataConst[k];
                    e = e + p[k] * s[k];
                }

                alpha = d / e;
                for (k = 0; k <= n - 1; k++)
                {
                    x[k] = x[k] + alpha * p[k];
                }

                Matrix mtxQ = mtxLECoef.Multiply(mtxResult);
                double[] q = mtxQ.GetData();

                d = 0.0;
                for (k = 0; k <= n - 1; k++)
                {
                    r[k] = pDataConst[k] - q[k];
                    d = d + r[k] * s[k];
                }

                beta = d / e;
                d = 0.0;
                for (k = 0; k <= n - 1; k++)
                {
                    d = d + r[k] * r[k];
                }
                if (Math.Abs(d) < float.Epsilon)
                {
                    break;
                }
                // 满足精度，求解结束
                d = Math.Sqrt(d);
                if (d < eps)
                {
                    break;
                }
                for (k = 0; k <= n - 1; k++)
                {
                    p[k] = r[k] - beta * p[k];
                }
                i = i + 1;
            }
        }
}

}