matrix.cs

在Visual 2008环境下

    public static Matrix MatrixDotMul(Matrix ma, Matrix mb)
    {
        if ((ma == null) || (mb == null))
            return null;

        if ((ma._Row != mb._Row) || (ma._Col != mb._Col))
            return null;

        int row = ma._Row, col = ma._Col;
        double[,] data_ma = ma._MatrixData;
        double[,] data_mb = mb._MatrixData;
        double[,] resultData = new double[row, col];

        for (int i = 0; i < row; i++)
            for (int j = 0; j < col; j++)
                resultData[i, j] = data_ma[i, j] * data_mb[i, j];

        return new Matrix(resultData);
    }

    /// <summary>
    /// 两矩阵点除
    /// </summary>
    public static Matrix MatrixDotDivide(Matrix ma, Matrix mb)
    {
        Matrix mb_temp = mb.ElementDividedBy1();
        return MatrixDotMul(ma, mb_temp);
    }
    #endregion //对两个矩阵元素一对一的操作

    #region //矩阵与标量之间的操作
    /// <summary>
    /// 矩阵与标量相加
    /// </summary>
    public static Matrix operator +(Matrix ma, double t)
    {
        if (ma == null)
            return null;

        int row = ma._Row, col = ma._Col;
        double[,] resultData = new double[row, col];

        for (int i = 0; i < row; i++)
            for (int j = 0; j < col; j++)
                resultData[i, j] = ma._MatrixData[i, j] + t;

        return new Matrix(resultData);
    }

    /// <summary>
    /// 标量与矩阵相加
    /// </summary>
    public static Matrix operator +(double t, Matrix ma)
    {
        return ma + t;
    }

    /// <summary>
    /// 矩阵与标量相减
    /// </summary>
    public static Matrix operator -(Matrix ma, double t)
    {
        return ma + (-t);
    }

    /// <summary>
    /// 标量与矩阵相减
    /// </summary>
    public static Matrix operator -(double t, Matrix ma)
    {
        Matrix ma_temp = ma.ElementNegative();
        return ma_temp + t;
    }

    /// <summary>
    /// 矩阵与标量相乘
    /// </summary>
    public static Matrix operator *(Matrix ma, double t)
    {
        if (ma == null)
            return null;

        int row = ma._Row, col = ma._Col;
        double[,] resultData = new double[row, col];

        for (int i = 0; i < row; i++)
            for (int j = 0; j < col; j++)
                resultData[i, j] = ma._MatrixData[i, j] * t;

        return new Matrix(resultData);
    }

    /// <summary>
    /// 标量与矩阵相乘
    /// </summary>
    public static Matrix operator *(double t, Matrix ma)
    {
        return ma * t;
    }

    /// <summary>
    /// 矩阵与标量相除
    /// </summary>
    public static Matrix operator /(Matrix ma, double t)
    {
        return ma * (1 / t);
    }

    /// <summary>
    /// 标量与矩阵相除
    /// </summary>
    public static Matrix operator /(double t, Matrix ma)
    {
        Matrix ma_temp = ma.ElementDividedBy1();
        return ma_temp * t;
    }
    #endregion //矩阵与标量之间的操作

    /// <summary>
    /// 两矩阵相乘
    /// </summary>
    public static Matrix operator *(Matrix ma, Matrix mb)
    {
        if ((ma == null) || (mb == null))
            return null;

        if (ma._Col != mb._Row)
            return null;

        double temp_sum = 0.0;
        int row_ma = ma._Row, col_mb = mb._Col;
        int temp_num = ma._Col;
        double[,] data_ma = ma._MatrixData;
        double[,] data_mb = mb._MatrixData;
        double[,] resultData = new double[row_ma, col_mb];

        for (int i = 0; i < row_ma; i++)
        {
            for (int j = 0; j < col_mb; j++)
            {
                for (int k = 0; k < temp_num; k++)
                    temp_sum += data_ma[i, k] * data_mb[k, j];
                resultData[i, j] = temp_sum;
                temp_sum = 0.0;
            }
        }

        return new Matrix(resultData);
    }
    #endregion //算术运算

    #region //统计
    /// <summary>
    /// 统计矩阵中零元素的个数
    /// </summary>

    /// <summary>
    /// 统计矩阵中位于某区域的元素的个数
    /// </summary>

    /// <summary>
    /// 计算矩阵元素之和
    /// </summary>
    public double Sum()
    {
        if (this == null)
            return 0;

        double sum = 0.0;
        for (int i = 0; i < col; i++)
            for (int j = 0; j < row; j++)
                sum += matrixData[j, i];

        return sum;
    }

    /// <summary>
    /// 计算矩阵元素绝对值之和
    /// </summary>
    public double AbsSum()
    {
        if (this == null)
            return 0;

        double absSum = 0.0;
        for (int i = 0; i < col; i++)
            for (int j = 0; j < row; j++)
                absSum += Math.Abs(matrixData[j, i]);

        return absSum;
    }

    /// <summary>
    /// 计算矩阵元素的平均值
    /// </summary>


    /// <summary>
    /// 查找矩阵元素的最大和最小值
    /// </summary>
    public double[] FindMaxAndMin()
    {
        if (this == null)
            return null;

        double max = 0.0, min = 0.0;

        for (int i = 0; i < row; i++)
        {
            for (int j = 0; j < col; j++)
            {
                if ((i == 0 && j == 0))
                    max = min = matrixData[i, j];

                if (max < matrixData[i, j])
                    max = matrixData[i, j];
                else if (min > matrixData[i, j])
                    min = matrixData[i, j];
            }
        }
        double[] resultData = new double[2] { max, min };

        return resultData;
    }

    #endregion //统计

    #region //线性代数
    public void SetIdentity(double t)
    {
        if (this == null || this.row != this.col)
            return;

        for (int i = 0; i < this.row; i++)
            this.matrixData[i, i] = t;
    }

    /// <summary>
    /// 矩阵转置
    /// </summary>
    public Matrix MatrixTranspose()
    {
        if (this == null)
            return null;

        double[,] resultData = new double[col, row];
        for (int i = 0; i < col; i++)
            for (int j = 0; j < row; j++)
                resultData[i, j] = matrixData[j, i];

        return new Matrix(resultData);
    }

    /// <summary>
    /// 求矩阵逆——参见《常用算法程序集(C语言描述)第三版P40》
    /// 用全选主元高斯—约当消去法求n阶实矩阵的逆
    /// </summary>
    public bool MatrixInvert()
    {
        int i, j, k;
        double d, p;
        int n = this._Row;
        double[,] a = this.matrixData;
        int[] Is = new int[n];
        int[] Js = 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++)
                {
                    p = Math.Abs(a[i, j]);
                    if (p > d)
                    {
                        d = p;
                        Is[k] = i;
                        Js[k] = j;
                    }
                }
            }
            if ((d + 1.0) == 1.0)
                return false;

            #region 为了数值计算的稳定性，整个求逆过程需要全选主元
            if (Is[k] != k)
            {
                for (j = 0; j <= n - 1; j++)
                {
                    p = a[k, j];
                    a[k, j] = a[Is[k], j];
                    a[Is[k], j] = p;
                }
            }

            if (Js[k] != k)
            {
                for (i = 0; i <= n - 1; i++)
                {
                    p = a[i, k];
                    a[i, k] = a[i, Js[k]];
                    a[i, Js[k]] = a[i, k];
                }
            }
            #endregion 为了数值计算的稳定性，整个求逆过程需要全选主元

            #region 归一化计算
            a[k, k] = 1.0 / a[k, k];
            for (j = 0; j <= n - 1; j++)
                if (j != k)
                    a[k, j] = a[k, j] * a[k, k];
            #endregion 归一化计算

            #region 消元计算
            for (i = 0; i <= n - 1; i++)
                if (i != k)
                    for (j = 0; j <= n - 1; j++)
                        if (j != k)
                            a[i, j] = a[i, j] - a[i, k] * a[k, j];

            for (i = 0; i <= n - 1; i++)
                if (i != k)
                    a[i, k] = -a[i, k] * a[k, k];
            #endregion 消元计算
        }

        #region 经全选主元后，矩阵求逆的计算过程是稳定的。但是最后需要对结果进行恢复
        for (k = n - 1; k >= 0; k--)
        {
            if (Js[k] != k)
            {
                for (j = 0; j <= n - 1; j++)
                {
                    p = a[k, j];
                    a[k, j] = a[Js[k], j];
                    a[Js[k], j] = p;
                }
            }

            if (Is[k] != k)
            {
                for (i = 0; i <= n - 1; i++)
                {
                    p = a[i, k];
                    a[i, k] = a[i, Is[k]];
                    a[i, Is[k]] = p;