Revision as of 23:42, 12 February 2011

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

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 [math]\displaystyle{ \mathbf{A} }[/math] with a unit (1's on the diagonal) lower triangular matrix, [math]\displaystyle{ \mathbf{L_1} }[/math], with non-zeros only in the first column, such as that schematically shown below:

[math]\displaystyle{ \begin{bmatrix} 1 & 0 & 0 \\ x & 1 & 0 \\ x & 0 & 1 \end{bmatrix} }[/math]

Following our example, [math]\displaystyle{ \mathbf{L_1 A=A'} }[/math] becomes:

[math]\displaystyle{ \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} }[/math]

The eliminations in the second column are achieved by left multiplying [math]\displaystyle{ \mathbf{A'=L_1 A} }[/math] with [math]\displaystyle{ \mathbf{L_2} }[/math] (which is unit lower triangular with non-zeros only in the second column):

[math]\displaystyle{ \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} }[/math]

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

[math]\displaystyle{ \mathbf{A=LU} }[/math].

Given L and U we can solve [math]\displaystyle{ \mathbf{Ax=b} }[/math] by first solving:

[math]\displaystyle{ \mathbf{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{ \mathbf{Ux=y} }[/math],

easily again, this time by backward substitution.

Computing the LU decomposition of an [math]\displaystyle{ n \times n }[/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]\displaystyle{ O(n^2) }[/math] 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:

do ii=2,4
  ...
  call dgetrs('N', N, NRHS, A, LDA, IPIV, B, LDB, INFO)
  ...
end do

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

do ii=2,4
  ...
  ! copy B into B1, so that we can recompute for
  ! the boundary conditions again.
  B1 = B

  ! Need to recalculate entries for the vector of known
  ! values--held in B--given the solution at the current
  ! timestep--held in B1.
  ! *NB* the loops below are only valid for the case r=1.
  ! LHS
  B(1,1) = 0 + B1(2,1)
  ! RHS
  B(N,1) = 0 + B1(N-1,1)
  ! middle bit
  do jj=2,N-1
     B(jj,1) = B1(jj-1,1) + B1(jj+1,1)
  end do
  ...
end do

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

Difference between revisions of "LinAlgebraPacks"