using System;
using System.Collections.Generic;
using System.Text;
namespace LinearAlgebra
{
public class Matrix
{
private double[,] _data;
public Matrix(int size)
{
this._data = new double[size, size];
}
public Matrix(int rows, int cols)
{
this._data = new double[rows, cols];
}
public Matrix(double[,] data)
{
this._data = data;
}
public static Matrix operator +(Matrix matrix1, Matrix matrix2)
{
int matrix1_Rows = matrix1.Data.GetLength(0);
int matrix1_Columns = matrix1.Data.GetLength(1);
int matrix2_Rows = matrix2.Data.GetLength(0);
int matrix2_Columns = matrix2.Data.GetLength(1);
if ((matrix1_Rows != matrix2_Rows) || (matrix1_Columns != matrix2_Columns))
{
throw new Exception("Matrix Dimensions Don't Agree!");
}
double[,] result = new double[matrix1_Rows, matrix1_Columns];
for (int i = 0; i < matrix1_Rows; i++)
{
for (int j = 0; j < matrix1_Columns; j++)
{
result[i, j] = matrix1.Data[i, j] + matrix2.Data[i, j];
}
}
return new Matrix(result);
}
public static Matrix operator -(Matrix matrix1, Matrix matrix2)
{
int matrix1_Rows = matrix1.Data.GetLength(0);
int matrix1_Columns = matrix1.Data.GetLength(1);
int matrix2_Rows = matrix2.Data.GetLength(0);
int matrix2_Columns = matrix2.Data.GetLength(1);
if ((matrix1_Rows != matrix2_Rows) || (matrix1_Columns != matrix2_Columns))
{
throw new Exception("Matrix Dimensions Don't Agree!");
}
double[,] result = new double[matrix1_Rows, matrix1_Columns];
for (int i = 0; i < matrix1_Rows; i++)
{
for (int j = 0; j < matrix1_Columns; j++)
{
result[i, j] = matrix1.Data[i, j] - matrix2.Data[i, j];
}
}
return new Matrix(result);
}
public static Matrix operator *(Matrix matrix1, Matrix matrix2)
{
int matrix1_Rows = matrix1.Data.GetLength(0);
int matrix1_Columns = matrix1.Data.GetLength(1);
int matrix2_Rows = matrix2.Data.GetLength(0);
int matrix2_Columns = matrix2.Data.GetLength(1);
if (matrix1_Columns != matrix2_Rows)
{
throw new Exception("Matrix Dimensions Don't Agree!");
}
double[,] result = new double[matrix1_Rows, matrix2_Columns];
for (int i = 0; i < matrix1_Rows; i++)
{
for (int j = 0; j < matrix2_Columns; j++)
{
for (int k = 0; k < matrix2_Rows; k++)
{
result[i, j] += matrix1.Data[i, k] * matrix2.Data[k, j];
}
}
}
return new Matrix(result);
}
public static Matrix operator /(double i, Matrix matrix)
{
return new Matrix(ScaleBy(i, INV(matrix.Data)));
}
public static bool operator ==(Matrix matrix1, Matrix matrix2)
{
bool result = true;
int matrix1_Rows = matrix1.Data.GetLength(0);
int matrix1_Columns = matrix1.Data.GetLength(1);
int matrix2_Rows = matrix2.Data.GetLength(0);
int matrix2_Columns = matrix2.Data.GetLength(1);
if ((matrix1_Rows != matrix2_Rows) || (matrix1_Columns != matrix2_Columns))
{
result = false;
}
else
{
for (int i = 0; i < matrix1_Rows; i++)
{
for (int j = 0; j < matrix1_Columns; j++)
{
if (matrix1.Data[i, j] != matrix2.Data[i, j]) result = false;
}
}
}
return result;
}
public static bool operator !=(Matrix matrix1, Matrix matrix2)
{
return !(matrix1 == matrix2);
}
public override int GetHashCode()
{
return base.GetHashCode();
}
public override bool Equals(object obj)
{
if (!(obj is Matrix)) return false;
return this == (Matrix)obj;
}
public void Display()
{
int rows = this.Data.GetLength(0);
int columns = this.Data.GetLength(1);
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < columns; j++)
{
Console.Write(this.Data[i, j].ToString("N2") + " ");
}
Console.WriteLine();
}
Console.WriteLine();
}
public void Display(string format)
{
int rows = this.Data.GetLength(0);
int columns = this.Data.GetLength(1);
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < columns; j++)
{
Console.Write(this.Data[i, j].ToString(format) + " ");
}
Console.WriteLine();
}
Console.WriteLine();
}
public void Display(System.Windows.Forms.Control container)
{
int rows = this.Data.GetLength(0);
int columns = this.Data.GetLength(1);
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < columns; j++)
{
if (container is System.Windows.Forms.RichTextBox)
{
((System.Windows.Forms.RichTextBox)container).AppendText(this.Data[i, j].ToString("N2") + " ");
}
else if (container is System.Windows.Forms.TextBox)
{
((System.Windows.Forms.TextBox)container).Text = this.Data[i, j].ToString("N2") + " ";
}
else if (container is System.Windows.Forms.ListBox)
{
((System.Windows.Forms.ListBox)container).Items.Add(this.Data[i, j].ToString("N2"));
}
}
}
}
public Matrix Inverse()
{
if ((this.IsSquare) && (!this.IsSingular))
{
return new Matrix(INV(this.Data));
}
else
{
throw new System.Exception(@"Cannot find inverse for non square /singular matrix");
}
}
public Matrix Transpose()
{
double[,] D = Transpose(this.Data);
return new Matrix(D);
}
public static Matrix Zeros(int size)
{
double[,] D = new double[size, size];
return new Matrix(D);
}
public static Matrix Zeros(int rows, int cols)
{
double[,] D = new double[rows, cols];
return new Matrix(D);
}
public Matrix LinearSolve(Matrix COF, Matrix CON)
{
return COF.Inverse() * CON;
}
public double Det()
{
if (this.IsSquare)
{
return Determinent(this.Data);
}
else
{
throw new System.Exception("Cannot Determine the DET for a non square matrix");
}
}
public bool IsSquare
{
get
{
return (this.Data.GetLength(0) == this.Data.GetLength(1));
}
}
public bool IsSingular
{
get
{
return ((int)this.Det() == 0);
}
}
public int Rows { get { return this.Data.GetLength(0); } }
public int Columns { get { return this.Data.GetLength(1); } }
private static double[,] INV(double[,] srcMatrix)
{
int rows = srcMatrix.GetLength(0);
int columns = srcMatrix.GetLength(1);
try
{
if (rows != columns)
throw new Exception();
}
catch
{
Console.WriteLine("Cannot find inverse for an non-square matrix");
}
int q;
double[,] desMatrix = new double[rows, columns];
double[,] unitMatrix = UnitMatrix(rows);
for (int p = 0; p < rows; p++)
{
for (q = 0; q < columns; q++)
{
desMatrix[p, q] = srcMatrix[p, q];
}
}
int i = 0;
double det = 1;
if (srcMatrix[0, 0] == 0)
{
i = 1;
while (i < rows)
{
if (srcMatrix[i, 0] != 0)
{
Matrix.InterRow(srcMatrix, 0, i);
Matrix.InterRow(unitMatrix, 0, i);
det *= -1;
break;
}
i++;
}
}
det *= srcMatrix[0, 0];
Matrix.RowDiv(unitMatrix, 0, srcMatrix[0, 0]);
Matrix.RowDiv(srcMatrix, 0, srcMatrix[0, 0]);
for (int p = 1; p < rows; p++)
{
q = 0;
while (q < p)
{
Matrix.RowSub(unitMatrix, p, q, srcMatrix[p, q]);
Matrix.RowSub(srcMatrix, p, q, srcMatrix[p, q]);
q++;
}
if (srcMatrix[p, p] != 0)
{
det *= srcMatrix[p, p];
Matrix.RowDiv(unitMatrix, p, srcMatrix[p, p]);
Matrix.RowDiv(srcMatrix, p, srcMatrix[p, p]);
}
if (srcMatrix[p, p] == 0)
{
for (int j = p + 1; j < columns; j++)
{
if (srcMatrix[p, j] != 0)
{
for (int k = 0; k < rows; k++)
{
for (q = 0; q < columns; q++)
{
srcMatrix[k, q] = desMatrix[k, q];
}
}
return Inverse(desMatrix);
}
}
}
}
for (int p = rows - 1; p > 0; p--)
{
for (q = p - 1; q >= 0; q--)
{
Matrix.RowSub(unitMatrix, q, p, srcMatrix[q, p]);
Matrix.RowSub(srcMatrix, q, p, srcMatrix[q, p]);
}
}
for (int p = 0; p < rows; p++)
{
for (q = 0; q < columns; q++)
{
srcMatrix[p, q] = desMatrix[p, q];
}
}
return unitMatrix;
}
/// <summary>
/// Inverse the Matrix
/// </summary>
/// <param name="srcMatrix">The Matrix to be Inversed</param>
/// <returns>The Inversed Matrix</returns>
private static double[,] Inverse(double[,] srcMatrix)
{
int rows = srcMatrix.GetLength(0);
int columns = srcMatrix.GetLength(1);
double[,] desMatrix = new double[rows, columns];
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < columns; j++)
{
desMatrix[i, j] = 0;
}
}
try
{
if (rows != columns) throw new Exception();
}
catch
{
Console.WriteLine("Cannot Find Inverse for an non-square Matrix");
}
double determine = Determinent(srcMatrix);
try
{
if (determine == 0)
{
throw new Exception("Cannot Perform Inversion. Matrix Singular");
}
}
catch (Exception ex)
{
Console.WriteLine(ex.Message);
}
for (int p = 0; p < rows; p++)
{
for (int q = 0; q < columns; q++)
{
double[,] tmp = FilterMatrix(srcMatrix, p, q);
double determineTMP = Determinent(tmp);
desMatrix[p, q] = Math.Pow(-1, p + q + 2) * determineTMP / determine;
}
}
desMatrix = Transpose(desMatrix);
return desMatrix;
}
/// <summary>
/// Calculate the Determinent
/// </summary>
/// <param name="srcMatrix">The Matrix used to calc</param>
/// <returns>The Determinent</returns>
private static double Determinent(double[,] srcMatrix)
{
int q = 0;
int rows = srcMatrix.GetLength(0);
int columns = srcMatrix.GetLength(1);
double[,] desMatrix = new double[rows, columns];
for (int p = 0; p < rows; p++)
{
for (q = 0; q < columns; q++)
{
desMatrix[p, q] = srcMatrix[p, q];
}
}
int i = 0;
double det = 1;
try
{
if (rows != columns)
{
throw new Exception("Error: Matrix not Square");
}
}
catch(Exception ex)
{
Console.WriteLine(ex.Message);
}
try
{
if (rows == 0)
{
throw new Exception("Dimension of the Matrix 0X0");
}
}
catch (Exception ex)
{
Console.WriteLine(ex.Message);
}
if (rows == 2)
{
return ((srcMatrix[0, 0] * srcMatrix[1, 1]) - (srcMatrix[0, 1] * srcMatrix[1, 0]));
}
if (srcMatrix[0, 0] == 0)
{
i = 1;
while (i < rows)
{
if (srcMatrix[i, 0] != 0)
{
Matrix.InterRow(srcMatrix, 0, i);
det *= -1;
break;
}
i++;
}
}
if (srcMatrix[0, 0] == 0) return 0;
det *= srcMatrix[0, 0];
Matrix.RowDiv(srcMatrix, 0, srcMatrix[0, 0]);
for (int p = 1; p < rows; p++)
{
q = 0;
while (q < p)
{
Matrix.RowSub(srcMatrix, p, q, srcMatrix[p, q]);
q++;
}
if (srcMatrix[p, p] != 0)
{
det *= srcMatrix[p, p];
Matrix.RowDiv(srcMatrix, p, srcMatrix[p, p]);
}
if (srcMatrix[p, p] == 0)
{
for (int j = p + 1; j < columns; j++)
{
if (srcMatrix[p, j] != 0)
{
Matrix.ColumnSub(srcMatrix, p, j, -1);
det *= srcMatrix[p, p];
Matrix.RowDiv(srcMatrix, p, srcMatrix[p, p]);
break;
}
}
}
if (srcMatrix[p, p] == 0) return 0;
}
for (int p = 0; p < rows; p++)
{
for (q = 0; q < columns; q++)
{
srcMatrix[p, q] = desMatrix[p, q];
}
}
return det;
}
/// <summary>
/// Scale the Matrix by a specified ratio
/// </summary>
/// <param name="scalar">The Ratio for the Scale</param>
/// <param name="srcMatrix">The Matrix which will be scaled</param>
/// <returns>A new Matrix which is scaled with a specified ratio</returns>
private static double[,] ScaleBy(double scalar, double[,] srcMatrix)
{
int rows = srcMatrix.GetLength(0);
int columns = srcMatrix.GetLength(1);
double[,] desMatrix = new double[rows, columns];
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < columns; j++)
{
desMatrix[i, j] = scalar * srcMatrix[i, j];
}
}
return desMatrix;
}
/// <summary>
/// To get a unit matrix
/// </summary>
/// <param name="dimension">Dimension</param>
/// <returns>The unit double array</returns>
private static double[,] UnitMatrix(int dimension)
{
double[,] a = new double[dimension, dimension];
for (int i = 0; i < dimension; i++)
{
for (int j = 0; j < dimension; j++)
{
if (i == j) a[i, j] = 1;
else a[i, j] = 0;
}
}
return a;
}
/// <summary>
/// Scale a specified Row
/// </summary>
/// <param name="srcMatrix">The Matrix to be scaled</param>
/// <param name="row">The Specified Row</param>
/// <param name="scaleRatio">The Scale Ratio</param>
private static void RowDiv(double[,] srcMatrix, int row, double scaleRatio)
{
int columns = srcMatrix.GetLength(1);
for (int i = 0; i < columns; i++)
{
srcMatrix[row, i] /= scaleRatio;
}
}
/// <summary>
/// Substract a specified Row from another Row
/// </summary>
/// <param name="srcMatrix">The Matrix to be substract</param>
/// <param name="row1">The Row Index to be substracted</param>
/// <param name="row2">The Row Index to substract</param>
/// <param name="scaleRatio">Scale Ratio</param>
private static void RowSub(double[,] srcMatrix, int row1, int row2, double scaleRatio)
{
int columns = srcMatrix.GetLength(1);
for (int q = 0; q < columns; q++)
{
srcMatrix[row1, q] -= scaleRatio * srcMatrix[row2, q];
}
}
/// <summary>
/// Substract a specified Column from another Column
/// </summary>
/// <param name="srcMatrix">The Matrix to be substracted</param>
/// <param name="column1">The Column Index to be Substracted</param>
/// <param name="column2">The Column Index to Substract</param>
/// <param name="scaleRatio">Scale Ratio</param>
private static void ColumnSub(double[,] srcMatrix, int column1, int column2, double scaleRatio)
{
int rows = srcMatrix.GetLength(0);
int columns = srcMatrix.GetLength(1);
for (int i = 0; i < rows; i++)
{
srcMatrix[i, column1] -= scaleRatio * srcMatrix[i, column2];
}
}
/// <summary>
/// Exchange a specified row's value
/// </summary>
/// <param name="srcMatrix">The Matrix to be exchanged</param>
/// <param name="row1">Row index</param>
/// <param name="row2">Row index</param>
/// <returns>Exchanged Matrix</returns>
private static double[,] InterRow(double[,] srcMatrix, int row1, int row2)
{
int rows = srcMatrix.GetLength(0);
int columns = srcMatrix.GetLength(1);
double tmp = 0;
for (int k = 0; k < columns; k++)
{
tmp = srcMatrix[row1, k];
srcMatrix[row1, k] = srcMatrix[row2, k];
srcMatrix[row2, k] = tmp;
}
return srcMatrix;
}
/// <summary>
/// To Filter the Matrix
/// </summary>
/// <param name="srcMatrix">The Matrix to be Filtered</param>
/// <param name="row">A specified Row</param>
/// <param name="column">A specified Column</param>
/// <returns>The Filtered Matrix</returns>
private static double[,] FilterMatrix(double[,] srcMatrix, int row, int column)
{
int rows = srcMatrix.GetLength(0);
double[,] desMatrix = new double[rows - 1, rows - 1];
int i = 0;
for (int p = 0; p < rows; p++)
{
int j = 0;
if (p != row)
{
for (int q = 0; q < rows; q++)
{
if (q != column)
{
desMatrix[i, j] = srcMatrix[p, q];
j++;
}
}
i++;
}
}
return desMatrix;
}
/// <summary>
/// Transpose Matrix
/// </summary>
/// <param name="srcMatrix">The Matrix to be Transposed</param>
/// <returns>The Transposed Matrix</returns>
private static double[,] Transpose(double[,] srcMatrix)
{
int rows = srcMatrix.GetLength(0);
int columns = srcMatrix.GetLength(1);
double[,] desMatrix = new double[rows, columns];
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < columns; j++)
{
desMatrix[j, i] = srcMatrix[i, j];
}
}
return desMatrix;
}
/// <summary>
/// Gets or Sets the specified value of the Matrix
/// </summary>
/// <param name="i">The row position in the Matrix</param>
/// <param name="j">The Column position in the Matrix</param>
/// <returns>Double, The specified value of the Matrix</returns>
public double this[int i, int j]
{
get { return this._data[i, j]; }
set { this._data[i, j] = value; }
}
/// <summary>
/// Gets the Matrix
/// </summary>
public double[,] Data
{
get { return _data; }
}
}
}
调用:
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
double[,] DATA = { { 2, 1, 3, 5, 9 }, { 2, 2, 1, 10, 1 }, { 7, 3, 5, 2, 1 }, { 3, 7, 7, 7, 9 }, { 12, 4, 1, 23, 5 } };
Matrix COF = new Matrix(DATA);
Matrix CON = new Matrix(5, 1);
CON[0, 0] = 198;
CON[1, 0] = 166;
CON[2, 0] = 115;
CON[3, 0] = 286;
CON[4, 0] = 429;
richTextBox1.AppendText("Determinent of COF matrix = " +
COF.Det().ToString() + " ");
Matrix Sol = COF.Inverse() * CON;
Sol.Display(this.richTextBox1);
}
}
DATA是五元一次方程的系数, Matrix CON是五元一次方程的结果, 即
2X+Y+3Z+5P+9Q=198
2X+2Y+Z+10P+Q=166
7X+3Y+5Z+2P+Q=115
3X+7Y+7Z+7P+9Q=286
12X+4Y+Z+23P+5Q=429
然后, 程序解方程得出的结果是:
X = 4
Y = 5
Z = 7
P = 13
Q = 11