LinAlgebraPacks
Packages for Linear Algebra: Solving your system of equations using (optimised, bug-free) code that someone else has already written!
Introduction
svn co https://svn.ggy.bris.ac.uk/subversion-open/num-methods1 ./num-methods1LAPACK:direct solutions for dense matrices
For example, consider the following system of linear equations:
- [math]\displaystyle{ \begin{alignat}{7} x &&\; + \;&& 3y &&\; - \;&& 2z &&\; = \;&& 5 & \\ 3x &&\; + \;&& 5y &&\; + \;&& 6z &&\; = \;&& 7 & \\ 2x &&\; + \;&& 4y &&\; + \;&& 3z &&\; = \;&& 8 & \end{alignat} }[/math]
(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:
- [math]\displaystyle{ \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 3 & 5 & 6 & 7 \\ 2 & 4 & 3 & 8 \end{array}\right] }[/math][math]\displaystyle{ \sim \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 0 & -4 & 12 & -8 \\ 2 & 4 & 3 & 8 \end{array}\right] }[/math][math]\displaystyle{ \sim \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 0 & -4 & 12 & -8 \\ 0 & -2 & 7 & -2 \end{array}\right] }[/math][math]\displaystyle{ \sim \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 0 & 1 & -3 & 2 \\ 0 & -2 & 7 & -2 \end{array}\right] }[/math][math]\displaystyle{ \sim \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 0 & 1 & -3 & 2 \\ 0 & 0 & 1 & 2 \end{array}\right] }[/math][math]\displaystyle{ \sim \left[\begin{array}{rrr|r} 1 & 3 & -2 & 5 \\ 0 & 1 & 0 & 8 \\ 0 & 0 & 1 & 2 \end{array}\right] }[/math][math]\displaystyle{ \sim \left[\begin{array}{rrr|r} 1 & 3 & 0 & 9 \\ 0 & 1 & 0 & 8 \\ 0 & 0 & 1 & 2 \end{array}\right] }[/math][math]\displaystyle{ \sim \left[\begin{array}{rrr|r} 1 & 0 & 0 & -15 \\ 0 & 1 & 0 & 8 \\ 0 & 0 & 1 & 2 \end{array}\right]. }[/math]
To solve this using LAPACK take a look inside dgesv-example.f90. First the declarations:
integer, parameter :: N = 3
integer, parameter :: LDA = N ! leading dimension of A
integer, parameter :: LDB = N ! leading dimension of B
integer, parameter :: NRHS = 1 ! no. of RHS, i.e columns in b
integer, dimension(N) :: IPIV
integer :: INFO
integer :: i
real(kind=8), dimension(LDA,N) :: A ! LDAxN matrix
real(kind=8), dimension(LDB,NRHS) :: B ! LDBxNRHS matrix
Next, we set up the matrix A:
! ( 1 3 -2)
! ( 3 5 6)
! ( 2 4 3)
A(1,1) = 1.0d+0
A(1,2) = 3.0d+0
A(1,3) = -2.0d+0
A(2,1) = 3.0d+0
A(2,2) = 5.0d+0
A(2,3) = 6.0d+0
A(3,1) = 2.0d+0
A(3,2) = 4.0d+0
A(3,3) = 3.0d+0
And the vector B:
! ( 5)
! ( 7)
! ( 8)
B(1,1) = 5.0d+0
B(2,1) = 7.0d+0
B(3,1) = 8.0d+0
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):
call dgesv(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
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. The reason for this is that for the same computational cost, we can factor A without knowing about b. As a result we can 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?
- [math]\displaystyle{ A'=L_1 A }[/math]
- [math]\displaystyle{ A = LU }[/math],
where U is an upper triangular matrix, and L is a unit (1's on the diagonal) lower triangular matrix.
Given L and U we can solve Ax=b by first solving:
- [math]\displaystyle{ Ly=b }[/math].
As L is unit lower triangular, we can do this by easily by forward substitution. Once we know y, we can solve:
- [math]\displaystyle{ Ux=y }[/math],
easily again, this time by backward substitution.
Computing the LU decomposition of an [math]\displaystyle{ n\timesn }[/math] system of equations has a computational cost of [math]\displaystyle{ O(n^3) }[/math] operations, whereas the forward and backward substitution steps have a cost of <math>O(n^2)<math> operations. Thus the cost of solution is dominated by the LU decomposition.
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.