gra_math.cs

矩阵运算, 包括全选主元Gusass求逆

using System;
using System.Windows.Forms;

namespace adjust
{
    /*****************************
     * 平差计算中用到的矩阵计算公式
     * 类名   gra_math
     * 
     * matrix_add()        矩阵加法
     * matrix_sub()        矩阵减法
     * matrix_mul()        矩阵乘法
     * matrix_mul_const()  矩阵乘常数
     * matrix_trans()      矩阵转置
     * matrix_inverse()    矩阵求逆
     * 
     *****************************/
    public class gra_math
    {
        /***********************
         * 函数名     matrix_add()
         * 参数：      
         *            A   类型  double型二维数组
         *            B   类型  double型二维数组
         * 返回值     double型二维数组
         * 函数功能   计算输入矩阵A，B的和
         * 
         **************************/

        public static double[,] matrix_add(double[,] A, double[,] B)
        {
            double[,] result;
            result = new double[A.GetLength(0), A.GetLength(1)];
            if (A.GetLength(0) != B.GetLength(0) || A.GetLength(1) != B.GetLength(1))
            {
                MessageBox.Show("矩阵大小不一样,无法相加!", "矩阵运算错误", MessageBoxButtons.OK, MessageBoxIcon.Warning);
                //添加程序运行结束代码  return false;
            }
            else
            {
                for (int i = 0; i < A.GetLength(0); i++)
                {
                    for (int j = 0; j < A.GetLength(1); j++)
                    {
                        result[i, j] = A[i, j] + B[i, j];
                    }
                }
            }

            return result;
        }

        /**************************
         * 函数名      matrix_sub()
         * 参数：
         *            A   类型  double型二维数组
         *            B   类型  double型二维数组
         * 返回值     double型二维数组
         * 函数功能   计算输入矩阵A，B的差
         ***************************/
        public static double[,] matrix_sub(double[,] A, double[,] B)
        {
            double[,] result;
            result = new double[A.GetLength(0), A.GetLength(1)];
            if (A.GetLength(0) != B.GetLength(0) || A.GetLength(1) != B.GetLength(1))
            {
                MessageBox.Show("矩阵大小不一样,无法相减!", "矩阵运算错误", MessageBoxButtons.OK, MessageBoxIcon.Warning);
                //return false;
            }
            else
            {
                for (int i = 0; i < A.GetLength(0); i++)
                {
                    for (int j = 0; j < A.GetLength(1); j++)
                    {
                        result[i, j] = A[i, j] - B[i, j];
                    }
                }
            }

            return result;
        }

        /**************************
        * 函数名     matrix_mul()
        * 参数：
        *            A   类型  double型二维数组
        *            B   类型  double型二维数组
        * 返回值     double型二维数组
        * 函数功能   计算输入矩阵A，B的乘积（A左乘B）
        ***************************/
        public static double[,] matrix_mul(double[,] A, double[,] B)
        {
            double[,] result;
            result = new double[A.GetLength(0), B.GetLength(1)];
            if (A.GetLength(1) != B.GetLength(0))
            {
                MessageBox.Show("矩阵大小不匹配,无法相乘!", "矩阵运算错误", MessageBoxButtons.OK, MessageBoxIcon.Warning);
                //return false;
            }
            else
            {
                for (int j = 0; j < B.GetLength(1); j++)
                {
                    for (int i = 0; i < A.GetLength(0); i++)
                    {
                        double sum = 0.0;
                        for (int k = 0; k < A.GetLength(1); k++)
                        {
                            sum += A[i, k] * B[k, j];
                        }
                        result[i, j] = sum;
                    }
                }
            }
            return result;
        }

        /*************************
         * 函数名    matrix_mul_const()
         * 参数：
         *           A   类型  double型二维数组
         *           c   类型  double
         * 返回值    double型二维数组
         * 函数功能  计算矩阵A与常数c的乘积
         ************************/
        public static double[,] matrix_mul_const(double[,] A, double c)
        {
            double[,] result;
            result = new double[A.GetLength(0), A.GetLength(1)];

            for (int i = 0; i < A.GetLength(0); i++)
            {
                for (int j = 0; j < A.GetLength(1); j++)
                {
                    result[i, j] = A[i, j] * c;
                }
            }
            return result;
        }

