Lecture 23 Linear Algebra

1. Lecture 23 Linear Algebra Comp 208 Computers in Engineering Yi Lin Fall, 2005 Comp208 Computers in Engineering

2. Lecture – 23 Learning goals • Linear algebra in C • Vector dot and cross products • Matrix operations • Addition, subtraction, multiplication, transpose • Solving matrix equations • [A]x = b matrix-vector equation • Useful for systems of linear equations Comp208 Computers in Engineering

3. Vectors in C As you know neither FORTRAN nor C have built-in knowledge of vectors • In C’s 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 Comp208 Computers in Engineering

4. Addition and subtraction of vectors Very simple, just add every corresponding cell add(u,v,w,N) for(i = 0;i < N;++i) w[i] = u[i] + v[i] sub(u,v,w,N) for(i = 0;i < N;++i) w[i] = u[i] - v[i] Comp208 Computers in Engineering

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 = 0;i < N;++i) dot += v[i]*u[i]; return dot; Comp208 Computers in Engineering

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)); Comp208 Computers in Engineering

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) w[0] = u[1]*v[2] – v[1]*u[2]; w[1] = v[0]*u[2] – u[0]*v[2]; w[2] = u[0]*v[1] – u[1]*v[0]; Comp208 Computers in Engineering

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 Comp208 Computers in Engineering

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 = 0;i < n;++i) for(j = 0;j < m;++j) m3[i][j] = m1[i][j] + m2[i][j] Comp208 Computers in Engineering

10. Matrix multiplication More difficult, but applying the definition directly just plain works matmult(m1,m2,m3,n,m,l) for(i = 0;i < n;++i) for(j = 0;j < l;++j) m3[i][j] = 0; for(k = 0;k < m;k++) m3[i][j] += m1[i][k] * m2[k][j]; Comp208 Computers in Engineering

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 Comp208 Computers in Engineering

12. 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 = 0;i < n ;++i) for(j = i;j < m;++j) m2[j][i] = m1[i][j]; Matrix transposition Comp208 Computers in Engineering

13. Other operations We would like to be able to do matrix inversions because engineers spend most of their time building problems like this: Comp208 Computers in Engineering

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… Comp208 Computers in Engineering

15. What GE does Comp208 Computers in Engineering

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 Comp208 Computers in Engineering

17. GE Not Pivoting // it will transform the matrix a to a upper right matrix. void GENP(double a[3][3], double b[], int n){ int r, row, col; double scale; int i, j; for(r=0; r<n-1; r++){ for(row=r+1; row < n ; row++){ // We start from the second row (i.e., row=1) since 1st row doesn't need to be changed scale = a[row][r] / a[r][r]; for(col=r; col < n; col++){ // Each element (row, col) must substract a value in order to make (row,0) to (row, row-1) as 0 a[row][col] = a[row][col] - a[r][col] * scale; } b[row] = b[row] - scale * b[r]; } } return; } Comp208 Computers in Engineering

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 aren’t strictly diagonally-dominant, it is neither optimal or guaranteeing an accurate solution Comp208 Computers in Engineering

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. Comp208 Computers in Engineering

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… Comp208 Computers in Engineering

21. // it will transform the matrix a to a upper right matrix. void GEPP(double a[3][3], double b[], int n){ int r, row, col; double scale; int i, tmpR; double biggestRR; for(r=0; r<n-1; r++){ // Before doing any upper trianglirization, swap the row whichever has the biggest element(r,r) with the row r. biggestRR = a[r][r]; for(i=r+1; i<n; i++){ if(a[i][r] > biggestRR) { tmpR = i; biggestRR = a[i][r]; } } if(tmpR!=r){ // swap if the row with biggest element (r,r) is not the row r for(i=0; i<n; i++){ biggestRR = a[r][i]; a[r][i] = a[tmpR][i]; a[tmpR][i] = biggestRR; } } // The rests are the same as GENP return; } GE Partial Pivoting Comp208 Computers in Engineering