微软公司内部培训程序员资料---求解线性方程组的类

/*
 * 求解线性方程组的类 LEquations
 * 
 * 周长发编制
 */

using System;

namespace MSAlgorithm
{
    public class LEquations
    {
        /// <summary>
        /// 系数矩阵
        /// </summary>
        private Matrix mtxLECoef;

        /// <summary>
        /// 常数矩阵
        /// </summary>
        private Matrix mtxLEConst;

        /// <summary>
        /// 基本构造函数
        /// </summary>
        public LEquations()
        { }

        /// <summary>
        /// 指定系数和常数构造函数
        /// </summary>
        /// <param name="mtxCoef">指定的系数矩阵</param>
        /// <param name="mtxConst">指定的常数矩阵</param>
        public LEquations(Matrix mtxCoef, Matrix mtxConst)
        {
            Init(mtxCoef, mtxConst);
        }

        /// <summary>
        /// 初始化函数
        /// </summary>
        /// <param name="mtxCoef">指定的系数矩阵</param>
        /// <param name="mtxConst">指定的常数矩阵</param>
        /// <returns>bool 型,初始化是否成功</returns>
        public bool Init(Matrix mtxCoef, Matrix mtxConst)
        {
            if (mtxCoef.GetNumRows() != mtxConst.GetNumRows())
                return false;

            mtxLECoef = new Matrix(mtxLECoef);
            mtxConst = new Matrix(mtxConst);

            return true;
        }

        /// <summary>
        /// 获取系数矩阵
        /// </summary>
        /// <returns>Matrix 型,返回系数矩阵</returns>
        public Matrix GetCoefMatrix()
        {
            return mtxLECoef;
        }

        /// <summary>
        /// 获取常数矩阵
        /// </summary>
        /// <returns>Matrix 型,返回常数矩阵</returns>
        public Matrix GetConstMatrix()
        {
            return mtxLEConst;
        }

        /// <summary>
        /// 获取方程个数
        /// </summary>
        /// <returns>int 型,返回方程组方程的个数</returns>
        public int GetNumEquations()
        {
            return GetCoefMatrix().GetNumRows();
        }

        /// <summary>
        /// 获取未知数个数
        /// </summary>
        /// <returns>int 型,返回方程组未知数的个数</returns>
        public int GetNumUnknowns()
        {
            return GetConstMatrix().GetNumColumns();
        }

        /// <summary>
        /// 全选主元高斯消去法
        /// </summary>
        /// <param name="mtxResult">Matrix对象,返回方程组的解</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetGauss(Matrix mtxResult)
        {
            int l, k, i, j, nIs = 0, p, q;
            double d, t;

            //方程组的属性,将常数矩阵赋给解矩阵
            mtxResult.SetValue(mtxLEConst);
            double[] pDataCoef = mtxLECoef.GetData();
            double[] pDataConst = mtxResult.GetData();
            int n = GetNumUnknowns();

            //临时缓冲区,存放列数
            int[] pnJs = new int[n];

            //消元
            l = 1;
            for (k = 0; k <= n - 2; k++)
            {
                d = 0.0;
                for (i = k; i <= n - 1; i++)
                {
                    for (j = k; j <= n - 1; j++)
                    {
                        t = Math.Abs(pDataCoef[i * n + j]);
                        if (t > d)
                        {
                            d = t;
                            pnJs[k] = j;
                            nIs = i;
                        }
                    }
                }

                if (d == 0.0)
                    l = 0;
                else
                {
                    if (pnJs[k] != k)
                    {
                        for (i = 0; i <= n - 1; i++)
                        {
                            p = i * n + k;
                            q = i * n + pnJs[k];
                            t = pDataCoef[p];
                            pDataCoef[p] = pDataCoef[q];
                            pDataCoef[q] = t;
                        }
                    }

                    if (nIs != k)
                    {
                        for (j = k; j <= n - 1; j++)
                        {
                            p = k * n + j;
                            q = nIs * n + j;
                            t = pDataCoef[p];
                            pDataCoef[p] = pDataCoef[q];
                            pDataCoef[q] = t;
                        }

                        t = pDataConst[k];
                        pDataConst[k] = pDataConst[nIs];
                        pDataConst[nIs] = t;
                    }
                }

                //求解失败
                if (l == 0)
                    return false;

                d = pDataCoef[k * n + k];
                for (j = k + 1; j <= n - 1; j++)
                {
                    p = k * n + j;
                    pDataCoef[p] = pDataCoef[p] / d;
                }

                pDataConst[k] = pDataConst[k] / d;
                for (i = k + 1; i <= n - 1; i++)
                {
                    for (j = k + 1; j <= n - 1; j++)
                    {
                        p = i * n + j;
                        pDataCoef[p] = pDataCoef[p] - pDataCoef[i * n + k] * pDataCoef[k * n + j];
                    }

                    pDataConst[i] = pDataConst[i] - pDataCoef[i * n + k] * pDataConst[k];
                }
            }

            //求解失败
            d = pDataCoef[(n - 1) * n + n - 1];
            if (d == 0.0)
            {
                return false;
            }

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

            //调整解的位置
            pnJs[n - 1] = n - 1;
            for (k = n - 1; k >= 0; k--)
            {
                if (pnJs[k] != k)
                {
                    t = pDataConst[k];
                    pDataConst[k] = pDataConst[pnJs[k]];
                    pDataConst[pnJs[k]] = t;
                }
            }

            return true;
        }

