### Homework 9.

Assigned: Mon Nov 8. Due: Mon Nov 15, at midnight

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:

• if the user requests only one output argument (nargout <= 1), the function returns the Cholesky factor (in the same diagonal format) if A is numerically positive definite (the factorization does not break down) and generates an error return (type "help error") if breakdown occurs
• if the user requests two output arguments (nargout = 2), the function returns the Cholesky factor (in the diagonal format) and sets the second output argument to zero if A is numerically positive definite, and, in the case of breakdown, returns a partial Cholesky factor (exactly what depends on the implementation) sand sets the second argument to a positive integer indicating where breakdown occurred.
Include enough comments so that the grader can understand what you are doing. Test your program on randomly generated examples (both positive definite and indefinite) as well as the discrete Laplacian matrix coded in discr_lapl.m. This returns a matrix in the standard full format and will need to be converted to your diagonal input format before being passed to your code.

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.