Adding Crout decomposition with partial pivoting to the matrix class

by Trent Guidry30. June 2009 13:35
This is the third part of a series about creating a matrix class in C#.  In this post I will add the solution of matrix equations using LU decomposition with partial pivoting and matrix inversion.  Matrix creation, cloning, transposing, creating an identity matrix, and converting a matrix to and from an array were previously covered in the first part of this series.  Matrix addition, subtraction, and multiplication were covered in the second part of this series.
More information on LU decomposition can be found on Wikipedia LU Decomposition and at Professor Autar K Kaw’s page LU DECOMPOSITION METHOD.
LU decomposition is a method for solving matrix equations.  It is more computationally efficient than naïve Gaussian Elimination in some situations, such a matrix inversion.
A system of linear equations can be written:
Ax = b
Where A is the NxM coefficient matrix, x is the Nx1 column vector of unknowns, and b is the Nx1 column solution vector.
For A, I will assume that there exists two matrices L and U such that LU=A.  Further, for L, every element above the matrix diagonal is 0 and for U, every element below the matrix diagonal is 0.
This can be substituted into the above equation to give:
LUx = b
Multiplying both sides by the inverse of L gives:
Ux = L-1b
d is defined to equal the left and right hand side of the above equation:
d = L-1b = Ux
Multiplying the above equation by L gives:
Ld = b
Using these equations, a matrix equation can be solved by:
1)      Decomposing A into L and U.
2)      Solving for d using Ld = b.
3)      Solving for x using Ux = d.
There is an additional complication that arises when doing this, however, and that is occasional division by 0 errors. In order to try and avoid these division by 0 issues and minimize floating point errors, matrix pivoting can be used.
To do this, I will introduce another matrix called P which is defined by the following equation:
LU = PA
P consists of a matrix with a single 1 in an element of each row and 0s for all of the other elements in that row.  For each column, there can only be one row that has an element value of 1 in that column.  An identity matrix is an example of this type of matrix, as is any matrix produced by interchanging rows of an identity matrix.
Multiplying the equation:
Ax = b
By P gives:
PAx = Pb
Substitution yields:
LUx = Pb
Defining d as above:
Ux = d
Followed by substitution gives:
Ld = Pb
Using these equations, I now  have the following algorithm for solving a matrix equation.
1)      Start with the identity matrix for P.  Solve for L and U, adjusting P as needed.
2)      Solve for d using Ld = Pb
3)      Solve for x using Ux = d
For the algorithm, rather than using a full matrix for P, I will represent it as an array of integers, since that is all that is actually needed.  Each value of the P array represents the column for that row that has a value of 1.  The other column values for that row are zero.
Now to starting coding this up.  I will add a method called SolveFor to the matrix class.  The stub of this function is shown below.
public Matrix SolveFor(Matrix rightMatrix)
{
}
This is the function that will solve the matrix equation.  For the equation,
Ax = b
The unknown matrix x can be solved by calling A.SolveFor(b).  In this case, the function is called on the coefficient matrix and the matrix of solutions, b is passed in.  The return value of the function is the matrix of unknowns.
The first thing I need to do to solve this equation is decompose A into L and U matrices.  For this I will add a private helper function called LUDecompose.   This function returns an instance of a class, LUDecompositionResults, that contains the results of the LU decomposition.  The code for this class is shown below.
public class LUDecompositionResults{private Matrix _lMatrix;private Matrix _uMatrix;private int[] _pivotArray;public LUDecompositionResults(){}public LUDecompositionResults(Matrix matL, Matrix matU, int[] nPivotArray){_lMatrix = matL;_uMatrix = matU;_pivotArray = nPivotArray;}public Matrix L{get { return _lMatrix; }set { _lMatrix = value; }}public Matrix U{get { return _uMatrix; }set { _uMatrix = value; }}public int[] PivotArray{get { return _pivotArray; }set { _pivotArray = value; }}}
This class is used to store the results of the LU decomposition and includes the L matrix, U matrix, pivot array, and the status of the LU decomposition operation.
The LUDecompose function is given below.
private LUDecompositionResults LUDecompose(){Debug.Assert(_columnCount == _rowCount);// Using Crout Decomp with P//// Ax = b //By definition of problem variables.//// LU = PA //By definition of L, U, and P.//// LUx = Pb //By substition for PA.//// Ux = d //By definition of d//// Ld = Pb //By subsitition for d.////For 4x4 with P = I// [l11 0 0 0 ] [1 u12 u13 u14] [a11 a12 a13 a14]// [l21 l22 0 0 ] [0 1 u23 u24] = [a21 a22 a23 a24]// [l31 l32 l33 0 ] [0 0 1 u34] [a31 a32 a33 a34]// [l41 l42 l43 l44] [0 0 0 1 ] [a41 a42 a43 a44]LUDecompositionResults result = new LUDecompositionResults();try{int[] pivotArray = new int[_rowCount]; //Pivot matrix.Matrix uMatrix = new Matrix(_rowCount, _columnCount);Matrix lMatrix = new Matrix(_rowCount, _columnCount);Matrix workingUMatrix = Clone();Matrix workingLMatrix = new Matrix(_rowCount, _columnCount);for(int i=0;i<_rowCount;i++){pivotArray[i] = i;}//Iterate down the number of rows in the U matrix.for (int i = 0; i < _rowCount; i++){//Do pivots first.//I want to make the matrix diagnolaly dominate.//Initialize the variables used to determine the pivot row.double maxRowRatio = double.NegativeInfinity;int maxRow = -1;int maxPosition = -1;//Check all of the rows below and including the current row//to determine which row should be pivoted to the working row position.//The pivot row will be set to the row with the maximum ratio//of the absolute value of the first column element divided by the//sum of the absolute values of the elements in that row.for (int j = i; j < _rowCount; j++){//Store the sum of the absolute values of the row elements in//dRowSum. Clear it out now because I am checking a new row.double rowSum = 0.0;//Go across the columns, add the absolute values of the elements in//that column to dRowSum.for (int k = i; k < _columnCount; k++){rowSum += Math.Abs(workingUMatrix[pivotArray[j], k]);}//Check to see if the absolute value of the ratio of the lead//element over the sum of the absolute values of the elements is larger//that the ratio for preceding rows. If it is, then the current row//becomes the new pivot candidate.if (rowSum == 0.0){throw new SingularMatrixException();}double dCurrentRatio = Math.Abs(workingUMatrix[pivotArray[j], i]) / rowSum;lock (this){if (dCurrentRatio > maxRowRatio){maxRowRatio = Math.Abs(workingUMatrix[pivotArray[j], i] / rowSum);maxRow = pivotArray[j];maxPosition = j;}}}//If the pivot candidate isn't the current row, update the//pivot array to swap the current row with the pivot row.if (maxRow != pivotArray[i]){int hold = pivotArray[i];pivotArray[i] = maxRow;pivotArray[maxPosition] = hold;}//Store the value of the left most element in the working U//matrix in dRowFirstElementValue.double rowFirstElementValue = workingUMatrix[pivotArray[i], i];//Update the columns of the working row. j is the column index.for(int j=0;j<_columnCount;j++){if (j < i){//If j<1, then the U matrix element value is 0.workingUMatrix[pivotArray[i], j] = 0.0;}else if (j == i){//If i == j, the L matrix value is the value of the//element in the working U matrix.workingLMatrix[pivotArray[i], j] = rowFirstElementValue;//The value of the U matrix for i == j is 1workingUMatrix[pivotArray[i], j] = 1.0;}else // j>i{//Divide each element in the current row of the U matrix by the//value of the first element in the rowworkingUMatrix[pivotArray[i], j] /= rowFirstElementValue;//The element value of the L matrix for j>i is 0workingLMatrix[pivotArray[i], j] = 0.0;}}//For the working U matrix, subtract the ratioed active row from the rows below it.//Update the columns of the rows below the working row. k is the row index.for (int k = i + 1; k < _rowCount; k++){//Store the value of the first element in the working row//of the U matrix.rowFirstElementValue = workingUMatrix[pivotArray[k], i];//Go accross the columns of row k.for(int j=0;j<_columnCount;j++){if (j < i){//If j<1, then the U matrix element value is 0.workingUMatrix[pivotArray[k], j] = 0.0;}else if (j == i){//If i == j, the L matrix value is the value of the//element in the working U matrix.workingLMatrix[pivotArray[k], j] = rowFirstElementValue;//The element value of the L matrix for j>i is 0workingUMatrix[pivotArray[k], j] = 0.0;}else //j>i{workingUMatrix[pivotArray[k], j] = workingUMatrix[pivotArray[k], j] - rowFirstElementValue * workingUMatrix[pivotArray[i], j];}}}}for (int i = 0; i < _rowCount; i++){for (int j = 0; j < _rowCount; j++){uMatrix[i, j] = workingUMatrix[pivotArray[i], j];lMatrix[i, j] = workingLMatrix[pivotArray[i], j];}}result.U = uMatrix;result.L = lMatrix;result.PivotArray = pivotArray;}catch (AggregateException ex2){if (ex2.InnerExceptions.Count > 0){throw ex2.InnerExceptions[0];}else{throw ex2;}}catch (Exception ex3){throw ex3;}return result;}
This function uses two matrices mUWorking and mLWorking for intermediate calculations.  The matrix mUWorking is initially set to be a clone of the matrix calling LUDecompose.  And it is reduced to upper triangular form during the calculations.
The first for loop of i iterates across all of the rows in the matrix. 
For each iteration, it finds which of the remaining rows to pivot to the pivoted row.  The pivot row is determined by which row has the largest ratio of the absolute value of the element in the working column divided by the sum of the absolute values of all of the elements in that row.  After that row is found, it then becomes the pivoted row. 
After the pivoted row is found the function then iterates over the columns of the pivoted row.  If the column number is less than the pivot column number, then the mUWorking element is set to 0.  If the column number is equal to the pivot column number then the element in mLWorking is set to the first element value in the pivoted row and the element in mUWorking is set to 1.  For column numbers greater than the pivot column number, the element value of mUWorking is divided by the original element value of the pivot column.
After processing the pivoted row, the function then process the other rows below the pivoted row in the pivot table.  For these rows, if the column number is less than the pivot column number, then the element value of  mUWorking is set to 0.  If the column number is equal to the pivot column number then element value of mLWorking is set to the pivot column element value of mUWorking and the mUWorking element value is set to 0.  For column numbers greater than the pivot column number, the pivoted row is multiplied by the first pivot column element of the row and then subtracted from the row.
After mUWorking and mLWorking have been solved, the element values of them are copied into mU and mL in the correct sort order.  mU, mL, and the pivot array are then packaged into the LUDecompositionResults, which is returned from the function.
The rest of the code for the SolveFor method is shown below.
public Matrix SolveFor(Matrix rightMatrix){Debug.Assert(rightMatrix.RowCount == _columnCount);Debug.Assert(_columnCount == _rowCount);Matrix resultMatrix = new Matrix(_columnCount, rightMatrix.ColumnCount);LUDecompositionResults resDecomp = LUDecompose();int[] nP = resDecomp.PivotArray;Matrix lMatrix = resDecomp.L;Matrix uMatrix = resDecomp.U;for (int k = 0; k < rightMatrix.ColumnCount; k++){//Solve for the corresponding d Matrix from Ld=Pbdouble sum = 0.0;Matrix dMatrix = new Matrix(_rowCount, 1);dMatrix[0, 0] = rightMatrix[nP[0], k] / lMatrix[0, 0];for (int i = 1; i < _rowCount; i++){sum = 0.0;for (int j = 0; j < i; j++){sum += lMatrix[i, j] * dMatrix[j, 0];}dMatrix[i, 0] = (rightMatrix[nP[i], k] - sum) / lMatrix[i, i];}//Solve for x using Ux = dresultMatrix[_rowCount - 1, k] = dMatrix[_rowCount - 1, 0];for (int i = _rowCount - 2; i >= 0; i--){sum = 0.0;for (int j = i + 1; j < _rowCount; j++){sum += uMatrix[i, j] * resultMatrix[j, k];}resultMatrix[i, k] = dMatrix[i, 0] - sum;}}return resultMatrix;}
This code decomposes the matrix into the L and U matrix and the pivot array by calling LUDecompose.  It then uses backsubstitution to solve for d.  After solving for d, it uses back substitution to solve for x.
Now that matrix equations can be solved with LU decomposition with partial pivoting, computing the inverse of a matrix is rather simple.  By definition, the inverse of a matrix A fulfills the following equation AA-1 =I.  This implies that to invert the matrix, all I need to do is call SolveFor on the matrix and pass in an identity matrix with the same number of rows and columns.  The code for that is shown below.
public Matrix Invert(){Debug.Assert(_rowCount == _columnCount);Matrix resultMatrix = SolveFor(Identity(_rowCount));Matrix matIdent = this * resultMatrix;return SolveFor(Identity(_rowCount));}
The full code for the matrix class at this point, including that from parts 1 and 2 is given below.
using System.Diagnostics;
using System;namespace NumericalMethods.ThirdBlog
{public class Matrix{#region ctorpublic Matrix(){_values = new double[_rowCount, _columnCount];}public Matrix(int rowCount, int columnCount){_rowCount = rowCount;_columnCount = columnCount;_values = new double[_rowCount, _columnCount];}#endregion#region Row Column valuespublic double this[int row, int column]{get { return _values[row, column]; }set { _values[row, column] = value; }}#endregion#region F&Pprivate double[,] _values;private int _rowCount = 3;public int RowCount{get { return _rowCount; }}private int _columnCount = 3;public int ColumnCount{get { return _columnCount; }}#endregion#region basic single matrix stuffpublic static Matrix Identity(int size){Matrix resultMatrix = new Matrix(size, size);for (int i = 0; i < size; i++){for (int j = 0; j < size; j++){resultMatrix[i, j] = (i == j) ? 1.0 : 0.0;}}return resultMatrix;}public Matrix Clone(){Matrix resultMatrix = new Matrix(_rowCount, _columnCount);for (int i = 0; i < _rowCount; i++){for (int j = 0; j < _columnCount; j++){resultMatrix[i, j] = this[i, j];}}return resultMatrix;}public Matrix Transpose(){Matrix resultMatrix = new Matrix(_columnCount, _rowCount);for (int i = 0; i < _rowCount; i++){for (int j = 0; j < _columnCount; j++){resultMatrix[j, i] = this[i, j];}}return resultMatrix;}#endregion#region Binary Mathpublic static Matrix Add(Matrix leftMatrix, Matrix rightMatrix){Debug.Assert(leftMatrix.ColumnCount == rightMatrix.ColumnCount);Debug.Assert(leftMatrix.RowCount == rightMatrix.RowCount);Matrix resultMatrix = new Matrix(leftMatrix.RowCount, rightMatrix.ColumnCount);for (int i = 0; i < leftMatrix.RowCount; i++){for (int j = 0; j < leftMatrix.ColumnCount; j++){resultMatrix[i, j] = leftMatrix[i, j] + rightMatrix[i, j];}}return resultMatrix;}public static Matrix operator +(Matrix leftMatrix, Matrix rightMatrix){return Matrix.Add(leftMatrix, rightMatrix);}public static Matrix Subtract(Matrix leftMatrix, Matrix rightMatrix){Debug.Assert(leftMatrix.ColumnCount == rightMatrix.ColumnCount);Debug.Assert(leftMatrix.RowCount == rightMatrix.RowCount);Matrix resultMatrix = new Matrix(leftMatrix.RowCount, rightMatrix.ColumnCount);for (int i = 0; i < leftMatrix.RowCount; i++){for (int j = 0; j < leftMatrix.ColumnCount; j++){resultMatrix[i, j] = leftMatrix[i, j] - rightMatrix[i, j];}}return resultMatrix;}public static Matrix operator -(Matrix leftMatrix, Matrix rightMatrix){return Matrix.Subtract(leftMatrix, rightMatrix);}public static Matrix Multiply(Matrix leftMatrix, Matrix rightMatrix){Debug.Assert(leftMatrix.ColumnCount == rightMatrix.RowCount);Matrix resultMatrix = new Matrix(leftMatrix.RowCount, rightMatrix.ColumnCount);for (int i = 0; i < resultMatrix.ColumnCount; i++){for (int j = 0; j < leftMatrix.RowCount; j++){double value = 0.0;for (int k = 0; k < rightMatrix.RowCount; k++){value += leftMatrix[j, k] * rightMatrix[k, i];}resultMatrix[j, i] = value;}}return resultMatrix;}public static Matrix operator *(Matrix leftMatrix, Matrix rightMatrix){return Matrix.Multiply(leftMatrix, rightMatrix);}public static Matrix Multiply(double left, Matrix rightMatrix){Matrix resultMatrix = new Matrix(rightMatrix.RowCount, rightMatrix.ColumnCount);for (int i = 0; i < resultMatrix.RowCount; i++){for (int j = 0; j < rightMatrix.ColumnCount; j++){resultMatrix[i, j] = left * rightMatrix[i, j];}}return resultMatrix;}public static Matrix operator *(double left, Matrix rightMatrix){return Matrix.Multiply(left, rightMatrix);}public static Matrix Multiply(Matrix leftMatrix, double right){Matrix resultMatrix = new Matrix(leftMatrix.RowCount, leftMatrix.ColumnCount);for (int i = 0; i < leftMatrix.RowCount; i++){for (int j = 0; j < leftMatrix.ColumnCount; j++){resultMatrix[i, j] = leftMatrix[i, j] * right;}}return resultMatrix;}public static Matrix operator *(Matrix leftMatrix, double right){return Matrix.Multiply(leftMatrix, right);}public static Matrix Divide(Matrix leftMatrix, double right){Matrix resultMatrix = new Matrix(leftMatrix.RowCount, leftMatrix.ColumnCount);for (int i = 0; i < leftMatrix.RowCount; i++){for (int j = 0; j < leftMatrix.ColumnCount; j++){resultMatrix[i, j] = leftMatrix[i, j] / right;}}return resultMatrix;}public static Matrix operator /(Matrix leftMatrix, double right){return Matrix.Divide(leftMatrix, right);}#endregion#region Assorted Castspublic static Matrix FromArray(double[] left){int length = left.Length;Matrix resultMatrix = new Matrix(length, 1);for (int i = 0; i < length; i++){resultMatrix[i, 0] = left[i];}return resultMatrix;}public static implicit operator Matrix(double[] left){return FromArray(left);}public static double[] ToArray(Matrix leftMatrix){Debug.Assert((leftMatrix.ColumnCount == 1 && leftMatrix.RowCount >= 1) || (leftMatrix.RowCount == 1 && leftMatrix.ColumnCount >= 1));double[] result = null;if (leftMatrix.ColumnCount > 1){int numElements = leftMatrix.ColumnCount;result = new double[numElements];for (int i = 0; i < numElements; i++){result[i] = leftMatrix[0, i];}}else{int numElements = leftMatrix.RowCount;result = new double[numElements];for (int i = 0; i < numElements; i++){result[i] = leftMatrix[i, 0];}}return result;}public static implicit operator double[](Matrix leftMatrix){return ToArray(leftMatrix);}public static Matrix FromDoubleArray(double[,] left){int length0 = left.GetLength(0);int length1 = left.GetLength(1);Matrix resultMatrix = new Matrix(length0, length1);for (int i = 0; i < length0; i++){for (int j = 0; j < length1; j++){resultMatrix[i, j] = left[i, j];}}return resultMatrix;}public static implicit operator Matrix(double[,] left){return FromDoubleArray(left);}public static double[,] ToDoubleArray(Matrix leftMatrix){double[,] result = new double[leftMatrix.RowCount, leftMatrix.ColumnCount];for (int i = 0; i < leftMatrix.RowCount; i++){for (int j = 0; j < leftMatrix.ColumnCount; j++){result[i, j] = leftMatrix[i, j];}}return result;}public static implicit operator double[,](Matrix leftMatrix){return ToDoubleArray(leftMatrix);}#endregionpublic Matrix SolveFor(Matrix rightMatrix){Debug.Assert(rightMatrix.RowCount == _columnCount);Debug.Assert(_columnCount == _rowCount);Matrix resultMatrix = new Matrix(_columnCount, rightMatrix.ColumnCount);LUDecompositionResults resDecomp = LUDecompose();int[] nP = resDecomp.PivotArray;Matrix lMatrix = resDecomp.L;Matrix uMatrix = resDecomp.U;for (int k = 0; k < rightMatrix.ColumnCount; k++){//Solve for the corresponding d Matrix from Ld=Pbdouble sum = 0.0;Matrix dMatrix = new Matrix(_rowCount, 1);dMatrix[0, 0] = rightMatrix[nP[0], k] / lMatrix[0, 0];for (int i = 1; i < _rowCount; i++){sum = 0.0;for (int j = 0; j < i; j++){sum += lMatrix[i, j] * dMatrix[j, 0];}dMatrix[i, 0] = (rightMatrix[nP[i], k] - sum) / lMatrix[i, i];}//Solve for x using Ux = dresultMatrix[_rowCount - 1, k] = dMatrix[_rowCount - 1, 0];for (int i = _rowCount - 2; i >= 0; i--){sum = 0.0;for (int j = i + 1; j < _rowCount; j++){sum += uMatrix[i, j] * resultMatrix[j, k];}resultMatrix[i, k] = dMatrix[i, 0] - sum;}}return resultMatrix;}private LUDecompositionResults LUDecompose(){Debug.Assert(_columnCount == _rowCount);// Using Crout Decomp with P//// Ax = b //By definition of problem variables.//// LU = PA //By definition of L, U, and P.//// LUx = Pb //By substition for PA.//// Ux = d //By definition of d//// Ld = Pb //By subsitition for d.////For 4x4 with P = I// [l11 0 0 0 ] [1 u12 u13 u14] [a11 a12 a13 a14]// [l21 l22 0 0 ] [0 1 u23 u24] = [a21 a22 a23 a24]// [l31 l32 l33 0 ] [0 0 1 u34] [a31 a32 a33 a34]// [l41 l42 l43 l44] [0 0 0 1 ] [a41 a42 a43 a44]LUDecompositionResults result = new LUDecompositionResults();try{int[] pivotArray = new int[_rowCount]; //Pivot matrix.Matrix uMatrix = new Matrix(_rowCount, _columnCount);Matrix lMatrix = new Matrix(_rowCount, _columnCount);Matrix workingUMatrix = Clone();Matrix workingLMatrix = new Matrix(_rowCount, _columnCount);for(int i=0;i<_rowCount;i++){pivotArray[i] = i;}//Iterate down the number of rows in the U matrix.for (int i = 0; i < _rowCount; i++){//Do pivots first.//I want to make the matrix diagnolaly dominate.//Initialize the variables used to determine the pivot row.double maxRowRatio = double.NegativeInfinity;int maxRow = -1;int maxPosition = -1;//Check all of the rows below and including the current row//to determine which row should be pivoted to the working row position.//The pivot row will be set to the row with the maximum ratio//of the absolute value of the first column element divided by the//sum of the absolute values of the elements in that row.for (int j = i; j < _rowCount; j++){//Store the sum of the absolute values of the row elements in//dRowSum. Clear it out now because I am checking a new row.double rowSum = 0.0;//Go across the columns, add the absolute values of the elements in//that column to dRowSum.for (int k = i; k < _columnCount; k++){rowSum += Math.Abs(workingUMatrix[pivotArray[j], k]);}//Check to see if the absolute value of the ratio of the lead//element over the sum of the absolute values of the elements is larger//that the ratio for preceding rows. If it is, then the current row//becomes the new pivot candidate.if (rowSum == 0.0){throw new SingularMatrixException();}double dCurrentRatio = Math.Abs(workingUMatrix[pivotArray[j], i]) / rowSum;lock (this){if (dCurrentRatio > maxRowRatio){maxRowRatio = Math.Abs(workingUMatrix[pivotArray[j], i] / rowSum);maxRow = pivotArray[j];maxPosition = j;}}}//If the pivot candidate isn't the current row, update the//pivot array to swap the current row with the pivot row.if (maxRow != pivotArray[i]){int hold = pivotArray[i];pivotArray[i] = maxRow;pivotArray[maxPosition] = hold;}//Store the value of the left most element in the working U//matrix in dRowFirstElementValue.double rowFirstElementValue = workingUMatrix[pivotArray[i], i];//Update the columns of the working row. j is the column index.for(int j=0;j<_columnCount;j++){if (j < i){//If j<1, then the U matrix element value is 0.workingUMatrix[pivotArray[i], j] = 0.0;}else if (j == i){//If i == j, the L matrix value is the value of the//element in the working U matrix.workingLMatrix[pivotArray[i], j] = rowFirstElementValue;//The value of the U matrix for i == j is 1workingUMatrix[pivotArray[i], j] = 1.0;}else // j>i{//Divide each element in the current row of the U matrix by the//value of the first element in the rowworkingUMatrix[pivotArray[i], j] /= rowFirstElementValue;//The element value of the L matrix for j>i is 0workingLMatrix[pivotArray[i], j] = 0.0;}}//For the working U matrix, subtract the ratioed active row from the rows below it.//Update the columns of the rows below the working row. k is the row index.for (int k = i + 1; k < _rowCount; k++){//Store the value of the first element in the working row//of the U matrix.rowFirstElementValue = workingUMatrix[pivotArray[k], i];//Go accross the columns of row k.for(int j=0;j<_columnCount;j++){if (j < i){//If j<1, then the U matrix element value is 0.workingUMatrix[pivotArray[k], j] = 0.0;}else if (j == i){//If i == j, the L matrix value is the value of the//element in the working U matrix.workingLMatrix[pivotArray[k], j] = rowFirstElementValue;//The element value of the L matrix for j>i is 0workingUMatrix[pivotArray[k], j] = 0.0;}else //j>i{workingUMatrix[pivotArray[k], j] = workingUMatrix[pivotArray[k], j] - rowFirstElementValue * workingUMatrix[pivotArray[i], j];}}}}for (int i = 0; i < _rowCount; i++){for (int j = 0; j < _rowCount; j++){uMatrix[i, j] = workingUMatrix[pivotArray[i], j];lMatrix[i, j] = workingLMatrix[pivotArray[i], j];}}result.U = uMatrix;result.L = lMatrix;result.PivotArray = pivotArray;}catch (AggregateException ex2){if (ex2.InnerExceptions.Count > 0){throw ex2.InnerExceptions[0];}else{throw ex2;}}catch (Exception ex3){throw ex3;}return result;}public Matrix Invert(){Debug.Assert(_rowCount == _columnCount);Matrix resultMatrix = SolveFor(Identity(_rowCount));Matrix matIdent = this * resultMatrix;return SolveFor(Identity(_rowCount));}}public class LUDecompositionResults{private Matrix _lMatrix;private Matrix _uMatrix;private int[] _pivotArray;public LUDecompositionResults(){}public LUDecompositionResults(Matrix matL, Matrix matU, int[] nPivotArray){_lMatrix = matL;_uMatrix = matU;_pivotArray = nPivotArray;}public Matrix L{get { return _lMatrix; }set { _lMatrix = value; }}public Matrix U{get { return _uMatrix; }set { _uMatrix = value; }}public int[] PivotArray{get { return _pivotArray; }set { _pivotArray = value; }}}public class SingularMatrixException : ArithmeticException{public SingularMatrixException(): base("Invalid operation on a singular matrix."){}}
}