5.2 LU Factorization

In the previous section, the actual procedure for solving the set of linear equations Ax = b uses the LU factorization technique.

If A is a square matrix, it has an LU factorization with a permutation matrix P, an upper triangular matrix U, and a lower triangular matrix L, such that A = PLU. This factorization can be used to solve linear equations Ax = b. It can also be used to compute the determinant of A from the equation det(A) = det(P)det(L)det(U), and to compute the inverse of A as A-1 = U-1L-1P-1.

Let us see how Math.h++ does this. The basic idea is that an LU factorization is actually implemented as a class in the Math.h++ class library. This is an example of the general notion in object-oriented programming that if you can think of something as an "it," you can probably make "it" into a class. The LU factorization class is called RWGenFact<T>, and it takes a matrix as an argument in its constructor:

template<class T>
RWGenFact<T>(const RWGenMat<T>&);

Once constructed, this LU factorization can be used for many things. For example, there is the global function inverse():

template<class T>
RWGenMat<T> inverse(const RWGenFact<T>&);

which takes an argument of type RWGenFact<T>. Hence, one could invert a matrix as follows:

RWGenMat<double> A;
.
.
.
// Construct the LU factorization:
RWGenFact<double> ALU(A); 

// Now invert it:
RWGenMat<double> Ainv = inverse(ALU);

The inverse() global function is overloaded to take either RWGenMat<T> or RWGenFact<T>. When inverse() is given an RWGenMat<T>, it creates a temporary RWGenFact<T> and passes that to inverse(). In the first example above, this is what is actually happening:

RWGenMat<double> a;

// Type conversion occurs to matrix a:
RWGenMat<double> ainv = inverse(a);

Similar to inverse(), the global function solve() is also overloaded:

RWMathVec<T> solve (const RWGenFact<T>&, const RWMathVec<T>&);
RWMathVec<T> solve (const RWGenMat<T>&, const RWMathVec<T>&);

When called with an argument of type RWGenMat<T>, a similar type conversion occurs.

The great advantage to explicitly creating an RWGenFact<T> is that you can create an LU factorization only once and use it many times, perhaps to solve many sets of linear equations:

RWGenMat<double> a;                       // Matrix a
RWMathVec<double> rhs[5];                 // Set of 5 vectors
RWMathVec<double> soln[5];                // Solution vectors
.     
.     // (Initialize the matrix a and vectors rhs)
.
// Construct the LU factorization of a:
RWGenFact<double> aLU(a);

// Solve 5 sets of equations using the same LU factorization:
for(int i = 0; i<5; i++)
   soln[i] = solve(aLU, rhs[i]);

This approach can result in great time savings because calculating the LU factorization takes most of the time, but need be done only once. Yet its calculation becomes transparent to the user who is not interested in the factorization (or wants to remain ignorant of it!)