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 <[email protected]> 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. > > > >