        /// <summary>
        /// 全选主元高斯-约当消去法
        /// </summary>
        /// <param name="mtxResult">Matrix对象,返回方程组的解</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetGaussJordan(Matrix mtxResult)
        {
            int l, k, i, j, nIs = 0, p, q;
            double d, t;

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

            //临时缓冲区,存放变换的列数
            int[] pnJs = new int[n];

            //消元
            l = 1;
            for (k = 0; k <= n - 1; k++)
            {
                d = 0.0;
                for (i = k; i <= n - 1; i++)
                {
                    for (j = k; j <= n - 1; j++)
                    {
                        t = Math.Abs(pDataCoef[i * n + j]);
                        if (t > d)
                        {
                            d = t;
                            pnJs[k] = j;
                            nIs = i;
                        }
                    }
                }

                if (d + 1.0 == 1.0)
                    l = 0;
                else
                {
                    if (pnJs[k] != k)
                    {
                        for (i = 0; i <= n - 1; i++)
                        {
                            p = i * n + k;
                            q = i * n + pnJs[k];
                            t = pDataCoef[p];
                            pDataCoef[p] = pDataCoef[q];
                            pDataCoef[q] = t;
                        }
                    }

                    if (nIs != k)
                    {
                        for (j = k; j <= n - 1; j++)
                        {
                            p = k * n + j;
                            q = nIs * n + j;
                            t = pDataCoef[p];
                            pDataCoef[p] = pDataCoef[q];
                            pDataCoef[q] = t;
                        }

                        for (j = 0; j <= m - 1; j++)
                        {
                            p = k * m + j;
                            q = nIs * m + j;
                            t = pDataConst[p];
                            pDataConst[p] = pDataConst[q];
                            pDataConst[q] = t;
                        }
                    }
                }

                //求解失败
                if (l == 0)
                    return false;

                d = pDataCoef[k * n + k];
                for (j = k + 1; j <= n - 1; j++)
                {
                    p = k * n + j;
                    pDataCoef[p] = pDataCoef[p] / d;
                }

                for (j = 0; j <= m - 1; j++)
                {
                    p = k * m + j;
                    pDataConst[p] = pDataConst[p] / d;
                }

                for (j = k + 1; j <= n - 1; j++)
                {
                    for (i = 0; i <= n - 1; i++)
                    {
                        p = i * n + j;
                        if (i != k)
                            pDataCoef[p] = pDataCoef[p] - pDataCoef[i * n + k] * pDataCoef[k * n + j];
                    }
                }

                for (j = 0; j <= m - 1; j++)
                {
                    for (i = 0; i <= n - 1; i++)
                    {
                        p = i * m + j;
                        if (i != k)
                            pDataConst[p] = pDataConst[p] - pDataCoef[i * n + k] * pDataConst[k * m + j];
                    }
                }
            }

            //调整
            for (k = n - 1; k >= 0; k--)
            {
                if (pnJs[k] != k)
                {
                    for (j = 0; j <= m - 1; j++)
                    {
                        p = k * m + j;
                        q = pnJs[k] * m + j;
                        t = pDataConst[p];
                        pDataConst[p] = pDataConst[q];
                        pDataConst[q] = t;
                    }
                }
            }

            return true;
        }

