using System; namespace Science.Mathematics.LinearAlgebra { /// /// Matrix /// public class Matrix { // Norm, Frobenius Norm protected int sr, sc; public int SizeOfRow { get{return sr;} } public int SizeOfColumn { get{return sc;} } public Matrix() { } protected double[,] el; public Matrix(int numberOfRows, int numberOfColumns) { sr = numberOfRows; sc = numberOfColumns; el = new double[numberOfRows, numberOfColumns]; } public Matrix(double[,] x) { sr = x.GetLength(0); sc = x.GetLength(1); el = new double[sr,sc]; for(int k = 0; k < sr; k++) for(int kk = 0; kk < sc; kk++) el[k,kk] = x[k,kk]; } public Matrix(Vector[] columnVector) { sr = columnVector[0].Size; sc = columnVector.Length; el = new double[sr, sc]; for (int k = 0; k < sr; k++) for (int kk = 0; kk < sc; kk++) el[k, kk] = columnVector[kk][k]; } public double this[int whereRow, int whereColumn] { get{return el[whereRow, whereColumn];} set{el[whereRow, whereColumn]=value;} } public int Rank { get { FactorizationEAeqR obj = new FactorizationEAeqR(this); obj.Compute(); int nofpivot = 0; for (int i = 0; i < this.SizeOfRow; i++) for (int j = 0; j < this.SizeOfColumn; j++) if (obj.ReducedRowEchelonForm[i, j] > 1.0E-12) { nofpivot++; break; } return nofpivot; } } public Matrix Transpose { get { Matrix m = new Matrix(sc,sr); for(int k = 0; k < sr; k++) for(int kk = 0; kk < sc; kk++) m[kk,k] = el[k,kk]; return m; } } public double Determinant { get { if (sr != sc) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("Square matrix should be considered."); } double[,] ell = new double[sr, sc]; for (int j = 0; j < sr; j++) for (int i = 0; i < sc; i++) ell[j, i] = el[j, i]; LUdcmp obj = new LUdcmp(ell); return obj.det(); } } public double Trace { get { if (sr != sc) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("Square matrix should be considered."); } double sum = 0.0; for (int k = 0; k < this.SizeOfRow; k++) sum += this.el[k, k]; return sum; } } public Matrix Inverse { get { if (sr != sc) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("Square matrix should be considered."); } double[,] ell = new double[sr, sc]; for (int j = 0; j < sr; j++) for (int i = 0; i < sc; i++) ell[j, i] = el[j, i]; LUdcmp obj = new LUdcmp(ell); double[,] a = obj.inverse(); Matrix im = new Matrix(sr, sc); for (int j = 0; j < sr; j++) for (int i = 0; i < sc; i++) im[j, i] = a[j, i]; return im; } } public Matrix Similar(Matrix M) { return M.Inverse * this * M; } public Matrix Pseudoinverse { get { FactorizationAeqUSVt obj = new FactorizationAeqUSVt(this); obj.Compute(); Matrix mi = new Matrix(sc, sc); for (int i = 0; i < sc; i++) { if (Math.Abs(obj.Diagonal[i, i]) < 1.0E-13) mi[i, i] = 0.0; else mi[i, i] = 1.0 / obj.Diagonal[i, i]; } Matrix pim = obj.Orthogonal * mi * obj.ColumnOrthogonal.Transpose; return pim; } } private Complex[] e; private VectorComplex[] v; public virtual void Diagonalize() { if (sr != sc) try { throw new Science.Error(); } catch (Science.Error ee) { ee.Write("Square matrix should be considered."); } FactorizationAeqSLSi obj = new FactorizationAeqSLSi(this); obj.Compute(); e = new Complex[sr]; v = new VectorComplex[sr]; for (int i = 0; i < sr; i++) { e[i] = obj.Eigenvalues[i, i]; v[i] = new VectorComplex(sr); for (int j = 0; j < sr; j++) v[i][j] = obj.Eigenvectors[j, i]; } } public Complex[] EigenvalueComplex { get { return e; } } public VectorComplex[] EigenvectorComplex { get { return v; } } public static Matrix operator *(double a, Matrix c) { Matrix b = new Matrix(c.SizeOfRow, c.SizeOfColumn); for (int k = 0; k < c.SizeOfRow; k++) for (int j = 0; j < c.SizeOfColumn; j++) b[k, j] = a * c[k, j]; return b; } public static Matrix operator *(Matrix c, double a) { return a*c; } public static Matrix operator +(Matrix a, Matrix c) { Matrix b = new Matrix(a.SizeOfRow, a.SizeOfColumn); for (int k = 0; k < a.SizeOfRow; k++) for (int j = 0; j < a.SizeOfColumn; j++) b[k, j] = a[k, j] + c[k, j]; return b; } public static Matrix operator -(Matrix a, Matrix c) { Matrix b = new Matrix(a.SizeOfRow, a.SizeOfColumn); for (int k = 0; k < a.SizeOfRow; k++) for (int j = 0; j < a.SizeOfColumn; j++) b[k, j] = a[k, j] - c[k, j]; return b; } public static Matrix operator -(Matrix a) { Matrix b = new Matrix(a.SizeOfRow, a.SizeOfColumn); for (int k = 0; k < a.SizeOfRow; k++) for (int j = 0; j < a.SizeOfColumn; j++) b[k, j] = - a[k, j]; return b; } public static Vector operator *(Matrix a, Vector c) { Vector b = new Vector(a.SizeOfRow); for(int k = 0; k < b.Size; k++) { double sum = 0.0; for(int s = 0; s < c.Size; s++) sum += a[k,s]*c[s]; b[k] = sum; } return b; } public static Vector operator *(Vector c, Matrix a) { Vector b = new Vector(a.SizeOfColumn); for (int k = 0; k < b.Size; k++) { double sum = 0.0; for (int s = 0; s < c.Size; s++) sum += c[s] * a[s, k]; b[k] = sum; } return b; } public static Matrix operator *(Matrix a, Matrix c) { Matrix b = new Matrix(a.SizeOfRow, c.SizeOfColumn); double sum; for (int k = 0; k < b.SizeOfRow; k++) for (int j = 0; j < b.SizeOfColumn; j++) { sum = 0.0; for (int s = 0; s < a.SizeOfColumn; s++) sum += a[k, s] * c[s, j]; b[k, j] = sum; } return b; } public void Chop() { for (int k = 0; k < sr; k++) { for (int kk = 0; kk < sc; kk++) { if (el[k, kk] < 0.00000001 && el[k, kk] > -0.00000001) el[k, kk] = 0.0; } } } public override string ToString() { string res = ""; for (int i = 0; i < this.SizeOfRow; i++) { for (int j = 0; j < this.SizeOfColumn; j++) res += this[i, j].ToString()+" "; res += "\r\n"; } res += "\r\n"; return res; } } }