在一些算法中需要用到矩阵,自然就需要用到矩阵的一些操作,比如行变换、列变换、最简式、求矩阵的秩等,下面是实现的代码
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
}
转载请注明出处