測量平差以及工科中常常用到矩陣的相關運算，因此自己寫了一個，同時考慮到了類庫的可用性，這次又重載了比較勻運算符，修正了一些問題
using System; using System.Collections.Generic; namespace CMath { [Serializable] public class Matrix { public double[] element; private int rows = 0; private int cols = 0; /// <summary> /// 獲取矩陣行數 /// </summary> public int Rows { get { return rows; } } /// <summary> /// 獲取矩陣列數 /// </summary> public int Cols { get { return cols; } } /// <summary> /// 獲取或設置第i行第j列的元素值 /// </summary> /// <param name="i">第i行</param> /// <param name="j">第j列</param> /// <returns>返回第i行第j列的元素值</returns> public double this[int i, int j] { get { if (i < Rows && j < Cols) { return element[i * cols + j]; } else { throw new Exception("索引越界"); } } set { element[i * cols + j] = value; } } /// <summary> /// 用二維數組初始化Matrix /// </summary> /// <param name="m">二維數組</param> public Matrix(double[][] m) { this.rows = m.GetLength(0); this.cols = m.GetLength(1); int count = 0; this.element=new double[Rows*Cols]; for (int i = 0; i < rows; i++) { for (int j = 0; j < cols; j++) { element[count++] = m[i][j]; } } } public Matrix(double[,] m) { this.rows = m.GetLength(0); this.cols = m.GetLength(1); this.element = new double[this.rows * this.cols]; int count = 0; for (int i = 0; i < rows; i++) { for (int j = 0; j < cols; j++) { element[count++] = m[i, j]; } } } public Matrix(List<List<double>> m) { this.rows = m.Count; this.cols = m[0].Count; this.element = new double[Rows * Cols]; for (int i = 0; i < rows; i++) { for (int j = 0; j < cols; j++) { this[i, j] = m[i][j]; } } } #region 矩陣數學運算 public static Matrix MAbs(Matrix a) { Matrix _thisCopy = a.DeepCopy(); for (int i = 0; i < a.Rows; i++) { for (int j = 0; j < a.Cols; j++) { _thisCopy[i, j] = Math.Abs(a[i, j]); } } return _thisCopy; } /// <summary> /// 矩陣相加 /// </summary> /// <param name="a">第一個矩陣,和b矩陣必須同等大小</param> /// <param name="b">第二個矩陣</param> /// <returns>返回矩陣相加后的結果</returns> public static Matrix operator +(Matrix a, Matrix b) { if (a.cols == b.cols && a.rows == b.rows) { double[,] res = new double[a.rows, a.cols]; for (int i = 0; i < a.Rows; i++) { for (int j = 0; j < a.Cols; j++) { res[i, j] = a[i, j] + b[i, j]; } } return new Matrix(res); } else { throw new Exception("兩個矩陣行列不相等"); } } /// <summary> /// 矩陣相減 /// </summary> /// <param name="a">第一個矩陣,和b矩陣必須同等大小</param> /// <param name="b">第二個矩陣</param> /// <returns>返回矩陣相減后的結果</returns> public static Matrix operator -(Matrix a, Matrix b) { if (a.cols == b.cols && a.rows == b.rows) { double[,] res = new double[a.rows, a.cols]; for (int i = 0; i < a.Rows; i++) { for (int j = 0; j < a.Cols; j++) { res[i, j] = a[i, j] - b[i, j]; } } return new Matrix(res); } else { throw new Exception("兩個矩陣行列不相等"); } } /// <summary> /// 對矩陣每個元素取相反數 /// </summary> /// <param name="a">二維矩陣</param> /// <returns>得到矩陣的相反數</returns> public static Matrix operator -(Matrix a) { Matrix res = a; for (int i = 0; i < a.rows; i++) { for (int j = 0; j < a.cols; j++) { res.element[i * a.cols + j] = -res.element[i * a.cols + j]; } } return res; } /// <summary> /// 矩陣相乘 /// </summary> /// <param name="a">第一個矩陣</param> /// <param name="b">第二個矩陣,這個矩陣的行要與第一個矩陣的列相等</param> /// <returns>返回相乘后的一個新的矩陣</returns> public static Matrix operator *(Matrix a, Matrix b) { if (a.cols == b.rows) { double[,] res = new double[a.rows, b.cols]; for (int i = 0; i < a.rows; i++) { for (int j = 0; j < b.cols; j++) { for (int k = 0; k < a.cols; k++) { res[i, j] += a[i, k] * b[k, j]; } } } return new Matrix(res); } else { throw new Exception("兩個矩陣行和列不等"); } } /// <summary> /// 矩陣與數相乘 /// </summary> /// <param name="a">第一個矩陣</param> /// <param name="num">一個實數</param> /// <returns>返回相乘后的新的矩陣</returns> public static Matrix operator *(Matrix a, double num) { Matrix res = a; for (int i = 0; i < a.rows; i++) { for (int j = 0; j < a.cols; j++) { res.element[i * a.cols + j] *= num; } } return res; } /// <summary> /// 矩陣轉置 /// </summary> /// <returns>返回當前矩陣轉置后的新矩陣</returns> public Matrix Transpose() { double[,] res = new double[cols, rows]; { for (int i = 0; i < cols; i++) { for (int j = 0; j < rows; j++) { res[i, j] = this[j, i]; } } } return new Matrix(res); } /// <summary> /// 矩陣求逆 /// </summary> /// <returns>返回求逆后的新的矩陣</returns> public Matrix Inverse() { //最后原始矩陣並不變，所以需要深拷貝一份 Matrix _thisCopy = this.DeepCopy(); if (cols == rows && this.Determinant() != 0) { //初始化一個同等大小的單位陣 Matrix res = _thisCopy.EMatrix(); for (int i = 0; i < rows; i++) { //首先找到第i列的絕對值最大的數，並將該行和第i行互換 int rowMax = i; double max = Math.Abs(_thisCopy[i, i]); for (int j = i; j < rows; j++) { if (Math.Abs(_thisCopy[j, i]) > max) { rowMax = j; max = Math.Abs(_thisCopy[j, i]); } } //將第i行和找到最大數那一行rowMax交換 if (rowMax != i) { _thisCopy.Exchange(i, rowMax); res.Exchange(i, rowMax); } //將第i行做初等行變換，將第一個非0元素化為1 double r = 1.0 / _thisCopy[i, i]; _thisCopy.Exchange(i, -1, r); res.Exchange(i, -1, r); //消元 for (int j = 0; j < rows; j++) { //到本行后跳過 if (j == i) continue; else { r = -_thisCopy[j, i]; _thisCopy.Exchange(i, j, r); res.Exchange(i, j, r); } } } return res; } else { throw new Exception("矩陣不是方陣無法求逆"); } } #region 重載比較運算符 public static bool operator <(Matrix a, Matrix b) { bool issmall = true; for (int i = 0; i < a.Rows; i++) { for (int j = 0; j < a.Cols; j++) { if (a[i, j] >= b[i, j]) issmall = false; } } return issmall; } public static bool operator >(Matrix a, Matrix b) { bool issmall = true; for (int i = 0; i < a.Rows; i++) { for (int j = 0; j < a.Cols; j++) { if (a[i, j] <= b[i, j]) issmall = false; } } return issmall; } public static bool operator <=(Matrix a, Matrix b) { bool issmall = true; for (int i = 0; i < a.Rows; i++) { for (int j = 0; j < a.Cols; j++) { if (a[i, j] > b[i, j]) issmall = false; } } return issmall; } public static bool operator >=(Matrix a, Matrix b) { bool issmall = true; for (int i = 0; i < a.Rows; i++) { for (int j = 0; j < a.Cols; j++) { if (a[i, j] < b[i, j]) issmall = false; } } return issmall; } public static bool operator !=(Matrix a, Matrix b) { bool issmall = true; issmall = ReferenceEquals(a, b); if (issmall) return issmall; for (int i = 0; i < a.Rows; i++) { for (int j = 0; j < a.Cols; j++) { if (a[i, j] == b[i, j]) issmall = false; } } return issmall; } public static bool operator ==(Matrix a, Matrix b) { bool issmall = true; issmall = ReferenceEquals(a, b); if (issmall) return issmall; for (int i = 0; i < a.Rows; i++) { for (int j = 0; j < a.Cols; j++) { if (a[i, j] != b[i, j]) issmall = false; } } return issmall; } public override bool Equals(object obj) { Matrix b = obj as Matrix; return this == b; } public override int GetHashCode() { return base.GetHashCode(); } #endregion public double Determinant() { if (cols == rows) { Matrix _thisCopy = this.DeepCopy(); //行列式每次交換行，都需要乘以-1 double res = 1; for (int i = 0; i < rows; i++) { //首先找到第i列的絕對值最大的數 int rowMax = i; double max = Math.Abs(_thisCopy[i, i]); for (int j = i; j < rows; j++) { if (Math.Abs(_thisCopy[j, i]) > max) { rowMax = j; max = Math.Abs(_thisCopy[j, i]); } } //將第i行和找到最大數那一行rowMax交換,同時將單位陣做相同初等變換 if (rowMax != i) { _thisCopy.Exchange(i, rowMax); res *= -1; } //消元 for (int j = i + 1; j < rows; j++) { double r = -_thisCopy[j, i] / _thisCopy[i, i]; _thisCopy.Exchange(i, j, r); } } //計算對角線乘積 for (int i = 0; i < rows; i++) { res *= _thisCopy[i, i]; } return res; } else { throw new Exception("不是行列式"); } } #endregion #region 初等變換 /// <summary> /// 初等變換：交換第r1和第r2行 /// </summary> /// <param name="r1">第r1行</param> /// <param name="r2">第r2行</param> /// <returns>返回交換兩行后的新的矩陣</returns> public Matrix Exchange(int r1, int r2) { if (Math.Min(r2, r1) >= 0 && Math.Max(r1, r2) < rows) { for (int j = 0; j < cols; j++) { double temp = this[r1, j]; this[r1, j] = this[r2, j]; this[r2, j] = temp; } return this; } else { throw new Exception("超出索引"); } } /// <summary> /// 初等變換：將r1行乘以某個數加到r2行 /// </summary> /// <param name="r1">第r1行乘以num</param> /// <param name="r2">加到第r2行，若第r2行為負，則直接將r1乘以num並返回</param> /// <param name="num">某行放大的倍數</param> /// <returns></returns> public Matrix Exchange(int r1, int r2, double num) { if (Math.Min(r2, r1) >= 0 && Math.Max(r1, r2) < rows) { for (int j = 0; j < cols; j++) { this[r2, j] += this[r1, j] * num; } return this; } else if (r2 < 0) { for (int j = 0; j < cols; j++) { this[r1, j] *= num; } return this; } else { throw new Exception("超出索引"); } } /// <summary> /// 得到一個同等大小的單位矩陣 /// </summary> /// <returns>返回一個同等大小的單位矩陣</returns> public Matrix EMatrix() { if (rows == cols) { double[,] res = new double[rows, cols]; for (int i = 0; i < rows; i++) { for (int j = 0; j < cols; j++) { if (i == j) res[i, j] = 1; else res[i, j] = 0; } } return new Matrix(res); } else throw new Exception("不是方陣，無法得到單位矩陣"); } #endregion /// <summary> /// 深拷貝，僅僅將值拷貝給一個新的對象 /// </summary> /// <returns>返回深拷貝后的新對象</returns> public Matrix DeepCopy() { double[,] ele = new double[rows, cols]; for (int i = 0; i < rows; i++) { for (int j = 0; j < cols; j++) { ele[i, j] = this[i, j]; } } return new Matrix(ele); } public override string ToString() { string str = ""; for (int i = 0; i < Rows; i++) { for (int j = 0; j < Cols; j++) { str += this[i, j].ToString(); if (j != Cols - 1) str += " "; else if (i != Rows - 1) str += Environment.NewLine; } } return str; } } }
矩陣的求秩過幾天補上。