矩阵变换

在一些算法中需要用到矩阵，自然就需要用到矩阵的一些操作，比如行变换、列变换、最简式、求矩阵的秩等，下面是实现的代码

public class Matrix
    {
        #region 属性
        /// <summary>
        /// 行数
        /// </summary>
        public int RowCount
        {
            get { return _rowCount; }
        }

        /// <summary>
        /// 列数
        /// </summary>
        public int ColumnCount
        {
            get { return _columnCount; }
        }
        #endregion

        #region 内部变量
        /// <summary>
        /// 行数
        /// </summary>
        /// <remarks>Using this instead of the RowCount property to speed up calculating
        /// a matrix index in the data array.</remarks>
        private readonly int _rowCount;

        /// <summary>
        /// 列数
        /// </summary>
        /// <remarks>Using this instead of the ColumnCount property to speed up calculating
        /// a matrix index in the data array.</remarks>
        private readonly int _columnCount;

        /// <summary>
        /// 矩阵数据
        /// </summary>
        private readonly double[] _data;

        /// <summary>
        /// 矩阵的元素总数
        /// </summary>
        private readonly int _count;
        #endregion

        #region 构造函数
        /// <summary>
        /// 矩阵构造函数
        /// </summary>
        /// <param name="rowCount">行数</param>
        /// <param name="columnCount">列数</param>
        public Matrix(int rowCount, int columnCount)
        {
            if (rowCount <= 0 || columnCount <= 0)
            {
                throw new Exception("矩阵的行数和列数必须大于0");
            }
            _data = new double[rowCount * columnCount];
            _rowCount = rowCount;
            _columnCount = columnCount;
            _count = _rowCount * _columnCount;
        }
        #endregion

        #region 矩阵元素的取值和赋值
        /// <summary>
        /// 得到指定位置的矩阵元素(索引从零开始)
        /// </summary>
        public double At(int row, int column)
        {
            if (!IsLegalIndex(row, column))
            {
                return 0xFFFFFFFF;
            }
            return _data[(column * _rowCount) + row];
        }

        /// <summary>
        ///给指定位置的矩阵元素赋值(索引从零开始)
        /// </summary>
        public void At(int row, int column, double value)
        {
            if (!IsLegalIndex(row, column))
            {
                return;
            }
            _data[(column * _rowCount) + row] = value;
        }

        #region AtRow
        /// <summary>
        /// 给指定的行赋值
        /// </summary>
        /// <param name="row"></param>
        /// <param name="values"></param>
        public void AtRow(int row, params double[] values)
        {
            //for (int column = 0; column < _columnCount; column++)
            //{
            //    At(row, column, values[column]);
            //}
            AtRow(row, values, 0, ColumnCount);
        }

        /// <summary>
        /// 给指定的行赋值
        /// </summary>
        /// <param name="row"></param>
        /// <param name="values"></param>
        /// <param name="valuesStartIndex">values数组的起始索引，通常是values的长度大于矩阵的列数时需要指定.</param>
        /// <param name="valuesCount">values从起始索引算起需要赋值的个数.</param>
        public void AtRow(int row, double[] values, int valuesStartIndex, int valuesCount)
        {
            for (int column = 0; column < valuesCount; column++)
            {
                At(row, column, values[valuesStartIndex + column]);
            }
        }

        /// <summary>
        /// 给指定的行赋值
        /// </summary>
        /// <param name="row"></param>
        /// <param name="matrix">将指定矩阵的对应行的数据赋值到当前矩阵的行</param>
        public void AtRow(int row, Matrix matrix)
        {
            AtRow(row, matrix, row);
        }

        /// <summary>
        /// 给指定的行赋值
        /// </summary>
        /// <param name="row">目标矩阵的行</param>
        /// <param name="rowSource">来源矩阵的行</param>
        /// <param name="matrixSource">将指定矩阵的对应行的数据赋值到当前矩阵的行</param>
        public void AtRow(int row, Matrix matrixSource, int rowSource)
        {
            if (this.ColumnCount != matrixSource.ColumnCount)
            {
                throw new Exception("指定矩阵的列数与当前矩阵的列数不相等");
            }
            for (int column = 0; column < ColumnCount; column++)
            {
                At(row, column, matrixSource.At(rowSource, column));
            }
        }
        #endregion

        /// <summary>
        /// 是否为合法的索引
        /// </summary>
        /// <param name="row"></param>
        /// <param name="column"></param>
        /// <returns></returns>
        private bool IsLegalIndex(int row, int column)
        {
            if (row < 0 || column < 0)
            {
                throw new Exception("行列索引必须大于或等于0");
            }

            if ((column * _rowCount) + row > _count - 1)
            {
                throw new Exception("行列索引超出范围");
            }

            return true;
        }
        #endregion

        #region Copy Function
        public Matrix Copy()
        {
            Matrix target = new Matrix(this._rowCount, this._columnCount);
            Array.Copy(this._data, target._data, _count);
            return target;
        }
        #endregion

        #region Clear
        /// <summary>
        /// 清空矩阵.所有元素都置为0.
        /// </summary>
        public void Clear()
        {
            Array.Clear(_data, 0, _data.Length);
        }
        #endregion

        #region IsMostSimpleRow
        /// <summary>
        /// 是否为最简行.
        /// 即除了指定列和最后一列为非零外，其他列都为零.
        /// </summary>
        /// <param name="row"></param>
        /// <param name="notZeroColumn"></param>
        /// <returns></returns>
        public bool IsMostSimpleRow(int row, int notZeroColumn)
        {
            if (At(row, notZeroColumn) == 0 || At(row, ColumnCount - 1) == 0)
            {//指定列和最后一列为零
                return false;
            }

            for (int column = 0; column < ColumnCount; column++)
            {
                if (column == notZeroColumn || column == ColumnCount - 1)
                {
                    continue;
                }

                if (At(row, column) != 0)
                {//除了指定列和最后一列为非零外，其他列都为零
                    return false;
                }
            }

            return true;
        }
        #endregion

        #region Rank Function
        /// <summary>
        /// 计算矩阵的秩
        /// </summary>
        /// <returns></returns>
        public int Rank()
        {
            //matrix为空则直接默认已经是最简形式
            if (_count == 0)
            {
                return 0;
            }

            //复制一个matrix到martrix_copy，之后因计算需要改动矩阵时并不改动matrix本身
            Matrix matrix_copy = this.Copy();
           
            matrix_copy.TransformMostSimple();
            //行最简矩阵的秩即为所求
            int rank = matrix_copy.RankOfMostSimpleMatrix();
            return rank;
        }

        #region SortByLeftNotZeroPosition
        /// <summary>
        /// 排序（按左侧最前非零位位置自上而下升序排列）
        /// </summary>
        /// <param name="matrix">矩阵</param>
        private void SortByLeftNotZeroPosition()
        {
            //统计每行第一个非零元素的出现位置
            int[] counter = new int[_rowCount];
            for (int r = 0; r < _rowCount; r++)
            {
                for (int c = 0; c < _columnCount; c++)
                {
                    if (At(r, c) == 0)
                    {
                        counter[r]++;
                    }
                    else
                    {
                        break;
                    }
                }
            }

            //按每行非零元素的出现位置升序排列
            for (int r = 0; r < _rowCount; r++)
            {
                for (int j = r; j < _rowCount; j++)
                {
                    if (counter[r] > counter[j])
                    {
                        ExchangeRow(r, j);
                    }
                }
            }
        }
        #endregion

        #region IsMostSimple
        /// <summary>
        /// 判断矩阵是否变换到最简形式（非零行数达到最少）
        /// </summary>
        /// <param name="matrix"></param>
        /// <returns>true:</returns>
        private bool IsMostSimple()
        {
            //统计每行第一个非零元素的出现位置
            int[] counter = new int[_rowCount];
            for (int r = 0; r < _rowCount; r++)
            {
                for (int c = 0; c < _columnCount; c++)
                {
                    if (At(r, c) == 0)
                    {
                        counter[r]++;
                    }
                    else break;
                }
            }

            //后面行的非零元素出现位置必须在前面行的后面，全零行除外
            for (int i = 1; i < counter.Length; i++)
            {
                if (counter[i] <= counter[i - 1] && counter[i] != _columnCount)
                {
                    return false;
                }
            }

            return true;
        }
        #endregion

        #region ElementaryTrasform Function
        /// <summary>
        /// 行初等变换（左侧最前非零位位置最靠前的行，只保留一个）
        /// </summary>
        /// <param name="matrix">矩阵</param>
        private void ElementaryTrasform()
        {
            //统计每行第一个非零元素的出现位置,数值从1开始.
            int[] firstNotZeroPosArr = new int[_rowCount];
            for (int r = 0; r < _rowCount; r++)
            {
                for (int c = 0; c < _columnCount; c++)
                {
                    if (At(r, c) == 0)
                    {
                        firstNotZeroPosArr[r]++;
                    }
                    else
                    {
                        break;
                    }
                }
            }

            for (int row = 1; row < firstNotZeroPosArr.Length; row++)
            {
                if (firstNotZeroPosArr[row] == firstNotZeroPosArr[row - 1] && firstNotZeroPosArr[row] != _columnCount)
                {//该行的非零位置与上面一行的非零位置相同，并且不是最后一列.

                    //上面一行非零位置的值
                    double upRowNotZeroValue = At(row - 1, firstNotZeroPosArr[row - 1]);
                    //当前行非零位置的值
                    double currentRowNotZeroValue = At(row, firstNotZeroPosArr[row]);

                    At(row, firstNotZeroPosArr[row], 0);//将当前位置的值设为0.
                    for (int j = firstNotZeroPosArr[row] + 1; j < _columnCount; j++)
                    {
                        //上面一行非零位置的下一个位置的值
                        double upRowNotZeroNextPosValue = At(row - 1, j);
                        double value = At(row, j) - (upRowNotZeroNextPosValue * currentRowNotZeroValue / upRowNotZeroValue);
                        At(row, j, value);
                    }
                    break;
                }
            }
        }

        #endregion

        #region AssignZeroAlmost
        /// <summary>
        /// 将和0非常接近的数字视为0
        /// </summary>
        /// <param name="matrix"></param>
        private void AssignZeroAlmost()
        {
            for (int i = 0; i < _count; i++)
            {
                if (Math.Abs(_data[i]) <= 0.00001)
                {
                    _data[i] = 0;
                }
            }
        }
        #endregion

        #region RankOfMostSimpleMatrix Function
        /// <summary>
        /// 计算行最简矩阵的秩
        /// </summary>
        /// <param name="matrix"></param>
        /// <returns></returns>
        public int RankOfMostSimpleMatrix()
        {
            int rank = -1;
            bool isAllZero = true;
            for (int r = 0; r < _rowCount; r++)
            {
                isAllZero = true;

                //查看当前行有没有0
                for (int j = 0; j < _columnCount; j++)
                {
                    if (At(r, j) != 0)
                    {
                        isAllZero = false;
                        break;
                    }
                }

                //若第i行全为0，则矩阵的秩为i
                if (isAllZero)
                {
                    rank = r;
                    break;
                }
            }
            //满秩矩阵的情况
            if (rank == -1)
            {
                rank = _rowCount;
            }

            return rank;
        }
        #endregion

        #endregion

        #region TransformMostSimple Function
        /// <summary>
        /// 变换为最简矩阵
        /// </summary>
        public void TransformMostSimple()
        {
            //先以最左侧非零项的位置进行行排序
            SortByLeftNotZeroPosition();

            //循环化简矩阵
            while (!IsMostSimple())
            {
                ElementaryTrasform();
                SortByLeftNotZeroPosition();
            }

            //过于趋近0的项，视作0，减小误差
            AssignZeroAlmost();
        }
        #endregion

        #region ExchangeRow Function
        /// <summary>
        /// 交换矩阵的行
        /// </summary>
        /// <param name="sourceRow">源行</param>
        /// <param name="targetRow">目标行</param>
        public void ExchangeRow(int sourceRow, int targetRow)
        {
            double sourceTemp = 0;
            double targetTemp = 0;
            for (int c = 0; c < _columnCount; c++)
            {
                sourceTemp = At(sourceRow, c);
                targetTemp = At(targetRow, c);
                At(sourceRow, c, targetTemp);
                At(targetRow, c, sourceTemp);
            }
        }
        #endregion

        #region ExchangeColumn Function
        /// <summary>
        /// 交换矩阵的列
        /// </summary>
        /// <param name="sourceColumn">源行</param>
        /// <param name="targetColumn">目标行</param>
        public void ExchangeColumn(int sourceColumn, int targetColumn)
        {
            double sourceTemp = 0;
            double targetTemp = 0;
            for (int row = 0; row < RowCount; row++)
            {
                sourceTemp = At(row, sourceColumn);
                targetTemp = At(row, targetColumn);
                At(row, sourceColumn, targetTemp);
                At(row, targetColumn, sourceTemp);
            }
        }
        #endregion        

        #region OfRowArray Function
        /// <summary>
        /// 得到行数组
        /// </summary>
        /// <param name="rowIndex"></param>
        /// <returns></returns>
        public double[] OfRowArray(int rowIndex)
        {
            double[] rowArray = new double[_columnCount];
            for (int c = 0; c < _columnCount; c++)
            {
                rowArray[c] = At(rowIndex, c);
            }
            return rowArray;
        }
        #endregion

        #region OfColumnArray Function
        /// <summary>
        /// 得到列数组
        /// </summary>
        /// <param name="rowIndex"></param>
        /// <returns></returns>
        public double[] OfColumnArray(int columnIndex)
        {
            double[] columnArray = new double[_rowCount];
            for (int r = 0; r < _rowCount; r++)
            {
                columnArray[r] = At(r, columnIndex);
            }
            return columnArray;
        }
        #endregion

        #region ToString
        public override string ToString()
        {
            //return base.ToString();
            StringBuilder sb = new StringBuilder();
            for (int r = 0; r < _rowCount; r++)
            {
                for (int c = 0; c < _columnCount; c++)
                {
                    //sb.Append(String.Format("{0:000.000},", At(r, c)));
                    //sb.Append(String.Format("{0:00.0},", At(r, c)));
                    sb.Append(At(r, c) + ",");
                }
                sb.Append(";" + Environment.NewLine);
            }
            return sb.ToString();
        }
        #endregion

        #region TransformRow
        /// <summary>
        /// 行变换(row1=row1*factor1-row2*factor2)
        /// </summary>
        /// <param name="row1"></param>
        /// <param name="row2"></param>
        /// <param name="factor1">行1乘以的系数</param>
        /// <param name="factor2">行2乘以的系数</param>
        public void TransformRow(int row1, int row2, double factor2)
        {
            double value1 = 0;
            double value2 = 0;
            for (int column = 0; column < ColumnCount; column++)
            {
                value1 = At(row1, column);
                value2 = At(row2, column) * factor2;
                At(row1, column, value1 - value2);
            }
        }
        #endregion

        #region TransformColumn
        /// <summary>
        /// 列变换(column1=column1*factor1-column2*factor2)
        /// </summary>
        /// <param name="column1"></param>
        /// <param name="column2"></param>
        /// <param name="factor1">列1乘以的系数</param>
        /// <param name="factor2">列2乘以的系数</param>
        public void TransformColumn(int column1, int column2, double factor2)
        {
            double value1 = 0;
            double value2 = 0;
            for (int row = 0; row < RowCount; row++)
            {
                value1 = At(row, column1);
                value2 = At(row, column2) * factor2;
                At(row, column1, value1 - value2);
            }
        }
        #endregion

        #region Assign
        /// <summary>
        /// 将源矩阵的值赋入当前矩阵中
        /// </summary>
        /// <param name="matrixSource"></param>
        public void Assign(Matrix matrixSource)
        {
            for (int row = 0; row < matrixSource.RowCount; row++)
            {
                AtRow(row, matrixSource);
            }
        }

        /// <summary>
        /// 将源矩阵的值赋入当前矩阵中
        /// </summary>
        /// <param name="matrixSource"></param>
        /// <param name="skipRowOfSource">源矩阵中需要跳过的行</param>
        public void Assign(Matrix matrixSource, int skipRowOfSource)
        {
            int rowNew = 0;
            for (int row = 0; row < matrixSource.RowCount; row++)
            {
                if (row == skipRowOfSource)
                {
                    continue;
                }
                AtRow(rowNew, matrixSource, row);
                rowNew++;
            }
        }
        #endregion

        #region EqualeValue
        public bool EqualeValue(Matrix compareMatrix)
        {
            if (compareMatrix == null
                || compareMatrix.RowCount != this.RowCount
                || compareMatrix.ColumnCount != this.ColumnCount)
            {
                return false;
            }

            for (int row = 0; row < RowCount; row++)
            {
                for (int column = 0; column < ColumnCount; column++)
                {
                    if (At(row, column) != compareMatrix.At(row, column))
                    {
                        return false;
                    }
                }
            }
            return true;
        }
        #endregion

        #region 重载运行符
        public static Matrix operator *(Matrix source, double factor)
        {
            for (int r = 0; r < source.RowCount; r++)
            {
                for (int c = 0; c < source.ColumnCount; c++)
                {
                    source.At(r, c, source.At(r, c) * factor);
                }
            }
            return source;
        }
        #endregion

        #region RowMulti
        /// <summary>
        /// 行乘以一个系数
        /// </summary>
        /// <param name="row"></param>
        /// <param name="factor"></param>
        public void RowMulti(int row, double factor)
        {
            for (int c = 0; c < ColumnCount; c++)
            {
                At(row, c, At(row, c) * factor);
            }
        }
        #endregion

        #region ColumnMulti
        /// <summary>
        /// 列乘以一个系数
        /// </summary>
        /// <param name="column"></param>
        /// <param name="factor"></param>
        public void ColumnMulti(int column, double factor)
        {
            for (int row = 0; row < RowCount; row++)
            {
                At(row, column, At(row, column) * factor);
            }
        }
        #endregion
    }

转载请注明出处