Computational Science I

1
Computational Science I

• Lecture 13
• CAAM 420, Fall 2004
• Instructor Tim Warburton

2
Class Exercise cont

• Download the BlasDemo.tar.gz tarball from the
  website and
• make sure you can do the following
• gzip d BlasDemo.tar.gz
• tar xvf BlasDemo.tar
• cd BlasDemo
• make blas1
• ./myblaslevel1
• Modify the myblaslevel1.c file to compute
• c a.(alphaa b) using ddot, and daxpy
• c b.(alphaa) using ddot, and dscal

3
cont

• How would we compute the followingusing BLAS

4
LAPACK

• Last class we discussed the BLAS subroutines for
  basic vector-vector, matrix-vector and
  matrix-matrix computations.
• BLAS is used as the building block for the Linear
  Algegra Package, LAPACK.
• LAPACK can be downloaded from
• http//www.netlib.org/lapack/
• An online LAPACK manual can be downloaded from
• http//www.netlib.org/lapack/lug/index.html
• This library consists of a set of higher level
  linear algebra functions with interface described
  at
• http//www.netlib.org/lapack/individualroutines.ht
  ml

5
LAPACK cont

• There are a very large number of linear algebra
  subroutines available in LAPACK
• I would like to focus on a subset of very
  commonly used subroutines
• dgetrf which is used to compute LU
  factorizations of a matrix
• dgetrs which uses an LU factorization from dgetrf
  to solve a system
• dgetri which uses the LU above to compute the
  inverse of a matrix
• dgesv, this is essentially a combined call to
  dgetrf and dgetrs
• dgeev, this computes the eigenvalues of a matrix.
• Note each of dgetrf, dgesv, and dgeev are prone
  to writing over your original matrix.

6
dgetrfdouble general triangular factorization

• This routine computes an LU factorization of a
  matrix with partial pivoting.
• There is model code at
• http//www.netlib.org/lapack/double/dgetrf.f
• This entails computing a lower triangular matrix
  L and an upper triangular matrix U and a
  permutation matrix P such that
• The two matrices L and U are stored in the space
  originally occupied by A
• The pivot matrix (i.e. permutation matrix) P is
  stored as a column vector of integers.
• See http//www.math.unm.edu/timwar/MA375S04/Lectu
  re19.pdf for description of the LU algorithm.

7
dgetrf cont
Sample invocation of dgetrf. Notice I have added
an int pointer, ipivot, to the dmat struct
It is important to stress that we now have the
ability to LU factorize, and it only took
theconstruction of this interface i.e. we did
not have to reinvent the wheel.
8
Philosophy Time

• Heres the difficult part for a lot of people
• We do not need to know how dgetrf performs the LU
  factorization.
• We only need to know the input and output for
  dgetrf
• Beyond that we just need to know how to interpret
  and/or use the output from dgetrf
• This should not be a big worry to you. Most of
  you drive a car, however I doubt most of you know
  much about the inner workings of a modern engine.
  However, you do know that to make the car go from
  A to B you make sure it has enough gas, take off
  the parking brake and hit the gas possibly
  checking your mirrors.
• Caveat if dgetrf breaks down then we may need
  to resort to the users manual. Specifically,
  dgetrf returns an info code and possibly error
  messages as we discussed before these are
  pivotal in resolving issues.

9
dgetrsdouble general triangular solve

• dgetrs uses the output from dgetrf , i.e. the LU
  factorization of a matrix to find the matrix X
  which satisfies (PLU)\X B for a given B
• Again, we do not need to know how dgetrs works.
• So I am not going to tell you how it works!.
• Think about that I am giving you the keys to a
  car, but not lecturing you on the inner-workings
  of the car.

10
dgetrs cont
Again, I have wrapped a complicated LAPACK call
in a more accessible dmat function
11
dgetrf and dgetrs in action

• In my dmat code I hide the dgetrf and dgetrs
  calls and instead use
• / compute the LU factorization of A /
• PLU dmatLUfactorize(A)
• / backsolve using the LU factorization of i.e.
  H(PLU)\A /
• H dmatLUbacksolve(PLU, A)

12
dgesvdouble general solve

• This one wraps the dgetrf and dgetrs routines
  together

Version with dgetrf,dgetrs
Version with dgesv
13
if else... endif

• Notice on the previous piece of code I used a
  conditional statement passed directly to the
  compiler.
• This is not C-code, it is a directive for the
  compiler.
• In this instance
• if 1
• task1
• else
• task2
• endif
• Will cause task1 to be compiled, and not task2
• This can be useful for instance you can use
  this to exclude some of your error checking
  routines for when the code is running and
  debugged, i.e. in what is called production mode.

14
Class Exercise

• Download LapackDemo.tar.gz from the website
• make mylapack
• ./mylapack
• Examine the mylapack.c file
• Based on mylapack.c write code to find all the
  eigenvalues of
• Is it obvious how to code this? Did you need to
  know about dgeev?