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 solution of dense matrices
For example, consider the following system:
- [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 of the above system of equations:
- [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:
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
Solving a 2D laplacian (a 4x4 domain):
./simple-laplace.exe
And also (using a more flexibly sized square domain):
./mxm-laplace.exe
Direct Solve using 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.