QR + Block Arnoldi

[Simkern]

Here comes the updated block versions of the algorithms.

QR algorithm

• qr_factorization: Fixed unnecessary step and added comments about the second pass.
• test_qr_breakdown: A special test case added to the QR test suite to check QR stability; Highlight when the single-pass GS-based QR algorithm fails.

Block Arnoldi

• block_arnoldi_factorization: Block version of the regular arnoldi factorization. Very condensed through the use of Krylov basis utilities and the QR factorization.
• test_block_arnoldi_full_factorization and test_block_arnoldi_basis_orthogonality are added to the arnoldi testsuite. They are essentially copies of the regular arnoldi tests adapted for non-unity blocks.

For now I have not touched the regular arnoldi algorithms. If we want we can either write an interface to choose the appropriate subroutine or replace the arnoldi entirely by the block version. The block size is an optional parameter for the subroutines.

[Simkern]

All the requests are implemented.

• initalize_kylov_subspace now takes an optional argument with the initial vector. Default is all zeros.
• arnoldi_factorization has been completely replaced with the block version. I have kept the separate tests for unit-size and larger blocks, it might come in handy for debugging later.
• All the examples have been refactored and streamlined, in particular regarding the initialisation.

I was thinking to replace all the orthonormalisation loops with the calls of mat_mult so that they would easily generalize to their block version but it turns out that this is quite messy unless the algorithms are actually written with a (possibly unit-valued) block size. This is mainly due to the fact that mat_mult is defined for matrix/Krylov basis-valued inputs and outputs and one would have to artificially include the extra dimension even if it is singleton.

There are 3 alternatives:

1. We simply leave the code as it is.
2. We fully generalize the algorithms to their block versions (this might require a few minor changes)
3. We generalize mat_mult with a specific procedure for rank 1 matrices, i.e. connect it to get_vec and the regular Kylov dot product.