LU decomposition

== Description ==
Consider the system of equations <math>A\vec{x}=\vec{b}</math>, where <math>A</math> is an <math>n\times n</math> nonsingular matrix.  <math>A</math> may be decomposed into an lower triangular part <math>L</math> and an upper triangular part <math>U</math> that will lead us to a direct procedure for the solution of the original system.  This decomposition procedure is quite useful when more than one right-hand side (more than one <math>\vec{b}</math>) is to be used.

The algorithm is relatively straightforward - first, we determine the upper and lower triangular parts:

:<math>A = LU.</math>

Then,

:<math> A \vec{x} = (LU) \vec{x} = L(U\vec{x})=L\vec{y},</math>

where <math>y=Ux</math>.  Once we solve the system 

:<math>L\vec{y} = \vec{b},</math> 

we will be able to find the solution to the original system by solving

:<math>U\vec{x} = \vec{y}</math>

The first solution is a foward substitution, while the second solution is a backward substitution.  Both can be done efficiently once the factorization is available.  The forward substitution may be expressed as

:<math>
y_i  = {1 \over {l_{ii} }}\left( {b_i  - \sum\limits_{j = 1}^{i-1} {l_{ij} y_j } } \right),
i = 1,2,\ldots,n</math> 

and the backward substitution process may be expressed as

:<math>
x_i  = {1 \over {u_{ii} }}\left( {y_i  - \sum\limits_{j = i + 1}^n {u_{ij} x_j } } \right),
i = n,n-1,\ldots,1.</math>


LU decomposition essentially stores the operations of Gaussian elimination in "higher-level" form (see [[#References|Golub and Van Loan]]), so repeated solutions using the same left-hand side are computed without repetition of operations that are independent of the right-hand side.

== Algorithm ==
''Add factorization here''<br>
Forward substitution <br>
:     for k:=1 step until n do <br>
::     for i:=1 step until k-1 <br>
:::     <math>b_k=b_k-l_{ki}b_{i}</math> <br>
::     end loop (i) <br>
::    <math>b_{k}=b_{k}/l_{kk}</math><br>
:      end loop (k)
Backward substitution
:     for k:=n stepdown until 1 do <br>
::     for i:=k+1 step until n <br>
:::     <math>b_k=b_k-u_{ki}b_{i}</math> <br>
::     end loop (i) <br>
::    <math>x_{k}=b_{k}/u_{kk}</math><br>
:      end loop (k)

== Important Considerations ==
As with [[Gaussian elimination]],  LU decomposition is probably best used for relatively small, relatively non-sparse systems of equations (with small and non-sparse open to some interpretation).  For larger and/or sparse problems, it would probably be best to either use an iterative method or use a direct solver package (e.g. [http://www.cse.psu.edu/~raghavan/software.html DSCPACK]) as opposed to writing one of your own.

If one has a single lefthand-side matrix and many right-hand side vectors, then LU decomposition would be a good solution procedure to consider.  At the very least, it should be faster than solving each system separately with Gaussian elimination.

Here a Fortran code fragment for LU decomposition written by  D. Partenov

 do j = 1, n
  do i = 1, j
    if (i == 1) then
      a(i,j) = a(i, j)
    else
       suma = 0.0
      do k = 1, i-1
         suma = suma + a(i,k)*a(k,j)
      end do
      a(i,j) = a(i,j) - suma
    end if
  end do
  if ( j < n) then
    do s = 1, n-j
      i = j + s
      if (j == 1) then
        a(i,j) = a(i,j)/a(j,j)
      else
        suma = 0.0
        do k = 1, j-1
          suma = suma + a(i,k)*a(k,j)
        end do
        a(i,j) = (a(i,j) - suma)/a(j,j)
      end if
    end do
  end if
 end do
 y(1) = b(1)
 do i = 2, n
  suma = 0.0
  do j = 1, i-1
    suma = suma + a(i,j)*y(j)
  end do
  y(i) = b(i) - suma
 end do
 i = n
 j = n
 x(i) = y(i)/a(i,j)
 do s = 1, n -1
  i = n - s
  suma = 0.0
  do j = i+1, n
    suma = suma + a(i,j)*x(j)
  end do
  x(i) = (y(i) - suma)/a(i,i)
 end do

== References ==
*{{reference-book|author=Golub and Van Loan|year=1996|title=Matrix Computations|rest= 3rd edition, The Johns Hopkins University Press, Baltimore}}
*[http://en.wikipedia.org/wiki/LU_decomposition Wikipedia article ''LU decomposition]