Notes Linear Algebra Assignment 1 will be out later today Last - PDF document

Notes Linear Algebra � Assignment 1 will be out later today � Last class: (look on the web) we reduced the problem of “optimally” interpolating scattered data to solving a system of linear equations � This week: start delving into numerical linear algebra � Often almost all of the computational work in a scientific computing code is linear algebra operations cs542g-term1-2006 1 cs542g-term1-2006 2 Basic Definitions Accuracy � Matrix/vector notation � How accurate can we expect a floating point matrix-vector multiply to be? � Dot product, outer product • Assume result is the exact answer to a � Vector norms perturbed problem � Matrix norms � How accurate are real implementations? cs542g-term1-2006 3 cs542g-term1-2006 4 BLAS Speed in BLAS � Many common matrix/vector operations have � In each level: been standardized into an API called the BLAS multithreading, prefetching, vectorization, (Basic Linear Algebra Subroutines) loop unrolling, etc. • Level 1: vector operations � In level 2, especially in level 3: blocking copy, scale, dot, add, norms, … • Operate on sub-blocks of the matrix that fit the • Level 2: matrix-vector operations multiply, triangular solve, … memory architecture well • Level 3: matrix-matrix operations � General goal: multiply, triangular solve, … if it � s easy to phrase an operation in terms � FORTRAN bias, but callable from other langs of BLAS, get speed+safety for free � Goals: • The higher the level better • As fast as possible, but still safe/accurate � www.netlib.org/blas cs542g-term1-2006 5 cs542g-term1-2006 6

LAPACK Specializations � The BLAS only solves triangular systems � When solving a linear system, first question to • Forward or backward substitution ask: what sort of system? � LAPACK is a higher level API for matrix � Many properties to consider: operations: • Single precision or double? • Solving linear systems • Real or complex? • Solving linear least squares problems • Invertible or (nearly) singular? • Solving eigenvalue problems • Symmetric/Hermitian? • Definite or Indefinite? � Built on the BLAS, with blocking in mind to keep • Dense or sparse or specially structured? high performance • Multiple right-hand sides? � Biggest advantage: safety � LAPACK/BLAS take advantage of many of these • Designed to handle difficult problems gracefully (sparse matrices the big exception…) � www.netlib.org/lapack cs542g-term1-2006 7 cs542g-term1-2006 8 Accuracy Condition Number � Before jumping into algorithms, how � Sometimes we can estimate the condition accurate can we hope to be in solving a number of a matrix a priori linear system? � Special case: for a symmetric matrix, � Key idea: backward error analysis 2-norm condition number is ratio of extreme eigenvalues � Assume calculated answer is the exact solution of a perturbed problem . � LAPACK also provides cheap estimates • Try to construct a vector ||x|| that comes close to maximizing ||A -1 x|| cs542g-term1-2006 9 cs542g-term1-2006 10 Gaussian Elimination LU Factorization � Each step of row reduction is multiplication � Let � s start with the simplest unspecialized by an elementary matrix algorithm: Gaussian Elimination � Gathering these together, we find GE is � Assume the matrix is invertible, but essentially a matrix factorization: otherwise nothing special known about it A=LU � GE simply is row-reduction to upper where triangular form, followed by backwards L is lower triangular (and unit diagonal), U is upper triangular substitution • Permuting rows if we run into a zero � Solving Ax=b by GE is then Ly=b Ux=y cs542g-term1-2006 11 cs542g-term1-2006 12

Block Approach to LU � Rather than get bogged down in details of GE (hard to see forest for trees) � Partition the equation A=LU � Gives natural formulas for algorithms � Extends to block algorithms cs542g-term1-2006 13