I am trying to use a custom preconditioner for an iterative solver (CG for instance) with Eigen. Specifically, I have to solve a similar problem multiple times: the matrix changes slightly but stays close to a mean matrix. I would like to compute a Cholesky decomposition of my mean matrix and then use this as a preconditioner.
What I had in mind is something like:
ConjugateGradient< SparseMatrix<double>, Lower, CholmodSupernodalLLT<SparseMatrix<double>> > solver(meanMatrix);
solver.preconditioner().compute(meanMatrix);
// Loop on n similar matrices
for(int i = 0; i < n; i++){
// create matrix: it is similar (in structure and in values) to meanMatrix
SparseMatrix<double> matrix = ...;
// create right-hand-side
VectorXd rhs = ...;
// update matrix reference for solver
solver.compute(matrix);
// solve using the preconditioned CG
solver.solve(rhs);
}
The problem is that calling solver.compute(matrix) actually causes ConjugateGradient (in fact IterativeSolverBase) to call compute on its preconditioner (see l. 111 of IterativeSolverBase.h, Eigen 3.2.9):
m_preconditioner.compute(*mp_matrix);
In other words, the preconditioner based on the mean matrix is replaced by the Cholesky decomposition of the new matrix, and thus the CG solve converges in 1 iteration. On the contrary, I would like to keep the same preconditioner (the Cholesky decomposition of the mean matrix, computed once and for all before the loop), and solve for the different matrices using the preconditioned CG.
Is there an easy way to achieve what I am trying to do?
Many thanks in advance for your help! I hope this makes sense. If not, please do not hesitate to ask me to clarify.