        /// <summary>
        /// 复系数方程组的全选主元高斯消去法
        /// </summary>
        /// <param name="mtxCoefImag">系数矩阵的虚部矩阵</param>
        /// <param name="mtxConstImag">常数矩阵的虚部矩阵</param>
        /// <param name="mtxResult">Matrix对象,返回方程组解矩阵的实部矩阵</param>
        /// <param name="mtxResultImag">Matrix对象,返回方程组解矩阵的虚部矩阵</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetGauss(Matrix mtxCoefImag, Matrix mtxConstImag, Matrix mtxResult, Matrix mtxResultImag)
        {
            int l, k, i, j, nIs = 0, u, v;
            double p, q, s, d;

            //方程组的属性,将常数矩阵赋给解矩阵
            mtxResult.SetValue(mtxLEConst);
            mtxResultImag.SetValue(mtxConstImag);
            double[] pDataCoef = mtxLECoef.GetData();
            double[] pDataConst = mtxResult.GetData();
            double[] pDataCoefImag = mtxCoefImag.GetData();
            double[] pDataConstImag = mtxResultImag.GetData();
            int n = GetNumUnknowns();
            int m = mtxLEConst.GetNumColumns();

            //临时缓冲区,存放变换的列数
            int[] pnJs = new int[n];

            //消元
            for (k = 0; k <= n - 2; k++)
            {
                d = 0.0;
                for (i = k; i <= n - 1; i++)
                {
                    for (j = k; j <= n - 1; j++)
                    {
                        u = i * n + j;
                        p = pDataCoef[u] * pDataCoef[u] + pDataCoefImag[u] * pDataCoefImag[u];
                        if (p > d)
                        {
                            d = p;
                            pnJs[k] = j;
                            nIs = i;
                        }
                    }
                }

                //求解失败
                if (d == 0.0)
                    return false;

                if (nIs != k)
                {
                    for (j = k; j <= n - 1; j++)
                    {
                        u = k * n + j;
                        v = nIs * n + j;
                        p = pDataCoef[u];
                        pDataCoef[u] = pDataCoef[v];
                        pDataCoef[v] = p;
                        p = pDataCoefImag[u];
                        pDataCoefImag[u] = pDataCoefImag[v];
                        pDataCoefImag[v] = p;
                    }

                    p = pDataConst[k];
                    pDataConst[k] = pDataConst[nIs];
                    pDataConst[nIs] = p;
                    p = pDataConstImag[k];
                    pDataConstImag[k] = pDataConstImag[nIs];
                    pDataConstImag[nIs] = p;
                }

                if (pnJs[k] != k)
                {
                    for (i = 0; i <= n - 1; i++)
                    {
                        u = i * n + k;
                        v = i * n + pnJs[k];
                        p = pDataCoef[u];
                        pDataCoef[u] = pDataCoef[v];
                        pDataCoef[v] = p;
                        p = pDataCoefImag[u];
                        pDataCoefImag[u] = pDataCoefImag[v];
                        pDataCoefImag[v] = p;
                    }
                }

                v = k * n + k;
                for (j = k + 1; j <= n - 1; j++)
                {
                    u = k * n + j;
                    p = pDataCoef[u] * pDataCoef[v];
                    q = -pDataCoefImag[u] * pDataCoefImag[v];
                    s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataCoef[u] + pDataCoefImag[u]);
                    pDataCoef[u] = (p - q) / d;
                    pDataCoefImag[u] = (s - p - q) / d;
                }

                p = pDataConst[k] * pDataCoef[v];
                q = -pDataConstImag[k] * pDataCoefImag[v];
                s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataConst[k] + pDataConstImag[k]);
                pDataConst[k] = (p - q) / d;
                pDataConstImag[k] = (s - p - q) / d;

                for (i = k + 1; i <= n - 1; i++)
                {
                    u = i * n + k;
                    for (j = k + 1; j <= n - 1; j++)
                    {
                        v = k * n + j;
                        l = i * n + j;
                        p = pDataCoef[u] * pDataCoef[v];
                        q = pDataCoefImag[u] * pDataCoefImag[v];
                        s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataCoef[v] + pDataCoefImag[v]);
                        pDataCoef[l] = pDataCoef[l] - p + q;
                        pDataCoefImag[l] = pDataCoefImag[l] - s + p + q;
                    }

