Notes Linear Algebra � Assignment 1 is out (questions?) � Last class: we reduced the problem of “optimally” interpolating scattered data to solving a system of linear equations � Typical: often almost all of the computational work in a scientific computing code is linear algebra operations cs542g-term1-2007 1 cs542g-term1-2007 2 Basic Definitions BLAS � Many common matrix/vector operations have � Matrix/vector notation been standardized into an API called the BLAS � Dot product, outer product (Basic Linear Algebra Subroutines) • Level 1: vector operations � Vector norms copy, scale, dot, add, norms, … � Matrix norms • Level 2: matrix-vector operations multiply, triangular solve, … • Level 3: matrix-matrix operations multiply, triangular solve, … � FORTRAN bias, but callable from other langs � Goals: • As fast as possible, but still safe/accurate � www.netlib.org/blas cs542g-term1-2007 3 cs542g-term1-2007 4 Speed in BLAS LAPACK � In each level: � The BLAS only solves triangular systems • Forward or backward substitution multithreading, prefetching, vectorization, loop unrolling, etc. � LAPACK is a higher level API for matrix operations: � In level 2, especially in level 3: blocking • Solving linear systems • Operate on sub-blocks of the matrix that fit the • Solving linear least squares problems memory architecture well • Solving eigenvalue problems � General goal: � Built on the BLAS, with blocking in mind to keep if it � s easy to phrase an operation in terms high performance of BLAS, get speed+safety for free � Biggest advantage: safety • The higher the level better • Designed to handle difficult problems gracefully � www.netlib.org/lapack cs542g-term1-2007 5 cs542g-term1-2007 6
Specializations Accuracy � When solving a linear system, first question to � Before jumping into algorithms, how ask: what sort of system? accurate can we hope to be in solving a � Many properties to consider: linear system? • Single precision or double? � Key idea: backward error analysis • Real or complex? � Assume calculated answer is the • Invertible or (nearly) singular? exact solution of a perturbed problem . • Symmetric/Hermitian? • Definite or Indefinite? • Dense or sparse or specially structured? � The condition number of a problem: • Multiple right-hand sides? how much errors in data get amplified in � LAPACK/BLAS take advantage of many of these solution (sparse matrices the big exception…) cs542g-term1-2007 7 cs542g-term1-2007 8 Condition Number Gaussian Elimination � Sometimes we can estimate the condition � Let � s start with the simplest unspecialized number of a matrix a priori algorithm: Gaussian Elimination � Special case: for a symmetric matrix, � Assume the matrix is invertible, but 2-norm condition number is ratio of otherwise nothing special known about it extreme eigenvalues � GE simply is row-reduction to upper triangular form, followed by backwards substitution � LAPACK also provides cheap estimates • Permuting rows if we run into a zero • Try to construct a vector ||x|| that comes close to maximizing ||A -1 x|| cs542g-term1-2007 9 cs542g-term1-2007 10 LU Factorization Block Approach to LU � Each step of row reduction is multiplication � Rather than get bogged down in details of by an elementary matrix GE (hard to see forest for trees) � Gathering these together, we find GE is � Partition the equation A=LU essentially a matrix factorization: � Gives natural formulas for algorithms A=LU where � Extends to block algorithms L is lower triangular (and unit diagonal), U is upper triangular � Solving Ax=b by GE is then Ly=b Ux=y cs542g-term1-2007 11 cs542g-term1-2007 12
Cholesky Factorization � If A is symmetric positive definite, can cut work in half: A=LL T • L is lower triangular � If A is symmetric but indefinite, possibly still have the Modified Cholesky factorization: A=LDL T • L is unit lower triangular • D is diagonal cs542g-term1-2007 13
Recommend
More recommend
