Write a MATLAB function to compute the Cholesky factorization of a symmetric positive definite band matrix A. Instead of working with Algorithm 23.1 in the text (which is a little hard to understand), you can modify my code cholDirect.m, presented in class. The program should work for bandwidth 2p+1 for any p from 0 to m-1, where m is the order of the matrix, and it should be efficient. The input is not A itself but an m by p+1 matrix, say Adiags, which stores the main diagonal and the upper diagonals of A. You do not need to explicitly pass the band parameter p since it is available by looking at size(Adiags). Write the code so that it optionally returns one or two output arguments, just like the built-in chol:
For the discrete Laplacian, are the zeros inside the band preserved in the Cholesky factor? Why or why not? (You will see better ways to solve the discrete Laplacian next semester.) Make timing comparisons of your routine with the built-in chol, both passing the Laplacian matrix as a full (dense) matrix , and also passing it as a sparse matrix (just type A = sparse(A) or A = full(A) to make the conversions). Make sure the matrix dimension is large enough that the timing comparisons are meaningful.
What is the flop count for your implementation, in terms of m and p? Verify this experimentally, by seeing what happens to the CPU time if you double m or p respectively. Answering this question carefully is an important part of the assignment.