LinAlgebraPacks

Packages for Linear Algebra: Solving your system of equations using (optimised, bug-free) code that someone else has already written!

=Introduction=

=LAPACK:direct solutions for dense matrices=

http://www.netlib.org/lapack

For example, consider the following system of linear equations:
 * $$\begin{alignat}{7}

x &&\; + \;&& 3y &&\; - \;&& 2z &&\; = \;&& 5 & \\ 3x &&\; + \;&& 5y &&\; + \;&& 6z &&\; = \;&& 7 & \\ 2x &&\; + \;&& 4y &&\; + \;&& 3z &&\; = \;&& 8 & \end{alignat}$$

(From http://en.wikipedia.org/wiki/System_of_linear_equations.)

The following computation shows Gauss-Jordan elimination applied to the matrix version (Ax=b) of the above system:
 * $$\left[\begin{array}{rrr|r}

1 & 3 & -2 & 5 \\ 3 & 5 & 6 & 7 \\ 2 & 4 & 3 & 8 \end{array}\right]$$$$\sim \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 0 & -4 & 12 & -8 \\ 2 & 4 & 3 & 8 \end{array}\right]$$$$\sim \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 0 & -4 & 12 & -8 \\ 0 & -2 & 7 & -2 \end{array}\right]$$$$\sim \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 0 & 1 & -3 & 2 \\ 0 & -2 & 7 & -2 \end{array}\right]$$$$\sim \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 0 & 1 & -3 & 2 \\ 0 & 0 & 1 & 2 \end{array}\right]$$$$\sim \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 0 & 1 & 0 & 8 \\ 0 & 0 & 1 & 2 \end{array}\right]$$$$\sim \left[\begin{array}{rrr|r} 1 & 3 & 0 & 9 \\ 0 & 1 & 0 & 8 \\ 0 & 0 & 1 & 2 \end{array}\right]$$$$\sim \left[\begin{array}{rrr|r} 1 & 0 & 0 & -15 \\ 0 & 1 & 0 & 8 \\ 0 & 0 & 1 & 2 \end{array}\right].$$

To solve this using LAPACK take a look inside dgesv-example.f90. First the declarations:

Next, we set up the matrix A:

And the vector B:

Once we have all that in place, all we need to do is call the solver routine. (Note that the solution vector replaces the values in B):

Run the program to satisfy yourself that it works:

cd num-methods1/examples/example1 make ./dgesv-example.exe

Next, take a look at simple-laplace.f90. This program solves the 2D laplace equation for a 4x4 domain, as described on NumMethodsPDEs. This program only differs from the previous one in way in which we initialise the matrix A and vector b. To run the program type:

./simple-laplace.exe

The last program in this first example is mxm-laplace.f90. In this case, we have just added a little more complexity, so that a 2D laplace equation can be solved on a more flexibly sized--mxm--square domain:

./mxm-laplace.exe

=Direct Solve using LU Decomposition=

In practice, the LAPACK routine dgesv does not perform Gaussian elimination, but rather computes an LU decomposition and then solves the equation using two substitution steps--one forward and one backward. This is because for the same computational cost, we can factor A without knowing about b. We can then find x cheaply, for many different b.  This is useful when solving parabolic equations, where x at time-step j forms the basis for b at time-step j+1.

OK, how to we obtain an LU decomposition?

The Guassian elimination process above can be expressed as a matrix multiplication. The elimination in the first column can be achieved by left multiplying A with a unit (1's on the diagonal) lower triangular matrix, $$L_1$$, with non-zeros only in the first column, such as that schematically shown below:

$$ \begin{bmatrix} 1 & 0 & 0 \\ x & 1 & 0 \\ x & 0 & 1 \end{bmatrix} $$

Following our example, $$L_1 A=A'$$ becomes:

$$ \begin{bmatrix} 1 & 0 & 0 \\ -3 & 1 & 0 \\ -2 & 0 & 1 \end{bmatrix} \times \begin{bmatrix} 1 & 3 & -2 \\ 3 & 5 & 6 \\ 2 & 4 & 3 \end{bmatrix} = \begin{bmatrix} 1 & 3 & -2 \\ 0 & -4 & 12 \\ 0 & -2 & 7 \end{bmatrix} $$

The eliminations in the second column are achieved by left multiplying $$A'=L_1 A$$ with $$L_2$$ (which is unit lower triangular with non-zeros only in the second column):

$$ \begin{bmatrix} 1 & 0           & 0 \\ 0 & 1            & 0 \\ 0 & -\frac{1}{2} & 1 \end{bmatrix} \times \begin{bmatrix} 1 & 0 & 0 \\ -3 & 1 & 0 \\ -2 & 0 & 1 \end{bmatrix} \times \begin{bmatrix} 1 & 3 & -2 \\ 3 & 5 & 6 \\ 2 & 4 & 3 \end{bmatrix} = \begin{bmatrix} 1 & 3 & -2 \\ 0 & -4 & 12 \\ 0 & 0 & 1 \end{bmatrix} $$

We note that the result $$L_2 L_1 A = A'' = U$$ is upper triangular. Rearranging we get, $$A=L_1^{-1} \cdots L^{-1}_{n-2} L^{-1}_{n-1} U$$, where the inverse matrices are easy to find and the product is itself unit lower triangular. Thus we may write:


 * $$A=LU$$.

Given L and U we can solve $$Ax=b$$ by first solving:


 * $$Ly=b$$.

As L is unit lower triangular, we can do this by easily by forward substitution. Once we know y, we can solve:


 * $$Ux=y$$,

easily again, this time by backward substitution.

Computing the LU decomposition of an $$n \times n$$ system of equations has a computational cost of $$O(n^3)$$ operations, whereas the forward and backward substitution steps have a cost of $$O(n^2)$$ operations. Thus the cost of solution is dominated by the LU decomposition.

Our second example uses a pre-computed LU decomposition when solving a parabolic equation using the Crank-Nicolson scheme.

cd ../example2

In simple-crank-nicolson.f90 we see that the initial call to dgesv replaces the entries in the matrix A with an encoding of the LU decompostion (the unit diagonal is assumed). For subsequent timesteps we can thus call dgetrs instead:

where we have remembered to take acount of the boundary conditions at each time-step too:

=PetSc: Iterative solutions and sparse matrices=

http://www.mcs.anl.gov/petsc/petsc-as/

cd ../example3 make ./ex1f

For further examples see INSTALL_DIR/petsc-3.1-p1/src/ksp/ksp/examples/tutorials.

=PLASMA: For Multicore Architectures=

http://icl.cs.utk.edu/plasma/index.html

cd ../example4 make ./mxm-laplace-plasma.exe

For further examples, see e.g. INSTALL_DIR/plasma-installer_2.3.1/build/plasma_2.3.1.

=MAGMA: To include GPUs in Heterogeneous Systems=

http://icl.cs.utk.edu/magma/index.html