Implement Point/Block Jacobi Preconditioning

[sslattery]

Per #25 block Jacobi preconditioning is the way to go for higher order SPN equations. This will now need to be implemented. First, we must choose either a Denovo or MCLS implementation followed by the actual implementation itself.

This task is done when the implementation is unit tested and utilized for a simple higher order SPN problem.

[sslattery]

Working on this now. I'm adding general full-matrix preconditioning capabilities in MCLS in a way that makes sense for Monte Carlo (where we need the entire preconditioner). This fits right in with the general left/right preconditioned Neumann-Ulam solver and MCSA/Sequential MC solvers I already have implemented. For now, only point and block Jacobi methods will be implemented.

The end result will be another set of Thyra/Stratimikos interfaces, this time for preconditioners. That way, this is still plug and play for users while providing us an extensible framework to add more preconditioners if necessary.

[sslattery]

The point/block Jacobi preconditioners have been implemented in MCLS for both Tpetra and Epetra under a general full-matrix preconditioner interface. These are unit tested and running in parallel.

The next step is to implement the Thyra preconditioner interfaces and write the Stratimikos adapter. The general preconditioner interface in MCLS mimics that of Thyra, so this will be really easy. I'll add these new preconditioners to the already existing Thyra and Stratimikos unit tests as well. Maybe the block preconditioner will get some more of those test matrices to converge.

[sslattery]

The Thyra interfaces have been implemented and tested for the new preconditioners. In addition, the Stratimikos adapter is written and also tested. Next step is to add the stratimikos code to Exnihilo and add preconditioning to the tests in my Exnihilo branch.

[sslattery]

The preconditioners have been incorporated into Exnihilo through Stratimikos and the kba_solvers tests are passing in parallel. I will run an SPN problem next to verify the preconditioning.

[sslattery]

I've run a couple of SPN problems with the new preconditioners through the pyspn interface and loading an xml file with MCLS solver parameters. SP1 is converging with point jacobi while SP3 is converging with block Jacobi. There are some issues with convergence for SP5 and SP7 but I think this is implementation related.

[sslattery]

I was incorrectly setting the pivot points generated by the LAPACK LU factorization in the matrix inversion code for inverting the Jacobi blocks. I have convergence now for the adjoint method and a simple Richardson iteration using the preconditioned iteration matrix and multigroup SP5 and SP7 problems.

I'm observing some issues with MCSA and this leads me to believe its related to how I'm generating the preconditioned system. I'll double check the solution from the two methods that are working as this will then rule out the iteration matrix being formed and the Monte Carlo solution.

[sslattery]

I've just verified that the adjoint Neumann-Ulam solver with block preconditioning is getting the same answer to machine precision as the GMRES solver in Aztec for a P7 problem. I'll check the Richardson iteration next.

[sslattery]

The Richardson iteration with the block preconditioned iteration matrix also gets the same answer as Aztec GMRES to machine precision. Must be in the MCSA implementation....

[sslattery]

So I discovered that my Exnihilo branch was not rebased to the master since early March, meaning that Steven's coefficient fixes in the SPN equations were not present, making a significant change in the final matrix generated and likely the erratic convergence behavior I have been observing. I'm having some rebase issues with my branch but hopefully I'll get this worked out soon to verify this hypothesis and close this issue.

[sslattery]

Just did the rebase and rebuilt with the correct SPN coefficients. Everything is converging as expected...and its fast :) We'll see just how fast in the next stage of this research.

I'm giving this to Paul to review.