### 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.