C#，码海拾贝（28）——求解“对称正定方程组”的“平方根法”之C#源代码

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 GetRootsetCholesky(Matrix mtxLECoef, Matrix mtxLEConst, Matrix mtxResult)
        {
            int i, j, k, u, v;

            // 方程组属性，将常数矩阵赋给解矩阵
            Matrix mtxCoef = new Matrix(mtxLECoef);
            mtxResult.SetValue(mtxLEConst);
            int n = mtxCoef.GetNumColumns();
            int m = mtxResult.GetNumColumns();
            double[] pDataCoef = mtxCoef.GetData();
            double[] pDataConst = mtxResult.GetData();

            // 非对称正定系数矩阵，不能用本方法求解
            if (pDataCoef[0] <= 0.0)
            {
                return false;
            }

            pDataCoef[0] = Math.Sqrt(pDataCoef[0]);
            for (j = 1; j <= n - 1; j++)
            {
                pDataCoef[j] = pDataCoef[j] / pDataCoef[0];
            }
            for (i = 1; i <= n - 1; i++)
            {
                u = i * n + i;
                for (j = 1; j <= i; j++)
                {
                    v = (j - 1) * n + i;
                    pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[v];
                }

                if (pDataCoef[u] <= 0.0)
                {
                    return false;
                }
                pDataCoef[u] = Math.Sqrt(pDataCoef[u]);
                if (i != (n - 1))
                {
                    for (j = i + 1; j <= n - 1; j++)
                    {
                        v = i * n + j;
                        for (k = 1; k <= i; k++)
                        {
                            pDataCoef[v] = pDataCoef[v] - pDataCoef[(k - 1) * n + i] * pDataCoef[(k - 1) * n + j];
                        }
                        pDataCoef[v] = pDataCoef[v] / pDataCoef[u];
                    }
                }
            }

            for (j = 0; j <= m - 1; j++)
            {
                pDataConst[j] = pDataConst[j] / pDataCoef[0];
                for (i = 1; i <= n - 1; i++)
                {
                    u = i * n + i;
                    v = i * m + j;
                    for (k = 1; k <= i; k++)
                    {
                        pDataConst[v] = pDataConst[v] - pDataCoef[(k - 1) * n + i] * pDataConst[(k - 1) * m + j];
                    }
                    pDataConst[v] = pDataConst[v] / pDataCoef[u];
                }
            }

            for (j = 0; j <= m - 1; j++)
            {
                u = (n - 1) * m + j;
                pDataConst[u] = pDataConst[u] / pDataCoef[n * n - 1];
                for (k = n - 1; k >= 1; k--)
                {
                    u = (k - 1) * m + j;
                    for (i = k; i <= n - 1; i++)
                    {
                        v = (k - 1) * n + i;
                        pDataConst[u] = pDataConst[u] - pDataCoef[v] * pDataConst[i * m + j];
                    }

                    v = (k - 1) * n + k - 1;
                    pDataConst[u] = pDataConst[u] / pDataCoef[v];
                }
            }

            return true;
        }
 

    }
}