        /*************************
        * 函数名    matrix_trans()
        * 参数：
        *           A   类型  double型二维数组
        * 返回值    double型二维数组
        * 函数功能  计算矩阵A的转置矩阵
        ************************/
        public static double[,] matrix_trans(double[,] A)
        {
            double[,] result;
            result = new double[A.GetLength(1), A.GetLength(0)];

            for (int i = 0; i < A.GetLength(1); i++)
            {
                for (int j = 0; j < A.GetLength(0); j++)
                {
                    result[i, j] = A[j, i];
                }
            }

            return result;
        }

        /****************************
         * 函数名    matrix_inverse()
         * 参数：
         *           A   类型  double型二维数组 
         * 返回值    double型二维数组
         * 函数功能  计算矩阵A的逆矩阵， 运用全选主元高斯-约当法求逆算法
         *****************************/
        public static double[,] matrix_inverse(double[,] A)
        {
            double[,] result;
            if (A.GetLength(0) != A.GetLength(1))
            {
                MessageBox.Show("矩阵不为方阵,无法求逆!", "矩阵运算错误", MessageBoxButtons.OK, MessageBoxIcon.Warning);
                //return false;
            }

            int[] IS = new int[A.GetLength(0)];                       //IS[] JS[]为工作数组
            int[] JS = new int[A.GetLength(0)];

            result = new double[A.GetLength(0), A.GetLength(1)];

            int i, j, k;
            double dMax;

            for (i = 0; i < A.GetLength(0); i++)                   //将A[,]的值付给result[,]
            {
                for (j = 0; j < A.GetLength(0); j++)
                {
                    result[i, j] = A[i, j];
                }
            }


            for (k = 0; k < result.GetLength(0); k++)
            {
                dMax = 0.0;
                for (i = k; i < result.GetLength(0); i++)
                {
                    for (j = k; j < result.GetLength(0); j++)
                    {
                        if (Math.Abs(result[i, j]) > dMax)
                        {
                            dMax = Math.Abs(result[i, j]);
                            IS[k] = i;
                            JS[k] = j;
                        }
                    }
                }

                if ((dMax + 1.0) == 1.0)
                {
                    MessageBox.Show("主元太小,无法求逆!", "矩阵运算错误", MessageBoxButtons.OK, MessageBoxIcon.Warning);
                    //return false;
                }

                for (j = 0; j < result.GetLength(0); j++)
                {
                    double t = result[k, j];
                    result[k, j] = result[IS[k], j];
                    result[IS[k], j] = t;
                }

                for (i = 0; i < result.GetLength(0); i++)
                {
                    double t = result[i, k];
                    result[i, k] = result[i, JS[k]];
                    result[i, JS[k]] = t;
                }

                result[k, k] = 1 / result[k, k];

                for (j = 0; j < result.GetLength(0); j++)
                {
                    if (j != k)
                    {
                        result[k, j] = result[k, j] * result[k, k];
                    }
                }

                for (i = 0; i < result.GetLength(0); i++)
                {
                    if (i != k)
                    {
                        for (j = 0; j < result.GetLength(0); j++)
                        {
                            if (j != k)
                            {
                                result[i, j] = result[i, j] - result[i, k] * result[k, j];
                            }
                        }
                    }
                }

                for (i = 0; i < result.GetLength(0); i++)
                {
                    if (i != k)
                    {
                        result[i, k] = -1.0 * result[i, k] * result[k, k];
                    }
                }
            }

            for (k = result.GetLength(0) - 1; k >= 0; k--)
            {
                for (j = 0; j < result.GetLength(0); j++)
                {
                    double t = result[k, j];
                    result[k, j] = result[JS[k], j];
                    result[JS[k], j] = t;
                }

                for (i = 0; i < result.GetLength(0); i++)
                {
                    double t = result[i, k];
                    result[i, k] = result[i, IS[k]];
                    result[i, IS[k]] = t;
                }
            }

            return result;

        }



    }
}