**LinearAlgebra****`****GaussianElimination****`**

The Mathematica function LinearSolve computes the solution to a system of linear equations. There are two distinct steps to this process, and LinearSolve does both of them. There are cases where you want to solve several systems of equations with the same left-hand side in each case, but with different right-hand sides. In these cases LinearSolve will certainly work, but since the first step of the solution process deals only with the left-hand side it will be doing much of the work repeatedly. To make matters worse, for large systems the first step is most of the work anyway.

Because the first step can be seen abstractly as factoring the matrix of coefficients into the product of a lower triangular and an upper triangular matrix, it is often referred to as LU factorization. LUFactor performs this factorization and returns a data object with a head of LU. This data object represents the factorization of the original matrix together with information regarding any permutation of the rows that might be necessary to maintain numerical stability in the computation.

The second step in the process is often referred to as back substitution or back-solving. LUSolve performs this operation.

Solving a linear system using Gaussian elimination and back substitution.

The data object generated by LUFactor is a compact way of storing the information contained in the upper and lower triangular matrices of the factorization. The zeros in the lower part of the upper triangular matrix, the upper part of the lower triangular matrix, and the ones along the diagonal of the lower triangular matrix are not stored explicitly. This allows both the upper and lower triangular matrices to be stored in the space of a single square matrix. The actual storage arrangement is even more complicated than this because the rows usually get permuted in an effort to make the solution process as numerically stable as possible.

This loads the package.
In[1]:= **<<LinearAlgebra`GaussianElimination`**

Here is the coefficient matrix of a system.
In[2]:= **MatrixForm[a = {{5, 3, 0}, {7, 9, 2},**

{-2, -8, -1}}]

Out[2]//MatrixForm=

This performs the LU factorization.
In[3]:= **lu = LUFactor[a]**

Out[3]=

This is the right-hand side of the system.
In[4]:= **b = {6, -3, 7}**

Out[4]=

LUSolve uses back substitution to get the solution to the system.
In[5]:= **LUSolve[lu, b]**

Out[5]=

We can also use back substitution to solve a system with a different right-hand side.
In[6]:= **LUSolve[lu, {1,2,3}]**

Out[6]=

This checks that the solution is correct.
In[7]:= **a . % - {1,2,3}**

Out[7]=

The functions LUFactor and LUSolve can be used in the context of other forms of arithmetic such as interval arithmetic or NumericalMath`ComputerArithmetic`.

This converts the coefficient matrix and the right-hand side to have elements that are intervals rather than numbers.
In[8]:= **a = Map[Interval[{#-.01,#+.01}]&, a, {2}]; b = Interval[{#-.01,#+.01}]& /@ b**

Out[8]=

Now we solve the system in the context of interval arithmetic.
In[9]:= **LUSolve[LUFactor[a], b]**

Out[9]=

This checks the solution.
In[10]:= **N[a . %] - b**

Out[10]=