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-methods1

LAPACK:direct solutions for dense matrices

http://www.netlib.org/lapack

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

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