Lecture 23 Linear Algebra

1
Lecture 23 Linear Algebra

• Comp 208 Computers in Engineering
• Yi Lin
• Fall, 2005

2
Lecture 23 Learning goals

• Linear algebra in C
• Vector dot and cross products
• Matrix operations
• Addition, subtraction, multiplication, transpose
• Solving matrix equations
• Ax b matrix-vector equation
• Useful for systems of linear equations

3
Vectors in C

• As you know neither FORTRAN nor C have built-in
  knowledge of vectors
• In Cs case this is partly because vectors are
  not important when making an operating system
• We will use arrays as vectors and write
  functions to perform vector operations

4
Addition and subtraction of vectors

• Very simple, just add every corresponding cell
• add(u,v,w,N)
• for(i 0i lt Ni)
• wi ui vi
• sub(u,v,w,N)
• for(i 0i lt Ni)
• wi ui - vi

5
Vector dot product

• Dot product is well known and easy to implement
  in any number of dimensions
• Let N be the number of dimensions for u ? v
• dotproduct(u,v,N)
• dot 0
• for(i 0i lt Ni)
• dot viui
• return dot

6
Obtaining the magnitude of a vector

• As done very often in engineering math use the
  dot product
• modulo(u,N)
• return sqrt(dot(u,u,N))

7
Vector cross product

• Cross product is well known for vectors in 3
  dimensions.
• Definition can get very hairy for more than 3
  dimensions good thing that we almost never want
  to take cross products of vectors in more than 3
  dimensions
• Simple to code, w u X v as
• cross(u,v,w)
• w0 u1v2 v1u2
• w1 v0u2 u0v2
• w2 u0v1 u1v0

8
Matrices in C

• Much like vectors, we must simulate matrices with
  multi-dimensional arrays and implement all the
  relevant functionality
• Determinants, adjoints and inversions are
  extremely difficult to do so we will not attempt
  them directly here

9
Matrix addition and subtraction

• Much like the vector versions but we must add
  and subtract multiple columns per row instead of
  one
• add(m1,m2,m3,n,m)
• for(i 0i lt ni)
• for(j 0j lt mj)
• m3ij m1ij m2ij

10
Matrix multiplication

• More difficult, but applying the definition
  directly just plain works
• matmult(m1,m2,m3,n,m,l)
• for(i 0i lt ni)
• for(j 0j lt lj)
• m3ij 0
• for(k 0k lt mk)
• m3ij m1ik m2kj

11
A note about dot products

• In math we actually treat all vectors as matrices
  which have one dimension set to 1
• Hence the duality between row vectors and column
  vectors
• If we were to use matmult on vectors, we would
  get the dot product out BUT one of them would
  have to be row vector and the other one a column
  vector

12
Matrix transposition

• We can save ourselves the pain of having to
  write routines for adjoints, determinants and
  then inverses if we know that a matrix is
  orthogonal
• Transposition is straightforward
• transpose(m1,m2,n,m)
• for(i 0i lt n i)
• for(j ij lt mj)
• m2ji m1ij

13
Other operations

• We would like to be able to do matrix inversions
  because engineers spend most of their time
  building problems like this

14
Is there another way?

• Instead of using matrix inversions we could use
  Gaussian Elimination. First lets talk about
  Gaussian Elimination by itself.
• Here the definition again just works.
• Lets write it in code

15
What GE does
16
Usage of GE

• Matrix must be square for a solution to exist (N
  equations N unknowns)
• GENP only triangularizes the matrix,
  back-substitution is need to fully solve the
  problem
• Only good if the matrix is strictly
  diagonally-dominant

17
GE Not Pivoting

• // it will transform the matrix a to a upper
  right matrix.
• void GENP(double a33, double b, int n)
• int r, row, col
• double scale
• int i, j
• for(r0 rltn-1 r)
• for(rowr1 row lt n row)
• // We start from the second row (i.e.,
  row1) since 1st row doesn't need to be changed
• scale arowr / arr
• for(colr col lt n col)
• // Each element (row, col) must substract
  a value in order to make (row,0) to (row, row-1)
  as 0
• arowcol arowcol - arcol
  scale
• brow brow - scale br

18
GENP

• Gaussian Elimination with No Pivoting
• This is the most basic gaussian elimination you
  learn in math courses
• Mathematically sound but what happens if there is
  a 0 on the diagonal?
• What happens if there is a very small number on
  the diagonal?
• GENP works, in theory, but for systems that
  arent strictly diagonally-dominant, it is
  neither optimal or guaranteeing an accurate
  solution

19
How can we make this better?

• Use pivoting!
• At i-th step, find the row that has the largest
  element in the i-th column
• Swap row i with that row
• Continue with elimination
• Guarantees that we always divide by the largest
  number available (stabilizes the machine results)
• Prevents dividing anything other than 0 by 0.

20
Usage of GEPP

• Gaussian Elimination with Partial Pivoting
• Partial because we only pivot with rows
• Full pivoting is harder, uses both columns and
  rows
• Can be used virtually any time that GENP blows up
• Even better than GENP on problems where GENP
  works
• Use GEPP whenever you can

21
GE Partial Pivoting

• // it will transform the matrix a to a upper
  right matrix.
• void GEPP(double a33, double b, int n)
• int r, row, col
• double scale
• int i, tmpR
• double biggestRR
• for(r0 rltn-1 r)
• // Before doing any upper trianglirization,
  swap the row whichever has the biggest
  element(r,r) with the row r.
• biggestRR arr
• for(ir1 iltn i)
• if(air gt biggestRR)
• tmpR i
• biggestRR air

• if(tmpR!r)
• // swap if the row with biggest element
  (r,r) is not the row r
• for(i0 iltn i)
• biggestRR ari
• ari atmpRi
• atmpRi biggestRR
• // The rests are the same as GENP
• return