                    p = pDataCoef[u] * pDataConst[k];
                    q = pDataCoefImag[u] * pDataConstImag[k];
                    s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataConst[k] + pDataConstImag[k]);
                    pDataConst[i] = pDataConst[i] - p + q;
                    pDataConstImag[i] = pDataConstImag[i] - s + p + q;
                }
            }

            u = (n - 1) * n + n - 1;
            d = pDataCoef[u] * pDataCoef[u] + pDataCoefImag[u] * pDataCoefImag[u];

            //求解失败
            if (d == 0.0)
                return false;

            //求解
            p = pDataCoef[u] * pDataConst[n - 1]; q = -pDataCoefImag[u] * pDataConstImag[n - 1];
            s = (pDataCoef[u] - pDataCoefImag[u]) * (pDataConst[n - 1] + pDataConstImag[n - 1]);
            pDataConst[n - 1] = (p - q) / d; pDataConstImag[n - 1] = (s - p - q) / d;

            for (i = n - 2; i >= 0; i--)
            {
                for (j = i + 1; j <= n - 1; j++)
                {
                    u = i * n + j;
                    p = pDataCoef[u] * pDataConst[j];
                    q = pDataCoefImag[u] * pDataConstImag[j];
                    s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataConst[j] + pDataConstImag[j]);
                    pDataConst[i] = pDataConst[i] - p + q;
                    pDataConstImag[i] = pDataConstImag[i] - s + p + q;
                }
            }

            //调整位置
            pnJs[n - 1] = n - 1;
            for (k = n - 1; k >= 0; k--)
            {
                if (pnJs[k] != k)
                {
                    p = pDataConst[k];
                    pDataConst[k] = pDataConst[pnJs[k]];
                    pDataConst[pnJs[k]] = p;
                    p = pDataConstImag[k];
                    pDataConstImag[k] = pDataConstImag[pnJs[k]];
                    pDataConstImag[pnJs[k]] = p;
                }
            }

            return true;
        }

        /// <summary>
        /// 复系数方程组的全选主元高斯-约当消去法
        /// </summary>
        /// <param name="mtxCoefImag">系数矩阵的虚部矩阵</param>
        /// <param name="mtxConstImag">常数矩阵的虚部矩阵</param>
        /// <param name="mtxResult">Matrix对象,返回方程组解矩阵的实部矩阵</param>
        /// <param name="mtxResultImag">Matrix对象,返回方程组解矩阵的虚部矩阵</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetGaussJordan(Matrix mtxCoefImag, Matrix mtxConstImag, Matrix mtxResult, Matrix mtxResultImag)
        {
            int l, k, i, j, nIs = 0, u, v;
            double p, q, s, d;

            //方程组的属性,将常数矩阵赋给解矩阵
            mtxResult.SetValue(mtxLEConst);
            mtxResultImag.SetValue(mtxConstImag);
            double[] pDataCoef = mtxLECoef.GetData();
            double[] pDataConst = mtxResult.GetData();
            double[] pDataCoefImag = mtxCoefImag.GetData();
            double[] pDataConstImag = mtxResultImag.GetData();
            int n = GetNumUnknowns();
            int m = mtxLEConst.GetNumColumns();

            //临时缓冲区,存放变换的列数
            int[] pnJs = new int[n];

            //消元
            for (k = 0; k <= n - 1; k++)
            {
                d = 0.0;
                for (i = k; i <= n - 1; i++)
                {
                    for (j = k; j <= n - 1; j++)
                    {
                        u = i * n + j;
                        p = pDataCoef[u] * pDataCoef[u] + pDataCoefImag[u] * pDataCoefImag[u];
                        if (p > d)
                        {
                            d = p;
                            pnJs[k] = j;
                            nIs = i;
                        }
                    }
                }

                //求解失败
                if (d == 0.0)
                    return false;

                if (nIs != k)
                {
                    for (j = k; j <= n - 1; j++)
                    {
                        u = k * n + j;
                        v = nIs * n + j;
                        p = pDataCoef[u];
                        pDataCoef[u] = pDataCoef[v];
                        pDataCoef[v] = p;
                        p = pDataCoefImag[u];
                        pDataCoefImag[u] = pDataCoefImag[v];
                        pDataCoefImag[v] = p;
                    }

                    for (j = 0; j <= m - 1; j++)
                    {
                        u = k * m + j;
                        v = nIs * m + j;
                        p = pDataConst[u];
                        pDataConst[u] = pDataConst[v];
                        pDataConst[v] = p;
                        p = pDataConstImag[u];
                        pDataConstImag[u] = pDataConstImag[v];
                        pDataConstImag[v] = p;
                    }
                }

                if (pnJs[k] != k)
                {
                    for (i = 0; i <= n - 1; i++)
                    {
                        u = i * n + k;
                        v = i * n + pnJs[k];
                        p = pDataCoef[u];
                        pDataCoef[u] = pDataCoef[v];
                        pDataCoef[v] = p;
                        p = pDataCoefImag[u];
                        pDataCoefImag[u] = pDataCoefImag[v];
                        pDataCoefImag[v] = p;
                    }
                }

                v = k * n + k;
                for (j = k + 1; j <= n - 1; j++)
                {
                    u = k * n + j;
                    p = pDataCoef[u] * pDataCoef[v];
                    q = -pDataCoefImag[u] * pDataCoefImag[v];
                    s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataCoef[u] + pDataCoefImag[u]);
                    pDataCoef[u] = (p - q) / d;
                    pDataCoefImag[u] = (s - p - q) / d;
                }

                for (j = 0; j <= m - 1; j++)
                {
                    u = k * m + j;
                    p = pDataConst[u] * pDataCoef[v];
                    q = -pDataConstImag[u] * pDataCoefImag[v];
                    s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataConst[u] + pDataConstImag[u]);
                    pDataConst[u] = (p - q) / d;
                    pDataConstImag[u] = (s - p - q) / d;
                }

                for (i = 0; i <= n - 1; i++)
                {
                    if (i != k)
                    {
                        u = i * n + k;
                        for (j = k + 1; j <= n - 1; j++)
                        {
                            v = k * n + j;
                            l = i * n + j;
                            p = pDataCoef[u] * pDataCoef[v];
                            q = pDataCoefImag[u] * pDataCoefImag[v];
                            s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataCoef[v] + pDataCoefImag[v]);
                            pDataCoef[l] = pDataCoef[l] - p + q;
                            pDataCoefImag[l] = pDataCoefImag[l] - s + p + q;
                        }

                        for (j = 0; j <= m - 1; j++)
                        {
                            l = i * m + j;
                            v = k * m + j;
                            p = pDataCoef[u] * pDataConst[v]; q = pDataCoefImag[u] * pDataConstImag[v];
                            s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataConst[v] + pDataConstImag[v]);
                            pDataConst[l] = pDataConst[l] - p + q;
                            pDataConstImag[l] = pDataConstImag[l] - s + p + q;
                        }
                    }
                }
            }

            //求解调整
            for (k = n - 1; k >= 0; k--)
            {
                if (pnJs[k] != k)
                {
                    for (j = 0; j <= m - 1; j++)
                    {
                        u = k * m + j;
                        v = pnJs[k] * m + j;
                        p = pDataConst[u];
                        pDataConst[u] = pDataConst[v];
                        pDataConst[v] = p;
                        p = pDataConstImag[u];
                        pDataConstImag[u] = pDataConstImag[v];
                        pDataConstImag[v] = p;
                    }
                }
            }

            return true;
        }

        /// <summary>
        /// 求解三对角线方程组的追赶法
        /// </summary>
        /// <param name="mtxResult">Matrix对象,返回方程组解矩阵</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetTriDiagonal(Matrix mtxResult)
        {
            int k, j;
            double s;

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

            int n = GetNumUnknowns();
            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) + 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 (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;
        }

        /// <summary>
        /// 一般带型方程组的求解
        /// </summary>
        /// <param name="nBandWidth">带宽</param>
        /// <param name="mtxResult">Matrix对象,返回方程组解矩阵</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetBand(int nBandWidth, Matrix mtxResult)
        {
            int ls, k, i, j, nis = 0, u, v;
            double p, t;

            //带宽必须为奇数
            if ((nBandWidth - 1) % 2 != 0)
                return false;

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

            //方程组属性
            int m = mtxLEConst.GetNumColumns();
            int n = GetNumUnknowns();
            if (mtxLECoef.GetNumRows() != n)
                return false;

            //带宽数组:带型矩阵的有效数据
            double[] pBandData = new double[nBandWidth * n];

            //半带宽
            ls = (nBandWidth - 1) / 2;

            //构造带宽数组
            for (i = 0; i < n; ++i)
            {
                j = 0;
                for (k = Math.Max(0, i - ls); k < Math.Max(0, i - ls) + nBandWidth; ++k)
                {
                    if (k < n)
                        pBandData[i * nBandWidth + j++] = mtxLECoef.GetElement(i, k);
                    else
                        pBandData[i * nBandWidth + j++] = 0;
                }
            }

            //求解
            for (k = 0; k <= n - 2; k++)
            {
                p = 0.0;
                for (i = k; i <= ls; i++)
                {
                    t = Math.Abs(pBandData[i * nBandWidth]);
                    if (t > p)
                    {
                        p = t;
                        nis = i;
                    }
                }

                if (p == 0.0)
                    return false;

                for (j = 0; j <= m - 1; j++)
                {
                    u = k * m + j;
                    v = nis * m + j;
                    t = pDataConst[u];
                    pDataConst[u] = pDataConst[v];
                    pDataConst[v] = t;
                }

                for (j = 0; j <= nBandWidth - 1; j++)
                {
                    u = k * nBandWidth + j;
                    v = nis * nBandWidth + j;
                    t = pBandData[u];
                    pBandData[u] = pBandData[v];
                    pBandData[v] = t;
                }

                for (j = 0; j <= m - 1; j++)
                {
                    u = k * m + j;
                    pDataConst[u] = pDataConst[u] / pBandData[k * nBandWidth];
                }

                for (j = 1; j <= nBandWidth - 1; j++)
                {
                    u = k * nBandWidth + j;
                    pBandData[u] = pBandData[u] / pBandData[k * nBandWidth];
                }

                for (i = k + 1; i <= ls; i++)
                {
                    t = pBandData[i * nBandWidth];
                    for (j = 0; j <= m - 1; j++)
                    {
                        u = i * m + j;
                        v = k * m + j;
                        pDataConst[u] = pDataConst[u] - t * pDataConst[v];
                    }

                    for (j = 1; j <= nBandWidth - 1; j++)
                    {
                        u = i * nBandWidth + j;
                        v = k * nBandWidth + j;
                        pBandData[u - 1] = pBandData[u] - t * pBandData[v];
                    }

                    u = i * nBandWidth + nBandWidth - 1; pBandData[u] = 0.0;
                }

                if (ls != (n - 1))
                    ls = ls + 1;
            }

            p = pBandData[(n - 1) * nBandWidth];
            if (p == 0.0)
            {
                return false;
            }

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

            ls = 1;
            for (i = n - 2; i >= 0; i--)
            {
                for (k = 0; k <= m - 1; k++)
                {
                    u = i * m + k;
                    for (j = 1; j <= ls; j++)
                    {
                        v = i * nBandWidth + j;
                        nis = (i + j) * m + k;
                        pDataConst[u] = pDataConst[u] - pBandData[v] * pDataConst[nis];
                    }
                }

                if (ls != (nBandWidth - 1))
                    ls = ls + 1;
            }

            return true;
        }

        /// <summary>
        /// 求解对称方程组的分解法
        /// </summary>
        /// <param name="mtxResult">CMatrix引用对象,返回方程组解矩阵</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetDjn(Matrix mtxResult)
        {
            int i, j, l, k, u, v, w, k1, k2, k3;
            double p;

            //方程组属性,将常数矩阵赋给解矩阵
            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;

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

            for (i = 1; i <= n - 2; i++)
            {
                u = i * n + i;
                for (j = 1; j <= i; j++)
                {
                    v = i * n + j - 1;
                    l = (j - 1) * n + j - 1;
                    pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[v] * pDataCoef[l];
                }

                p = pDataCoef[u];
                if (p == 0.0)
                    return false;

                for (k = i + 1; k <= n - 1; k++)
                {
                    u = k * n + i;
                    for (j = 1; j <= i; j++)
                    {
                        v = k * n + j - 1;
                        l = i * n + j - 1;
                        w = (j - 1) * n + j - 1;
                        pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[l] * pDataCoef[w];
                    }

                    pDataCoef[u] = pDataCoef[u] / p;
                }
            }

            u = n * n - 1;
            for (j = 1; j <= n - 1; j++)
            {
                v = (n - 1) * n + j - 1;
                w = (j - 1) * n + j - 1;
                pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[v] * pDataCoef[w];
            }

            p = pDataCoef[u];
            if (p == 0.0)
                return false;

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

            for (i = 1; i <= n - 1; i++)
            {
                u = (i - 1) * n + i - 1;
                for (j = i; j <= n - 1; j++)
                {
                    v = (i - 1) * n + j;
                    w = j * n + i - 1;
                    pDataCoef[v] = pDataCoef[u] * pDataCoef[w];
                }
            }

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

                    pDataConst[u] = pDataConst[u] / pDataCoef[k3 * n + k3];
                }
            }

            return true;
        }

        /// <summary>
        /// 求解对称正定方程组的平方根法
        /// </summary>
        /// <param name="mtxResult">CMatrix引用对象,返回方程组解矩阵</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetCholesky(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;
        }

        /// <summary>
        /// 求解大型稀疏方程组的全选主元高斯-约去消去法
        /// </summary>
        /// <param name="mtxResult">CMatrix引用对象,返回方程组解矩阵</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetGgje(Matrix mtxResult)
        {
            int i, j, k, nIs = 0, u, v;
            double d, t;

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

            //临时缓冲区,存放变换的列数
            int[] pnJs = new int[n];

            //消元
            for (k = 0; k <= n - 1; k++)
            {
                d = 0.0;
                for (i = k; i <= n - 1; i++)
                {
                    for (j = k; j <= n - 1; j++)
                    {
                        t = Math.Abs(pDataCoef[i * n + j]);
                        if (t > d)
                        {
                            d = t;
                            pnJs[k] = j;
                            nIs = i;
                        }
                    }
                }

                if (d == 0.0)
                    return false;

                if (nIs != k)
                {
                    for (j = k; j <= n - 1; j++)
                    {
                        u = k * n + j;
                        v = nIs * n + j;
                        t = pDataCoef[u];
                        pDataCoef[u] = pDataCoef[v];
                        pDataCoef[v] = t;
                    }

                    t = pDataConst[k];
                    pDataConst[k] = pDataConst[nIs];
                    pDataConst[nIs] = t;
                }

                if (pnJs[k] != k)
                {
                    for (i = 0; i <= n - 1; i++)
                    {
                        u = i * n + k;
                        v = i * n + pnJs[k];
                        t = pDataCoef[u];
                        pDataCoef[u] = pDataCoef[v];
                        pDataCoef[v] = t;
                    }
                }

                t = pDataCoef[k * n + k];
                for (j = k + 1; j <= n - 1; j++)
                {
                    u = k * n + j;
                    if (pDataCoef[u] != 0.0)
                        pDataCoef[u] = pDataCoef[u] / t;
                }

                pDataConst[k] = pDataConst[k] / t;
                for (j = k + 1; j <= n - 1; j++)
                {
                    u = k * n + j;
                    if (pDataCoef[u] != 0.0)
                    {
                        for (i = 0; i <= n - 1; i++)
                        {
                            v = i * n + k;
                            if ((i != k) && (pDataCoef[v] != 0.0))
                            {
                                nIs = i * n + j;
                                pDataCoef[nIs] = pDataCoef[nIs] - pDataCoef[v] * pDataCoef[u];
                            }
                        }
                    }
                }

                for (i = 0; i <= n - 1; i++)
                {
                    u = i * n + k;
                    if ((i != k) && (pDataCoef[u] != 0.0))
                        pDataConst[i] = pDataConst[i] - pDataCoef[u] * pDataConst[k];
                }
            }

            //调整
            for (k = n - 1; k >= 0; k--)
            {
                if (k != pnJs[k])
                {
                    t = pDataConst[k];
                    pDataConst[k] = pDataConst[pnJs[k]];
                    pDataConst[pnJs[k]] = t;
                }
            }

            return true;
        }

        /// <summary>
        /// 求解托伯利兹方程组的列文逊方法
        /// </summary>
        /// <param name="mtxResult">CMatrix引用对象,返回方程组解矩阵</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetTlvs(Matrix mtxResult)
        {
            int i, j, k;
            double a, beta, q, c, h;

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

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

            //常数数组
            double[] pDataConst = mtxLEConst.GetData();

            //建立T数组
            double[] t = new double[n];

            //构造T数组
            for (i = 0; i < n; ++i)
                t[i] = mtxLECoef.GetElement(0, i);

            //临时数组
            double[] s = new double[n];
            double[] y = new double[n];

            //非托伯利兹方程组,不能用本方法求解
            a = t[0];
            if (a == 0.0)
                return false;

            //列文逊方法求解
            y[0] = 1.0;
            x[0] = pDataConst[0] / a;
            for (k = 1; k <= n - 1; k++)
            {
                beta = 0.0;
                q = 0.0;
                for (j = 0; j <= k - 1; j++)
                {
                    beta = beta + y[j] * t[j + 1];
                    q = q + x[j] * t[k - j];
                }

                if (a == 0.0)
                    return false;

                c = -beta / a;
                s[0] = c * y[k - 1];
                y[k] = y[k - 1];
                if (k != 1)
                {
                    for (i = 1; i <= k - 1; i++)
                        s[i] = y[i - 1] + c * y[k - i - 1];
                }

                a = a + c * beta;
                if (a == 0.0)
                    return false;

                h = (pDataConst[k] - q) / a;
                for (i = 0; i <= k - 1; i++)
                {
                    x[i] = x[i] + h * s[i];
                    y[i] = s[i];
                }

                x[k] = h * y[k];
            }

            return true;
        }

        /// <summary>
        /// 高斯-赛德尔迭代法
        /// </summary>
        /// <param name="mtxResult">CMatrix引用对象,返回方程组解矩阵</param>
        /// <param name="eps">控制精度</param>
        /// <returns>bool 型,方程组求解是否成功</returns>
        public bool GetRootsetGaussSeidel(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;
        }

        /// <summary>
        /// 求解对称正定方程组的共轭梯度法
        /// </summary>
        /// <param name="mtxResult">CMatrix引用对象,返回方程组解矩阵</param>
        /// <param name="eps">控制精度</param>
        public void GetRootsetGrad(Matrix mtxResult, double eps)
        {
            int i, k;
            double alpha, beta, d, e;

            //未知数个数
            int n = GetNumUnknowns();

            //初始化解向量
            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];

                //满足精度,求解结束
                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;
            }
        }

        /// <summary>
        /// 求解线性最小二乘问题的豪斯荷尔德变换法
        /// </summary>
        /// <param name="mtxResult">Matrix对象,返回方程组解矩阵</param>
        /// <param name="mtxQ">Matrix对象,返回豪斯荷尔德变换的Q矩阵</param>
        /// <param name="mtxR">Matrix对象,返回豪斯荷尔德变换的R矩阵</param>
        /// <returns>bool型,方程组求解是否成功</returns>
        public bool GetRootsetMqr(Matrix mtxResult, Matrix mtxQ, Matrix mtxR)
        {
            int i, j;
            double d;

            //方程组的方程数和未知数个数
            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 (!mtxR.SplitQR(mtxQ))
                return false;

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

            //求解
			for (i=0; i<=n-1; i++)
			{ 
				d=0.0;
				for (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 (i = n - 2; i >= 0; i--)
            {
                d = 0.0;
                for (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;
        }

        /// <summary>
        /// 求解线性最小二乘问题的广义逆法
        /// </summary>
        /// <param name="mtxResult">Matrix对象,返回方程组解矩阵</param>
        /// <param name="mtxAP">Matrix对象,返回系数矩阵的广义逆矩阵</param>
        /// <param name="mtxU">Matrix对象,返回U矩阵</param>
        /// <param name="mtxV">Matrix对象,返回V矩阵</param>
        /// <param name="eps">控制精度</param>
        /// <returns>bool型,方程组求解是否成功</returns>
        public bool GetRootsetGinv(Matrix mtxResult, Matrix mtxAP, Matrix mtxU, Matrix mtxV, double eps)
        {
            int i, j;

            //方程个数和未知数个数
            int m = mtxLECoef.GetNumRows();
            int n = mtxLECoef.GetNumColumns();

            //初始化解向量
            mtxResult.Init(n, 1);

            double[] pDataConst = mtxLEConst.GetData();
            double[] x = mtxResult.GetData();

            //临时矩阵
            Matrix mtxA = new Matrix(mtxLECoef);

            //求广义逆矩阵
            if (!mtxA.InvertUV(mtxAP, mtxU, mtxV, eps))
                return false;

            double[] pAPData = mtxAP.GetData();

            //求解
            for (i = 0; i <= n - 1; i++)
            {
                x[i] = 0.0;
                for (j = 0; j <= m - 1; j++)
                    x[i] = x[i] + pAPData[i * m + j] * pDataConst[j];
            }

            return true;
        }

        /// <summary>
        /// Matrix对象,返回方程组解矩阵
        /// </summary>
        /// <param name="mtxResult">Matrix对象,返回方程组解矩阵</param>
        /// <param name="nMaxIt">叠加次数</param>
        /// <param name="eps">控制精度</param>
        /// <returns>bool型,方程组求解是否成功</returns>
        public bool GetRootsetMorbid(Matrix mtxResult, int nMaxIt, double eps)
        {
            int i, k;
            double q, qq;

            //方程的阶数
            int n = GetNumUnknowns();

            //设定迭代次数, 缺省为60
            i = nMaxIt;

            //用全选主元高斯消元法求解
            LEquations leqs = new LEquations(mtxLECoef, mtxLEConst);
            if (!leqs.GetRootsetGauss(mtxResult))
                return false;
            double[] x = mtxResult.GetData();

            q = 1.0 + eps;
            while (q >= eps)
            {
                // 迭代次数已达最大值,仍为求得结果,求解失败
                if (i == 0)
                    return false;

                // 迭代次数减1
                i = i - 1;

                // 矩阵运算
                Matrix mtxE = mtxLECoef.Multiply(mtxResult);
                Matrix mtxR = mtxLEConst.Substract(mtxE);

                // 用全选主元高斯消元法求解
                leqs = new LEquations(mtxLECoef, mtxR);
                Matrix mtxRR = new Matrix();
                if (!leqs.GetRootsetGauss(mtxRR))
                    return false;

                double[] r = mtxRR.GetData();

                q = 0.0;
                for (k = 0; k <= n - 1; k++)
                {
                    qq = Math.Abs(r[k]) / (1.0 + Math.Abs(x[k] + r[k]));
                    if (qq > q)
                        q = qq;
                }

                for (k = 0; k <= n - 1; k++)
                    x[k] = x[k] + r[k];
            }

            return true;
        }
    }
}

注:本代码为原作者所有,本人只是摘录,如有侵犯作者,请与本人联系

Email:359194966@qq.com
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值