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 MatrixNumeric SolveFor(MatrixNumeric mRight)
{
}
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 MatrixNumeric m_matL;
private MatrixNumeric m_matU;
private int[] m_nPivotArray;
private LUDecompositionResultStatus m_enuStatus = LUDecompositionResultStatus.OK;
public LUDecompositionResults()
{
}
public LUDecompositionResults(MatrixNumeric matL, MatrixNumeric matU, int[] nPivotArray, LUDecompositionResultStatus enuStatus)
{
m_matL = matL;
m_matU = matU;
m_nPivotArray = nPivotArray;
m_enuStatus = enuStatus;
}
public MatrixNumeric L
{
get { return m_matL; }
set { m_matL = value; }
}
public MatrixNumeric U
{
get { return m_matU; }
set { m_matU = value; }
}
public int[] PivotArray
{
get { return m_nPivotArray; }
set { m_nPivotArray = value; }
}
public LUDecompositionResultStatus Status
{
get { return m_enuStatus; }
set { m_enuStatus = value; }
}
}
public enum LUDecompositionResultStatus
{
OK = 0,
LinearlyDependent = 1
}
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(m_nNumColumns == m_nNumRows);
// 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 resRet = new LUDecompositionResults();
int[] nP = new int[m_nNumRows]; //Pivot matrix.
MatrixNumeric mU = new MatrixNumeric(m_nNumRows, m_nNumColumns);
MatrixNumeric mL = new MatrixNumeric(m_nNumRows, m_nNumColumns);
MatrixNumeric mUWorking = Clone();
MatrixNumeric mLWorking = new MatrixNumeric(m_nNumRows, m_nNumColumns);
for (int i = 0; i < m_nNumRows; i++)
{
nP[i] = i;
}
//Iterate down the number of rows in the U matrix.
for (int i = 0; i < m_nNumRows; i++)
{
//Do pivots first.
//I want to make the matrix diagnolaly dominate.
//Initialize the variables used to determine the pivot row.
double dMaxRowRatio = double.NegativeInfinity;
int nMaxRow = -1;
int nMaxPosition = -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 < m_nNumRows; 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 dRowSum = 0.0;
//Go across the columns, add the absolute values of the elements in
//that column to dRowSum.
for (int k = i; k < m_nNumColumns; k++)
{
dRowSum += Math.Abs(mUWorking[nP[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 (Math.Abs(mUWorking[nP[j], i] / dRowSum) > dMaxRowRatio)
{
dMaxRowRatio = Math.Abs(mUWorking[nP[j], i] / dRowSum);
nMaxRow = nP[j];
nMaxPosition = j;
}
}
//If the pivot candidate isn't the current row, update the
//pivot array to swap the current row with the pivot row.
if (nMaxRow != nP[i])
{
int nHold = nP[i];
nP[i] = nMaxRow;
nP[nMaxPosition] = nHold;
}
//Store the value of the left most element in the working U
//matrix in dRowFirstElementValue.
double dRowFirstElementValue = mUWorking[nP[i], i];
//Update the columns of the working row. j is the column index.
for (int j = 0; j < m_nNumRows; j++)
{
if (j < i)
{
//If j<1, then the U matrix element value is 0.
mUWorking[nP[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.
mLWorking[nP[i], j] = dRowFirstElementValue;
//The value of the U matrix for i == j is 1
mUWorking[nP[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 row
mUWorking[nP[i], j] /= dRowFirstElementValue;
//The element value of the L matrix for j>i is 0
mLWorking[nP[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 < m_nNumRows; k++)
{
//Store the value of the first element in the working row
//of the U matrix.
dRowFirstElementValue = mUWorking[nP[k], i];
//Go accross the columns of row k.
for (int j = 0; j < m_nNumRows; j++)
{
if (j < i)
{
//If j<1, then the U matrix element value is 0.
mUWorking[nP[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.
mLWorking[nP[k], j] = dRowFirstElementValue;
//The element value of the L matrix for j>i is 0
mUWorking[nP[k], j] = 0.0;
}
else //j>i
{
mUWorking[nP[k], j] = mUWorking[nP[k], j] - dRowFirstElementValue * mUWorking[nP[i], j];
}
}
}
}
for (int i = 0; i < m_nNumRows; i++)
{
for (int j = 0; j < m_nNumRows; j++)
{
mU[i, j] = mUWorking[nP[i], j];
mL[i, j] = mLWorking[nP[i], j];
}
}
resRet.U = mU;
resRet.L = mL;
resRet.PivotArray = nP;
return resRet;
}
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 MatrixNumeric SolveFor(MatrixNumeric mRight)
{
Debug.Assert(mRight.NumRows == m_nNumColumns);
Debug.Assert(m_nNumColumns == m_nNumRows);
MatrixNumeric mRet =
new MatrixNumeric(m_nNumColumns, mRight.NumColumns);
LUDecompositionResults resDecomp = LUDecompose();
int[] nP = resDecomp.PivotArray;
MatrixNumeric mL = resDecomp.L;
MatrixNumeric mU = resDecomp.U;
double dSum = 0.0;
for (int k = 0; k < mRight.NumColumns; k++)
{
//Solve for the corresponding d Matrix from Ld=Pb
MatrixNumeric D = new MatrixNumeric(m_nNumRows, 1);
D[0, 0] = mRight[nP[0], k] / mL[0, 0];
for (int i = 1; i < m_nNumRows; i++)
{
dSum = 0.0;
for (int j = 0; j < i; j++)
{
dSum += mL[i, j] * D[j, 0];
}
D[i, 0] = (mRight[nP[i], k] - dSum) / mL[i, i];
}
//Solve for x using Ux = d
//DMatrix X = new DMatrix(m_nNumRows, 1);
mRet[m_nNumRows - 1, k] = D[m_nNumRows - 1, 0];
for (int i = m_nNumRows - 2; i >= 0; i--)
{
dSum = 0.0;
for (int j = i + 1; j < m_nNumRows; j++)
{
dSum += mU[i, j] * mRet[j, k];
}
mRet[i, k] = D[i, 0] - dSum;
}
}
return mRet;
}
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 MatrixNumeric Invert()
{
System.Diagnostics.Debug.Assert(m_nNumRows == m_nNumColumns);
return SolveFor(Identity(m_nNumRows));
}
The full code for the matrix class at this point, including that from parts 1 and 2 is given below.
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;
namespace NumericalMethods
{
public class MatrixNumeric
{
private int m_nNumRows = 3;
private int m_nNumColumns = 3;
private double[,] m_dValues;
public MatrixNumeric()
{
m_dValues = new double[m_nNumRows, m_nNumColumns];
}
public MatrixNumeric(int nNumRows, int nNumColumns)
{
m_nNumRows = nNumRows;
m_nNumColumns = nNumColumns;
m_dValues = new double[m_nNumRows, m_nNumColumns];
}
public double this[int nRow, int nColumn]
{
get { return m_dValues[nRow, nColumn]; }
set { m_dValues[nRow, nColumn] = value; }
}
public int NumRows
{
get { return m_nNumRows; }
}
public int NumColumns
{
get { return m_nNumColumns; }
}
public MatrixNumeric Clone()
{
MatrixNumeric mRet = new MatrixNumeric(m_nNumRows, m_nNumColumns);
for (int i = 0; i < m_nNumRows; i++)
{
for (int j = 0; j < m_nNumColumns; j++)
{
mRet[i, j] = this[i, j];
}
}
return mRet;
}
public static MatrixNumeric Identity(int nSize)
{
MatrixNumeric mRet = new MatrixNumeric(nSize, nSize);
for (int i = 0; i < nSize; i++)
{
for (int j = 0; j < nSize; j++)
{
mRet[i, j] = (i == j) ? 1.0 : 0.0;
}
}
return mRet;
}
public MatrixNumeric Transpose()
{
MatrixNumeric mRet = new MatrixNumeric(m_nNumColumns, m_nNumRows);
for (int i = 0; i < m_nNumRows; i++)
{
for (int j = 0; j < m_nNumColumns; j++)
{
mRet[j, i] = this[i, j];
}
}
return mRet;
}
public static MatrixNumeric FromArray(double[] dLeft)
{
int nLength = dLeft.Length;
MatrixNumeric mRet = new MatrixNumeric(nLength, 1);
for (int i = 0; i < nLength; i++)
{
mRet[i, 0] = dLeft[i];
}
return mRet;
}
public static implicit operator MatrixNumeric(double[] dLeft)
{
return FromArray(dLeft);
}
public static double[] ToArray(MatrixNumeric mLeft)
{
Debug.Assert((mLeft.NumColumns == 1 && mLeft.NumRows >= 1) || (mLeft.NumRows == 1 && mLeft.NumColumns >= 1));
double[] dRet = null;
if (mLeft.NumColumns > 1)
{
int nNumElements = mLeft.NumColumns;
dRet = new double[nNumElements];
for (int i = 0; i < nNumElements; i++)
{
dRet[i] = mLeft[0, i];
}
}
else
{
int nNumElements = mLeft.NumRows;
dRet = new double[nNumElements];
for (int i = 0; i < nNumElements; i++)
{
dRet[i] = mLeft[i, 0];
}
}
return dRet;
}
public static implicit operator double[](MatrixNumeric mLeft)
{
return ToArray(mLeft);
}
public static MatrixNumeric FromDoubleArray(double[,] dLeft)
{
int nLength0 = dLeft.GetLength(0);
int nLength1 = dLeft.GetLength(1);
MatrixNumeric mRet = new MatrixNumeric(nLength0, nLength1);
for (int i = 0; i < nLength0; i++)
{
for (int j = 0; j < nLength1; j++)
{
mRet[i, j] = dLeft[i, j];
}
}
return mRet;
}
public static implicit operator MatrixNumeric(double[,] dLeft)
{
return FromDoubleArray(dLeft);
}
public static double[,] ToDoubleArray(MatrixNumeric mLeft)
{
double[,] dRet = new double[mLeft.NumRows, mLeft.NumColumns];
for (int i = 0; i < mLeft.NumRows; i++)
{
for (int j = 0; j < mLeft.NumColumns; j++)
{
dRet[i, j] = mLeft[i, j];
}
}
return dRet;
}
public static implicit operator double[,](MatrixNumeric mLeft)
{
return ToDoubleArray(mLeft);
}
public static MatrixNumeric Add(MatrixNumeric mLeft, MatrixNumeric mRight)
{
Debug.Assert(mLeft.NumColumns == mRight.NumColumns);
Debug.Assert(mLeft.NumRows == mRight.NumRows);
MatrixNumeric mRet = new MatrixNumeric(mLeft.NumRows, mRight.NumColumns);
for (int i = 0; i < mLeft.NumRows; i++)
{
for (int j = 0; j < mLeft.NumColumns; j++)
{
mRet[i, j] = mLeft[i, j] + mRight[i, j];
}
}
return mRet;
}
public static MatrixNumeric operator +(MatrixNumeric mLeft, MatrixNumeric mRight)
{
return MatrixNumeric.Add(mLeft, mRight);
}
public static MatrixNumeric Subtract(MatrixNumeric mLeft, MatrixNumeric mRight)
{
Debug.Assert(mLeft.NumColumns == mRight.NumColumns);
Debug.Assert(mLeft.NumRows == mRight.NumRows);
MatrixNumeric mRet = new MatrixNumeric(mLeft.NumRows, mRight.NumColumns);
for (int i = 0; i < mLeft.NumRows; i++)
{
for (int j = 0; j < mLeft.NumColumns; j++)
{
mRet[i, j] = mLeft[i, j] - mRight[i, j];
}
}
return mRet;
}
public static MatrixNumeric operator -(MatrixNumeric mLeft, MatrixNumeric mRight)
{
return MatrixNumeric.Subtract(mLeft, mRight);
}
public static MatrixNumeric Multiply(MatrixNumeric mLeft, MatrixNumeric mRight)
{
Debug.Assert(mLeft.NumColumns == mRight.NumRows);
MatrixNumeric mRet = new MatrixNumeric(mLeft.NumRows, mRight.NumColumns);
for (int i = 0; i < mRight.NumColumns; i++)
{
for (int j = 0; j < mLeft.NumRows; j++)
{
double dValue = 0.0;
for (int k = 0; k < mRight.NumRows; k++)
{
dValue += mLeft[j, k] * mRight[k, i];
}
mRet[j, i] = dValue;
}
}
return mRet;
}
public static MatrixNumeric operator *(MatrixNumeric mLeft, MatrixNumeric mRight)
{
return MatrixNumeric.Multiply(mLeft, mRight);
}
public static MatrixNumeric Multiply(double dLeft, MatrixNumeric mRight)
{
MatrixNumeric mRet = new MatrixNumeric(mRight.NumRows, mRight.NumColumns);
for (int i = 0; i < mRight.NumRows; i++)
{
for (int j = 0; j < mRight.NumColumns; j++)
{
mRet[i, j] = dLeft * mRight[i, j];
}
}
return mRet;
}
public static MatrixNumeric operator *(double dLeft, MatrixNumeric mRight)
{
return MatrixNumeric.Multiply(dLeft, mRight);
}
public static MatrixNumeric Multiply(MatrixNumeric mLeft, double dRight)
{
MatrixNumeric mRet = new MatrixNumeric(mLeft.NumRows, mLeft.NumColumns);
for (int i = 0; i < mLeft.NumRows; i++)
{
for (int j = 0; j < mLeft.NumColumns; j++)
{
mRet[i, j] = mLeft[i, j] * dRight;
}
}
return mRet;
}
public static MatrixNumeric operator *(MatrixNumeric mLeft, double dRight)
{
return MatrixNumeric.Multiply(mLeft, dRight);
}
public static MatrixNumeric Divide(MatrixNumeric mLeft, double dRight)
{
MatrixNumeric mRet = new MatrixNumeric(mLeft.NumRows, mLeft.NumColumns);
for (int i = 0; i < mLeft.NumRows; i++)
{
for (int j = 0; j < mLeft.NumColumns; j++)
{
mRet[i, j] = mLeft[i, j] / dRight;
}
}
return mRet;
}
public static MatrixNumeric operator /(MatrixNumeric mLeft, double dRight)
{
return MatrixNumeric.Divide(mLeft, dRight);
}
public MatrixNumeric SolveFor(MatrixNumeric mRight)
{
Debug.Assert(mRight.NumRows == m_nNumColumns);
Debug.Assert(m_nNumColumns == m_nNumRows);
MatrixNumeric mRet = new MatrixNumeric(m_nNumColumns, mRight.NumColumns);
LUDecompositionResults resDecomp = LUDecompose();
int[] nP = resDecomp.PivotArray;
MatrixNumeric mL = resDecomp.L;
MatrixNumeric mU = resDecomp.U;
double dSum = 0.0;
for (int k = 0; k < mRight.NumColumns; k++)
{
//Solve for the corresponding d Matrix from Ld=Pb
MatrixNumeric D = new MatrixNumeric(m_nNumRows, 1);
D[0, 0] = mRight[nP[0], k] / mL[0, 0];
for (int i = 1; i < m_nNumRows; i++)
{
dSum = 0.0;
for (int j = 0; j < i; j++)
{
dSum += mL[i, j] * D[j, 0];
}
D[i, 0] = (mRight[nP[i], k] - dSum) / mL[i, i];
}
//Solve for x using Ux = d
//DMatrix X = new DMatrix(m_nNumRows, 1);
mRet[m_nNumRows - 1, k] = D[m_nNumRows - 1, 0];
for (int i = m_nNumRows - 2; i >= 0; i--)
{
dSum = 0.0;
for (int j = i + 1; j < m_nNumRows; j++)
{
dSum += mU[i, j] * mRet[j, k];
}
mRet[i, k] = D[i, 0] - dSum;
}
}
return mRet;
}
private LUDecompositionResults LUDecompose()
{
Debug.Assert(m_nNumColumns == m_nNumRows);
// 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 resRet = new LUDecompositionResults();
int[] nP = new int[m_nNumRows]; //Pivot matrix.
MatrixNumeric mU = new MatrixNumeric(m_nNumRows, m_nNumColumns);
MatrixNumeric mL = new MatrixNumeric(m_nNumRows, m_nNumColumns);
MatrixNumeric mUWorking = Clone();
MatrixNumeric mLWorking = new MatrixNumeric(m_nNumRows, m_nNumColumns);
for (int i = 0; i < m_nNumRows; i++)
{
nP[i] = i;
}
//Iterate down the number of rows in the U matrix.
for (int i = 0; i < m_nNumRows; i++)
{
//Do pivots first.
//I want to make the matrix diagnolaly dominate.
//Initialize the variables used to determine the pivot row.
double dMaxRowRatio = double.NegativeInfinity;
int nMaxRow = -1;
int nMaxPosition = -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 < m_nNumRows; 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 dRowSum = 0.0;
//Go across the columns, add the absolute values of the elements in
//that column to dRowSum.
for (int k = i; k < m_nNumColumns; k++)
{
dRowSum += Math.Abs(mUWorking[nP[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 (Math.Abs(mUWorking[nP[j], i] / dRowSum) > dMaxRowRatio)
{
dMaxRowRatio = Math.Abs(mUWorking[nP[j], i] / dRowSum);
nMaxRow = nP[j];
nMaxPosition = j;
}
}
//If the pivot candidate isn't the current row, update the
//pivot array to swap the current row with the pivot row.
if (nMaxRow != nP[i])
{
int nHold = nP[i];
nP[i] = nMaxRow;
nP[nMaxPosition] = nHold;
}
//Store the value of the left most element in the working U
//matrix in dRowFirstElementValue.
double dRowFirstElementValue = mUWorking[nP[i], i];
//Update the columns of the working row. j is the column index.
for (int j = 0; j < m_nNumRows; j++)
{
if (j < i)
{
//If j<1, then the U matrix element value is 0.
mUWorking[nP[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.
mLWorking[nP[i], j] = dRowFirstElementValue;
//The value of the U matrix for i == j is 1
mUWorking[nP[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 row
mUWorking[nP[i], j] /= dRowFirstElementValue;
//The element value of the L matrix for j>i is 0
mLWorking[nP[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 < m_nNumRows; k++)
{
//Store the value of the first element in the working row
//of the U matrix.
dRowFirstElementValue = mUWorking[nP[k], i];
//Go accross the columns of row k.
for (int j = 0; j < m_nNumRows; j++)
{
if (j < i)
{
//If j<1, then the U matrix element value is 0.
mUWorking[nP[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.
mLWorking[nP[k], j] = dRowFirstElementValue;
//The element value of the L matrix for j>i is 0
mUWorking[nP[k], j] = 0.0;
}
else //j>i
{
mUWorking[nP[k], j] = mUWorking[nP[k], j] - dRowFirstElementValue * mUWorking[nP[i], j];
}
}
}
}
for (int i = 0; i < m_nNumRows; i++)
{
for (int j = 0; j < m_nNumRows; j++)
{
mU[i, j] = mUWorking[nP[i], j];
mL[i, j] = mLWorking[nP[i], j];
}
}
resRet.U = mU;
resRet.L = mL;
resRet.PivotArray = nP;
return resRet;
}
public MatrixNumeric Invert()
{
Debug.Assert(m_nNumRows == m_nNumColumns);
return SolveFor(Identity(m_nNumRows));
}
}
public class LUDecompositionResults
{
private MatrixNumeric m_matL;
private MatrixNumeric m_matU;
private int[] m_nPivotArray;
private LUDecompositionResultStatus m_enuStatus = LUDecompositionResultStatus.OK;
public LUDecompositionResults()
{
}
public LUDecompositionResults(MatrixNumeric matL, MatrixNumeric matU, int[] nPivotArray, LUDecompositionResultStatus enuStatus)
{
m_matL = matL;
m_matU = matU;
m_nPivotArray = nPivotArray;
m_enuStatus = enuStatus;
}
public MatrixNumeric L
{
get { return m_matL; }
set { m_matL = value; }
}
public MatrixNumeric U
{
get { return m_matU; }
set { m_matU = value; }
}
public int[] PivotArray
{
get { return m_nPivotArray; }
set { m_nPivotArray = value; }
}
public LUDecompositionResultStatus Status
{
get { return m_enuStatus; }
set { m_enuStatus = value; }
}
}
public enum LUDecompositionResultStatus
{
OK = 0,
LinearlyDependent = 1
}
}
Subscribe