MA/CS 371 - Lab 7 (Two Weeks)
Gaussian Elimination with Scaled Partial Pivoting, LU Factorization, and Solving Linear Systems of Equations


Section 1: Introduction

In this lab you will write subroutines to

This week there are no supplied files; you will be writing everything from scratch. You should make a directory for this lab and work in it:

	mkdir ~/cs371/lab7
	cd ~/cs371/lab7
and when you are ready to submit files make sure they are in your lab7 directory.

Note that this lab is tougher than the previous ones have been; that's why you're getting two weeks to do it. I will be available to answer questions about this lab during my office hours and will take e-mail questions and random visits throughout the week.

IMPORTANT: It is OK for you to work in small groups on this lab, but if you do so, you must let me know who you're working with in order to get credit. Put your partners' names in a comment header in your code or e-mail me the names.

Section 2: Topics for Lecture

MATLAB does not have a function to compute Gaussian elimination, but it does have functions which directly compute the other functions in this lab: The lu function calculates the LU factorization of a matrix A with scaled partial pivoting. You can either get the raw answer with
	[L,U] = lu(A)
or you can see the permutation matrix (a representation of the pivoting) with
	[L,U,P] = lu(A)
The permutation matrix is another way to represent the information in the pivot vector. One way to interpret the pivot vector is to consider the ith element l(i) and say "the row that is now the ith element of the matrix A used to be the l(i)th element." So for example if you have the pivot vector
	l = [ 4 1 3 2 5 ]
then you know that the row that used to be the fourth row of A is now the first row, the row that used to be the first row is now the second, and the row that used to be the second row is now the fourth.

The permutation matrix P has the same information as the pivot vector l. P(i,j) == 1 exactly when l(i) = j. In other words, the above pivot vector l corresponds to the permutation matrix

		| 0 0 0 1 0 |
		| 1 0 0 0 0 |
	  P  =  | 0 0 1 0 0 |
		| 0 1 0 0 0 |
		| 0 0 0 0 1 |
The advantage of using the permutation matrix is that you can swap the rows of A by premultiplying the permutation matrix -- in other words the matrix product PA is A with the rows swapped around. You should try a few examples just to be sure you understand what lu is doing. Try it also with a bad matrix (for example, one with two repeated rows or with a row of zeros).

The backslash function calculates the answer to a system of equations. When you use your LU factorization to compute the answer to a system of equations, you can use backslash to check your answer. If you have a matrix A and a column vector b and you want to solve the system of equations Ax = b, you can do

	x = A\b
where that slash is the 'backslash' and not the regular forward slash.

Section 3: Gaussian Elimination with Scaled Partial Pivoting

Implement Gaussian Elimination with Scaled Partial Pivoting using the pseudocode on pp.236-237 as a template. Remember that you can often vectorize code in MATLAB to effectively perform an entire loop in one step when that loop involves doing the same operation to every element of a vector.

Your subroutine should differ from the pseudocode in that it should not overwrite the input matrix. Your first line should be

	function [U,l] = gausselim(A);
where U is the output matrix and l is the pivot vector. (It is traditional, but not required, MATLAB usage for matrices to be capitalized and vectors to be in lower case.)

Hint: One way to avoid overwriting A is to create the matrix U right away -- maybe make it a matrix full of zeros the same size as A by doing

	U = zeros(size(A));
and then build U up as the algorithm progresses.

If your code runs and gets a correct answer for square input matrices, you will get most of the credit for this problem. For full credit your program must reject non-square input matrices and contain at least one fewer for loop than the pseudocode on pp236-237.

Section 4: LU Factorization

Section 6.4 in the text includes a discussion of how to augment Gaussian Elimination code to compute the LU factorization of a matrix. Extend your gausselim function to compute the LU factorization of a matrix with scaled partial pivoting. Your new file will be called lufact.m and the first line should read
	function [L,U,l] = lufact(A);
If your code runs and gets a correct answer for square input matrices, you will get most of the credit for this problem. For full credit your program must reject non-square input matrices.

Section 5: Solving Linear Systems of Equations with LU

Use your LU routine to solve linear systems of equations. Given an input matrix A and a vector b, solve the system of equations Ax = b by computing the LU factorization of A, using forward substitution to solve Lz = b, and back substitution to solve Ux = z. Formulas for this method are given on p.266 of the text. The first line of your program should read
	function x = solve(L,U,l,b)
where L, U, and l are the output of your lufact program and b is the right hand side of the original linear system Ax = b. As before, you will get most of the credit for a program which runs correctly on square input matrices, but for full credit you will need to reject non-square inputs. The output vector x should be permuted so that it corresponds to the original order of the rows in A -- that's the purpose of passing the pivot vector l into the function.

To test your functions, start with simple examples and work up to more complicated ones. For example, try doing the Gaussian elimination of the following matrices to start with:

	| 1 0 0 |	| 1 0 0 |
	| 0 1 0 |	| 0 1 1 |
	| 0 0 1 |	| 0 1 2 |
and work up to more complicated examples. Checking your answers for lufact and solve will be easier because you will be able to check them against MATLAB's lu and backslash functions, respectively. I can suggest some good test matrices in lab if you wish.

Section 6: Submitting Your Work

When you are ready, run
	~cs371/submit lab7
to submit your gausselim.m, lufact.m, and solve.m programs. Please be sure to use those exact names or else submit won't find them.

These programs will be due by 11:59pm on Thursday, March 14. Please budget your time appropriately, especially if you don't feel like you got halfway there by the end of lab on the 7th.

There are no questions for the lab period on March 7th -- I will write up some questions for the 14th and link them to this page a few days before then.