Re: [eigen] Reusing preconditioners

[Rasmus Munk Larsen <rmlarsen@xxxxxxxxxx>]

Hi Fernando,

I think this is a very useful feature to have, as you say for FEM it's quite common. Would you mind submitting this information as a feature request at: https://gitlab.com/libeigen/eigen/-/issues?sort=created_date&state=opened
I'm afraid we don't have a lot of cycles at the moment, so a merge request with the necessary changes would be even better :-)

Best,

   Rasmus

On Tue, Mar 8, 2022 at 12:13 PM <fpp2005st@xxxxxxxxxxx> wrote:

Dear all,

While upgrading from Eigen 3.2.x to 3.4.0 I faced a issue in my code,
related with the reuse of preconditioners. Say for example that you have
initially a given sparse matrix A and the corresponding iterative solver
and preconditioner, that you use to solve for a given rhs vector b.
Then, the non-zero coefficients of matrix A get slightly modified, such
that I would like to keep reusing the same preconditioner to solve for
the same or different vector b. With Eigen 3.2.x, I could directly call
for solve() with the second matrix (without a previous call to
compute(), as this would create a new preconditioner) and it would work.
With 3.4.0, it only works if the second matrix is not reset with a call
to setFromTriplets(), for example if the second matrix results from the
first one multiplied by a factor.

 From what I have read in the documentation, there are several warnings
stating that the solver creates a link to matrix A and if it changes,
then one needs to compute() or factorize() before a new solve(), i.e. it
shall not be possible to reuse the preconditioner. However, that was
allowed in 3.2.x and it keeps working in 3.4 as long as there is no call
to setFromTriplets(). That was a quite handy feature, as it allowed
reusing the preconditioner (factorization) in similar matrices, i.e.
sparse matrices with the same filling pattern, but different
coefficients. This happens a lot in finite element analysis.

My question is: is it possible to reuse a preconditioner in Eigen 3.4
after calling setFromTriplets() in the matrix originally used to build
the preconditioner? I think that this should be allowed as long as we
ensure that the non-zero patterns keeps unchanged. Please find bellow a
snippet showing what I described:

//****

#include <Eigen/IterativeLinearSolvers>
#include <vector>
#include <iostream>

typedef Eigen::SparseMatrix<double> SpMat;
typedef Eigen::Triplet<double> T;

void popTri
(
   std::vector<T>& tripletList,
   double f
)
{
   tripletList.push_back(T(0,0,-4*f));
   tripletList.push_back(T(0,1,2*f));

   tripletList.push_back(T(1,0,1*f));
   tripletList.push_back(T(1,1,-4*f));
   tripletList.push_back(T(1,2,1*f));

   tripletList.push_back(T(2,1,1*f));
   tripletList.push_back(T(2,2,-4*f));
}

int main()
{
   SpMat A(3,3);
   std::vector<T> tripletList;
   Eigen::Vector3d b1(1,2,3);
   Eigen::BiCGSTAB<SpMat, Eigen::IncompleteLUT<double> > solver;


   // Solve 1 (ok in both 3.2 and 3.4)
   popTri(tripletList,1);
   A.setFromTriplets(tripletList.begin(), tripletList.end());
   solver.analyzePattern(A);
   solver.factorize(A);
   std::cout << "A1\n" << solver.solve(b1) << std::endl;

   // Solve 2 (ok in both 3.2 and 3.4)
   A*=2.;
   std::cout << "A2\n" << solver.solve(b1) << std::endl;

   // Solve 3 (ok in 3.2 but not in 3.4)
   // Note: solve 3 makes the same modification in A as solve 2,
   // with the only difference that it uses setFromTriplets instead of
   // directly multiplying by the factor
   tripletList.clear();
   popTri(tripletList,2);
   A.setFromTriplets(tripletList.begin(), tripletList.end());
   std::cout << "A3\n" << solver.solve(b1) << std::endl;

   return 0;
}

//***

Thank you very much for your time.

Best regards,

Fernando.