Run Code
|
API
|
Code Wall
|
Misc
|
Feedback
|
Login
|
Theme
|
Privacy
|
Patreon
Square Real Matrix
//Rextester.Program.Main is the entry point for your code. Don't change it. //Compiler version 4.0.30319.17929 for Microsoft (R) .NET Framework 4.5 using System; using System.Collections.Generic; using System.Linq; using System.Text.RegularExpressions; using System.Text; namespace Rextester { public class HelpFunctions { public static int[,] IntArrayFromDouble(double[][] input) { int rows = input.GetLength(0); int columns = input.GetLength(1); int[,] ret = new int[rows, columns]; for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { ret[i, j] = (int)input[i][j]; } } return ret; } public static double[][] DoubleArrayFromInt(int[,] input) { int rows = input.GetLength(0); int columns = input.GetLength(1); double[][] ret = MatrixCreate(rows, columns); for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { ret[i][j] = (double)input[i, j]; } } return ret; } public static double[,] DoubleFromDoubleArray(double[][] input) { int rows = input.GetLength(0); int columns = input.GetLength(0); double[,] ret = new double[rows, columns]; for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { ret[i, j] = input[i][j]; } } return ret; } public static double[][] DoubleArrayFromDouble(double[,] input) { int rows = input.GetLength(0); int columns = input.GetLength(1); double[][] ret = MatrixCreate(rows, columns); for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { ret[i][j] = input[i, j]; } } return ret; } public static void Test() { Console.WriteLine("\nBegin matrix inverse using Crout LU decomp demo \n"); double[][] m = MatrixCreate(4, 4); m[0][0] = 3.0; m[0][1] = 7.0; m[0][2] = 2.0; m[0][3] = 5.0; m[1][0] = 1.0; m[1][1] = 8.0; m[1][2] = 4.0; m[1][3] = 2.0; m[2][0] = 2.0; m[2][1] = 1.0; m[2][2] = 9.0; m[2][3] = 3.0; m[3][0] = 5.0; m[3][1] = 4.0; m[3][2] = 7.0; m[3][3] = 1.0; Console.WriteLine("Original matrix m is "); Console.WriteLine(MatrixAsString(m)); double d = MatrixDeterminant(m); if (Math.Abs(d) < 1.0e-5) Console.WriteLine("Matrix has no inverse"); double[][] inv = MatrixInverse(m); Console.WriteLine("Inverse matrix inv is "); Console.WriteLine(MatrixAsString(inv)); double[][] prod = MatrixProduct(m, inv); Console.WriteLine("The product of m * inv is "); Console.WriteLine(MatrixAsString(prod)); Console.WriteLine("========== \n"); double[][] lum; int[] perm; int toggle = MatrixDecompose(m, out lum, out perm); Console.WriteLine("The combined lower-upper decomposition of m is"); Console.WriteLine(MatrixAsString(lum)); double[][] lower = ExtractLower(lum); double[][] upper = ExtractUpper(lum); Console.WriteLine("The lower part of LUM is"); Console.WriteLine(MatrixAsString(lower)); Console.WriteLine("The upper part of LUM is"); Console.WriteLine(MatrixAsString(upper)); Console.WriteLine("The perm[] array is"); ShowVector(perm); double[][] lowTimesUp = MatrixProduct(lower, upper); Console.WriteLine("The product of lower * upper is "); Console.WriteLine(MatrixAsString(lowTimesUp)); Console.WriteLine("\nEnd matrix inverse demo \n"); Console.ReadLine(); } public static void ShowVector(int[] vector) { Console.Write(" "); for (int i = 0; i < vector.Length; ++i) Console.Write(vector[i] + " "); Console.WriteLine("\n"); } public static double[][] MatrixInverse(double[][] matrix) { // assumes determinant is not 0 // that is, the matrix does have an inverse int n = matrix.Length; double[][] result = MatrixCreate(n, n); // make a copy of matrix for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) result[i][j] = matrix[i][j]; double[][] lum; // combined lower & upper int[] perm; int toggle; toggle = MatrixDecompose(matrix, out lum, out perm); double[] b = new double[n]; for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) if (i == perm[j]) b[j] = 1.0; else b[j] = 0.0; double[] x = Helper(lum, b); // for (int j = 0; j < n; ++j) result[j][i] = x[j]; } return result; } // MatrixInverse public static int MatrixDecompose(double[][] m, out double[][] lum, out int[] perm) { // Crout's LU decomposition for matrix determinant and inverse // stores combined lower & upper in lum[][] // stores row permuations into perm[] // returns +1 or -1 according to even or odd number of row permutations // lower gets dummy 1.0s on diagonal (0.0s above) // upper gets lum values on diagonal (0.0s below) int toggle = +1; // even (+1) or odd (-1) row permutatuions int n = m.Length; // make a copy of m[][] into result lu[][] lum = MatrixCreate(n, n); for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) lum[i][j] = m[i][j]; // make perm[] perm = new int[n]; for (int i = 0; i < n; ++i) perm[i] = i; for (int j = 0; j < n - 1; ++j) // process by column. note n-1 { double max = Math.Abs(lum[j][j]); int piv = j; for (int i = j + 1; i < n; ++i) // find pivot index { double xij = Math.Abs(lum[i][j]); if (xij > max) { max = xij; piv = i; } } // i if (piv != j) { double[] tmp = lum[piv]; // swap rows j, piv lum[piv] = lum[j]; lum[j] = tmp; int t = perm[piv]; // swap perm elements perm[piv] = perm[j]; perm[j] = t; toggle = -toggle; } double xjj = lum[j][j]; if (xjj != 0.0) { for (int i = j + 1; i < n; ++i) { double xij = lum[i][j] / xjj; lum[i][j] = xij; for (int k = j + 1; k < n; ++k) lum[i][k] -= xij * lum[j][k]; } } } // j return toggle; } // MatrixDecompose public static double[] Helper(double[][] luMatrix, double[] b) // helper { int n = luMatrix.Length; double[] x = new double[n]; b.CopyTo(x, 0); for (int i = 1; i < n; ++i) { double sum = x[i]; for (int j = 0; j < i; ++j) sum -= luMatrix[i][j] * x[j]; x[i] = sum; } x[n - 1] /= luMatrix[n - 1][n - 1]; for (int i = n - 2; i >= 0; --i) { double sum = x[i]; for (int j = i + 1; j < n; ++j) sum -= luMatrix[i][j] * x[j]; x[i] = sum / luMatrix[i][i]; } return x; } // Helper public static double MatrixDeterminant(double[][] matrix) { double[][] lum; int[] perm; int toggle = MatrixDecompose(matrix, out lum, out perm); double result = toggle; for (int i = 0; i < lum.Length; ++i) result *= lum[i][i]; return result; } // ---------------------------------------------------------------- public static double[][] MatrixCreate(int rows, int cols) { double[][] result = new double[rows][]; for (int i = 0; i < rows; ++i) result[i] = new double[cols]; return result; } public static double[][] MatrixProduct(double[][] matrixA, double[][] matrixB) { int aRows = matrixA.Length; int aCols = matrixA[0].Length; int bRows = matrixB.Length; int bCols = matrixB[0].Length; if (aCols != bRows) throw new Exception("Non-conformable matrices"); double[][] result = MatrixCreate(aRows, bCols); for (int i = 0; i < aRows; ++i) // each row of A for (int j = 0; j < bCols; ++j) // each col of B for (int k = 0; k < aCols; ++k) // could use k < bRows result[i][j] += matrixA[i][k] * matrixB[k][j]; return result; } public static string MatrixAsString(double[][] matrix) { string s = ""; for (int i = 0; i < matrix.Length; ++i) { for (int j = 0; j < matrix[i].Length; ++j) s += matrix[i][j].ToString("F3").PadLeft(8) + " "; s += Environment.NewLine; } return s; } public static double[][] ExtractLower(double[][] lum) { // lower part of an LU Doolittle decomposition (dummy 1.0s on diagonal, 0.0s above) int n = lum.Length; double[][] result = MatrixCreate(n, n); for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { if (i == j) result[i][j] = 1.0; else if (i > j) result[i][j] = lum[i][j]; } } return result; } public static double[][] ExtractUpper(double[][] lum) { // upper part of an LU (lu values on diagional and above, 0.0s below) int n = lum.Length; double[][] result = MatrixCreate(n, n); for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { if (i <= j) result[i][j] = lum[i][j]; } } return result; } } public class SquareRealMatrix { private double[,] InternalRep = null; public int Rows = 0; public int Columns = 0; private List<double> Vector = null; private void Zero() { for (int i = 0; i < Rows; i++) { for (int j = 0; j < Columns; j++) { InternalRep[i, j] = 0; } } } private void FromVector() { int cnt = 0; for (int i = 0; i < Rows; i++) { for (int j = 0; j < Columns; j++) { InternalRep[i, j] = Vector[cnt++]; } } } private int SignOfElement(int i, int j) { if ((i + j) % 2 == 0) { return 1; } else { return -1; } } public SquareRealMatrix(double[,] RC) { int rows = RC.GetLength(0); int columns = RC.GetLength(1); if (rows != columns) { throw new Exception("rows and columns must be equal for square matrix"); } this.Rows = rows; this.Columns = columns; InternalRep = RC; } //this method determines the sub matrix corresponding to a given element private double[,] CreateSmallerMatrix(double[,] input, int i, int j) { int order = int.Parse(System.Math.Sqrt(input.Length).ToString()); double[,] output = new double[order - 1, order - 1]; int x = 0, y = 0; for (int m = 0; m < order; m++, x++) { if (m != i) { y = 0; for (int n = 0; n < order; n++) { if (n != j) { output[x, y] = input[m, n]; y++; } } } else { x--; } } return output; } public double Determinant() { return Determinant(this.InternalRep); } private double Determinant(double[,] input) { int order = int.Parse(System.Math.Sqrt(input.Length).ToString()); if (order > 2) { double value = 0; for (int j = 0; j < order; j++) { double[,] Temp = CreateSmallerMatrix(input, 0, j); value = value + input[0, j] * (SignOfElement(0, j) * Determinant(Temp)); } return value; } else if (order == 2) { return ((input[0, 0] * input[1, 1]) - (input[1, 0] * input[0, 1])); } else { return (input[0, 0]); } } public SquareRealMatrix(int rows, int columns) { if (rows != columns) { throw new Exception("rows and columns must be equal for square matrix"); } this.Rows = rows; this.Columns = columns; InternalRep = new double[this.Rows, this.Columns]; Zero(); } public SquareRealMatrix(int rows, int columns, List<double> V) { if (rows != columns) { throw new Exception("rows and columns must be equal for square matrix"); } if (V.Count % rows != 0) { throw new Exception("Vector does not contain even row count"); } Vector = V; this.Rows = rows; this.Columns = columns; InternalRep = new double[this.Rows, this.Columns]; FromVector(); } public SquareRealMatrix Inverse() { double[][] m = HelpFunctions.DoubleArrayFromDouble(this.InternalRep); double[][] inv = HelpFunctions.MatrixInverse(m); double[,] mi = HelpFunctions.DoubleFromDoubleArray(inv); return new SquareRealMatrix(mi); } public static SquareRealMatrix operator +(SquareRealMatrix a, SquareRealMatrix b) { SquareRealMatrix ret = new SquareRealMatrix(a.Rows, a.Columns); for (int i = 0; i < ret.Rows; i++) { for (int j = 0; j < ret.Columns; j++) { ret.InternalRep[i, j] = a.InternalRep[i, j] + b.InternalRep[i, j]; } } return ret; } public static SquareRealMatrix operator *(SquareRealMatrix a, SquareRealMatrix b) { SquareRealMatrix ret = new SquareRealMatrix(a.Rows, a.Columns); for (int i = 0; i < ret.Rows; i++) { for (int j = 0; j < ret.Columns; j++) { for (int k = 0; k < ret.Columns; k++) { ret.InternalRep[i, j] += a.InternalRep[i, k] * b.InternalRep[k, j]; double sens = ret.InternalRep[i, j]; if (Math.Abs(sens) < 1.0e-8d) //f**ked up value set to zero { ret.InternalRep[i, j] = 0; } } } } return ret; } public double this[int r, int c] { get { if (!(r < Rows && c < Columns)) { throw new Exception("rows and columns out of range of square matrix"); } return (InternalRep[r, c]); } set { InternalRep[r, c] = value; } } public override string ToString() { StringBuilder sb = new StringBuilder(); for (int i = 0; i < Rows; i++) { for (int j = 0; j < Columns; j++) { sb.AppendFormat("{0}\t", InternalRep[i, j]); } sb.Append("\r\n"); } return sb.ToString(); } } public class Program { public static void Main(string[] args) { List<double> v = new List<double> { 1, 2, 3, 4 }; List<double> v2 = new List<double> { 7, 8, 9, 10 }; SquareRealMatrix m = new SquareRealMatrix(2, 2, v); SquareRealMatrix m2 = new SquareRealMatrix(2, 2, v2); SquareRealMatrix mm = m * m2; Console.WriteLine(mm.ToString()); SquareRealMatrix s = new SquareRealMatrix(2, 2, v); SquareRealMatrix i = s.Inverse(); Console.WriteLine(s.ToString()); Console.WriteLine(i.ToString()); SquareRealMatrix id = s * i; Console.WriteLine(id.ToString()); } } }
run
|
edit
|
history
|
help
0
araray
Practice run
SquareRealMatrixWithMathJax
Stacks: Balanced Brackets
Order of Ops 4.2
Check if string is palindrome
Lambda Expressions Are Cool 2
Palindrome Number
Date Diff
test code