Exercise "matrix diagonalization"

1. (6 points) Jacobi diagonalization with cyclic sweeps

   1. Implement the Jacobi eigenvalue method for real symmetric matrices using the cyclic sweep algorithm: sequentially eliminate all off-diagonal elements row after row.

   2. Make some interesting examples proving your implementation works as intended.

   Typically one stores the diagonal elements of the matrix in a separate vector and only updates the strict upper triangular part of the matrix. The user can then always restore their symmetric matrix from the remaining lower triangle.

2. (3 points) Jacobi diagonalization with classic sweeps and indexing of largest elements in rows

   Implement the classic Jacobi diagonalization algorithm with indexing:

   • Before starting diagonalization find out the largest off-diagonal elements in each row and store their indices in a separate index array. It's O(n²) operations—expensive—but you only do it only once in the beginning.
   • Find the largest off-diagonal element from the pre-found largest elements in each row (only O(n) operations) and zero it by a Jacobi rotation.
   • Update your index. Since the Jacobi rotation only affects elements in rows and columns with indices p and q the update is also only O(n) operations.
   • Repeat the last two steps until converged

   Compare the speed of the classical-sweep-with-indexing and the cyclic-sweep algorithms on, say, random (real symmetric) matrices of different sizes. If you do everything right the classic must do fewer rotations, but the total speeds must be close.

3. (1 points) Consider a sparse real symmetric matrix (of modest size) which is diagonal except for the first column (and first row). Such matrices appear often in quantum-mechanical variational calculations. Find out whether diagonalization of this sort of matrix by Jacobi method is generally faster than diagonalization of a random (symmetric) matrix.
4. (0 point) Comparison with library routines

   Compare the speed of your implementations with library implementations (e.g. GSL, Armadillo, and NumPy).