Lab assignment L1

LU factorization

For due dates see the calendar web page. For hand in details the lab main page.

Objective

To practically implement the routines involved in LU factorization and thereby get a hands-on knowledge of the numerics of factorization, how to write numerical routines in Matlab, particularly how to access sub-parts of vectors and matrices, and how to write clean, correct vercorized code.

Procedure

The goal is to compose a linear equation solver from the following parts:
  1. Write a function [M,L] = elimMat(A,k) which given matrix A, index k computes the elimination matrix M and L as specified on page 67
  2. Write an myLU function which combines M1'M2' ... M2 M1 A into L and U
  3. Do a backSubst function. (a fwdSubst function is provided below as an example)
  4. Show how to use your programs to solve Ax=b by method of ( p. 68 )
    1. Factorization: L U = A
    2. fwdSubst to solve for y in: Ly=b
    3. backSubst to solve for x in: Ux = y
  5. Test your program with e.g. example 2.13 and one test of your own design. Write a m-file script to execute the text and comment on what's happening in each step.
It is not required to implement pivoting in this lab, but think about: What happens in your program if diagonal element 0? If the system is ill conditioned?

(NOT required) In particular, for those who want some challenge, think about the following problem: Given a matrix A and vector b, decide if the linear system Ax=b is consistent or not; if yes, then deliver one solution (an arbitrary one if there are multiple ones).

Test your program on the matrix, say,

A = [1, 2, 3, 4, 0; 0, 0, 0, 2, 3; 0, 0, 0, 2, 1; 0, 0, 0, 0, 0; 0, 0, 0, 0, 0];

b = [1; 2; 3; 0; 0];

Style

Show how to use matlab matrix operations to write the above in a clear compact form with a minimum number of for loops. (e.g. copying the for-loop based algorithms 2.2 and 2.3 will not yield good marks). Write your code in a readable, commented and suitably indented fashion. (There is a matlab mode for emacs you can download)

Hints

First understand and experiment with the supplied algorithm on the forward substitution part alone. Modify it to do back substitution. Then write the elimMat and myLU routines. myLU can be written recursively much the same way as the fwdSubst, or written with one "for" loop. You need no loop in computing the elimination matrix and at most one loop (or recursion) in computing the LU and forward, backward substitutions. Achieving this depends mostly on skills in linear algebra to do operations on vectors and matrices as a whole instead of element by element and on understanding Matlab's matrix elements adressing as described in the first part of the tutorial. This way each function (elim matrix, lu and substitutions) can each be written in just 5-10 lines of code. 
% fwdSubst  Forward substitution to solve Ly=b by tail recursion

% cmput 340, 2004

function y = fwdSubst(L,b,k)

  [m,n]=size(L);
  if ~exist('k')  % If first call no k param given, but k=1
    k=1;
  end
  
  y=b(k)/L(k,k)
  if k < n % Recursion step
    l = [zeros(k,1);L(k+1:m,k)]
    y = [y;fwdSubst(L,b-y*l,k+1)]
  end

Misc. Q and A: 


Newsgroups: ualberta.courses.cmput.340
Subject: Re: L1 myLu(...)
References: 
--text follows this line--
Paul Cartledge  writes:

> If I could get some clarification on the myLU(...) function, that would be
> great.
> 
> I believe it is wrong to assume that we will only be getting 3x3 matrices
> for this lab, so I am wondering how we are to send in parameters to
> myLU().

Obvioulsy your function should not just handle 3x3 matrices.
Look back to the reading list for lab 1, and make sure you
have read and understood the parts of the Matlab tutorial that
was covered. In particular, ch 3 has the answer for how to extract
properties (such as size) from a passed in matrix. Especially pay
attention to table 4. Then take a look at the sample fwdSubst function
supplied and see how parameters are used.

> 
> If there's a varying number of matrices to be sent in, is there a

Why would you like to pass in a varying number of matrices?

> way to store the intermediate matrices we need for computation obtained
> from elimMat() so that we can send them into myLU().

If you use recursion you have the benefit of a new set of intermediate
variables at each call. Again study the tutorial, then the sample routine.

/Martin

----------------------------------------

Newsgroups: ualberta.courses.cmput.340
Subject: Re: Leading Diag Entries of Unreduced Portions being Zero
References: 
--text follows this line--
Don Rivers  writes:

> The textbook talks about the possibility of there being a leading zero in
> a part we've yet to reduce.  Is there an elegant way for us to interchange
> rows, 

1. Yes, this is essentially the same as piviting

2. You are not required to implement this in the lab, but in preparation
for the lab test you should be aware of the issue from your readings 
(as you obvioulsy are) so you can answer oral questions.


> or will all of the matrices that we will be tested on not have
> leading zeros in parts that we are working on?

On a test matrix with a leading 0 you will have to point out why 
your code doesn't handle it, identify the problem, and the particular
place in the code.

/Martin

------------------------------------

Newsgroups: ualberta.courses.cmput.340
Subject: Re: NxN, or MxN?
References: 
--text follows this line--
Don Rivers  writes:

> I don't know if this has been answered, but are only going to be tested on
> Square Matrices, or on nonSquares as well?

Not required to handle non-square matrices.

> 
> I notice that matlab's lu() complains if you pass it a nonsquare matrix,
> so do we have a similar requirement?

Actually, it is possible to compute an L U factorization through
elimination even for an mxn matrix. It is just not implemented
in matlab 5 (but was fixed in matlab 6). 

/Martin