[eigen] Reusing preconditioners

[fpp2005st@xxxxxxxxxxx]

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.