## LU Matrix Decomposition in C++

To download:

1 | bzr branch http://bazaar.enseed.com/lib/Generic/ look under Geometry/LUDecomposition.h |

A LU decomposition decomposes a matrix **M** into a **L**ower and **U**pper matrix pair such that **L**×**U**=**M**. This representation is useful in many ways – for example, it’s easy to solve (**L**×**U**)x = **b**. (for more info, see LU Decomposition on Wikipedia or LU Decomposition on MathWorld.

This particular implementation uses the Matrix class and the same philosophy. Since our matrix class has the size specified as a template parameter, the LU decomposition knows the size of the matrix at compile time. Like most LU decomposition algorithms, it will use a recursive method (For a NxN matrix, first solve for a (N^{-1})xN matrix…) except that this is done here using recursive templates. The compiler can thus generate the entire unrolled decomposition.

## Decomposing a Matrix

Decomposing a matrix in its **LU** pair is done like this:

1 2 3 4 | Matrix4d m; Matrix4d l, u; LUDecomposition::decompose(m, &l, &u); // l is now a lower matrix, u is an upper matrix, l*u = m |

## Solving Mx=b

You once you have a **LU** pair, you can solve **M**x = b = **(L×U)**x = b:

1 2 3 | Matrix xb; // set this to b initially, you'll get x after LUDecomposition::solve(l, u, &xb); // xb now holds x |

## Matrix Inverse

In the previous example, x and b were vectors (or 1×4 matrices), but they could be nx4 matrices. For example, say we have **MX** =**B** where **X** = **M ^{-1}** and

**B**=

**I**, solving for

**X**would be solving for the inverse of our matrix

**M**:

1 2 3 | Matrix4d inverse = Matrix4d::identityMatrix(); // set this to b initially, you'll get x after LUDecomposition::solve(l, u, &xb); // xb now holds x, which is the inverse of m (or l*u) |

## Being More Robust

Sometimes, matrices won’t be well formed and the LU decomposition will generate large numerical errors. One way to prevent this is to solve for **LU** = **PM** where **P** is a pivot matrix (it changes the orders of the rows in **M**. When the ultimate goal is to use the decomposition to solve a system, we can live with this solution and obtain better precision:

1 2 3 4 5 | Matrix4d m; Matrix4d l, u; Point4i p; // pivot LUDecomposition::decompose(m, &l, &u, &p); // we now have l*u = p*m. |

Now, instead of solving **M**x = b, we’ll have to solve for **PM**x = **P**b = **LU**x. In other words, we have to solve for **P**b instead of just b. Say we want to inverse our matrix like in the previous example:

1 2 3 4 5 6 7 8 9 | Matrix4d pivotMatrix = Matrix4d::zeroMatrix(); for (int i = 0; i < 4; ++i) pivotMatrix[i][p[i]] = 1; Matrix4d xb = Matrix4d::itentityMatrix(); // this is b xb = pivotMatirx * xb; // this is Pb LUDecomposition::solve(l, u, &xb); // xb is now x, which is the inverse of our matrix |

Hello,

can anybody help me out? I’m a beginner in programming. I’m using MS Visual Studio 2008.

My purpose (example):

Solve M * x = b,

where M is a 7*7 matrix with complex elements (except column 7) and b is a 1*7 matrix with real elements.

How to use LUDecomposition.h?

1. copy this header to header files

2. write a source file (“test.cpp)

3. run

Right prosedure? Did I forgot anything? How to include header to test.cpp? How to use complex numbers?

Any hints for me?

Sorry for the language, I’m german…

Never mind…

I found a solution that works for me on myself.

Cu