Using Amgcl with Eigen::SparseMatrix

[SomayehHesabi]

Dear Denis, I came accross your open-source library while looking for a fastest solver for my sparse system Ax=b. In fact, I already used Eigen Solvers for my case, however, as they are a bit slow, I would like to apply your method. Actually, I am a bit confused. If I pass my Eigen::SparseMatrix directly to my Amgcl solver, do I have to implement all the stuff regrading sparse matrices, such as row_iterator, etc? Or, probably I am wrong and there is another simple way to use Eigen::SparseMatrix methods.

I really appriceite it if you could help me. Somayeh

[ddemidov]

You should be able to #include <amgcl/backend/eigen.hpp> and use amgcl::backend::eigen<double> instead of amgcl::backend::builtin<double>.

[ddemidov]

Here is an example: https://gist.github.com/ddemidov/ac9a060507edea7586d0516b174ccd13.

~EDIT: note that the example is rather old and boost::tie should be changed to std::tie in line 50.~

EDIT 2: I've just updated the example to be compatible with the latest amgcl version.

[ddemidov]

The latest commit (41d479cdc021a3453c72c622e089f4027057a775) splits eigen backend into adapter and backend parts, so that eigen types can be used with builtin backend, see example here:

https://gist.github.com/ddemidov/9265647f0e18aeead3bbd3a072c10d18

The important differences with the example above are (the example above still works without any changes):

• Include <amgcl/adapter/eigen.hpp> instead of amgcl/backend/eigen.hpp, and include <amgcl/backend/builtin.hpp> (lines 8-9),
• Call AMGCL_USE_EIGEN_VECTORS_WITH_BUILTIN_BACKEND macro to tell amgcl that it is ok to use Eigen vectors as RHS and solution with builtin backend (line 17)
• Replace amgcl::backend::eigen<double> with amgcl::backend::builtin<double> (lines 42-49).

The advantages are:

• Faster setup time: no need to copy internal matrices into Eigen format. Actually, in case of eigen backend the improvement is negligible, since MappedSparseMatrix is used and no memory is actually copied.
• Potentially faster solve time: builtin backend uses OpenMP and Eigen does not.

[xgarnaud]

Dear Denis,

Thank you for your library and for the example above. When I try to use it to solve with a matrix different from that used to build the preconditioner, I get the following error

~/Documents/Perso/minifem/./external/amgcl/amgcl/backend/interface.hpp: In instantiation of 'struct amgcl::backend::residual_impl<Eigen::SparseMatrix<double, 1>, Eigen::Matrix<double, -1, 1>, Eigen::Matrix<double, -1, 1>, amgcl::backend::numa_vector<double>, void>':
~/Documents/Perso/minifem/./external/amgcl/amgcl/backend/interface.hpp:315:60:   required from 'void amgcl::backend::residual(const Vector1&, const Matrix&, const Vector2&, Vector3&) [with Matrix = Eigen::SparseMatrix<double, 1>; Vector1 = Eigen::Matrix<double, -1, 1>; Vector2 = Eigen::Matrix<double, -1, 1>; Vector3 = amgcl::backend::numa_vector<double>]'
~/Documents/Perso/minifem/./external/amgcl/amgcl/solver/cg.hpp:156:30:   required from 'std::tuple<long unsigned int, typename amgcl::math::scalar_of<typename Backend::value_type, void>::type> amgcl::solver::cg<Backend, InnerProduct>::operator()(const Matrix&, const Precond&, const Vec1&, Vec2&&) const [with Matrix = Eigen::SparseMatrix<double, 1>; Precond = amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>; Vec1 = Eigen::Matrix<double, -1, 1>; Vec2 = Eigen::Matrix<double, -1, 1>&; Backend = amgcl::backend::builtin<double>; InnerProduct = amgcl::solver::detail::default_inner_product; typename amgcl::math::scalar_of<typename Backend::value_type>::type = double]'
~/Documents/Perso/minifem/./external/amgcl/amgcl/make_solver.hpp:132:34:   required from 'std::tuple<long unsigned int, typename amgcl::math::scalar_of<typename IterativeSolver::backend_type::value_type, void>::type> amgcl::make_solver<Precond, IterativeSolver>::operator()(const Matrix&, const Vec1&, Vec2&&) const [with Matrix = Eigen::SparseMatrix<double, 1>; Vec1 = Eigen::Matrix<double, -1, 1>; Vec2 = Eigen::Matrix<double, -1, 1>&; Precond = amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>; IterativeSolver = amgcl::solver::cg<amgcl::backend::builtin<double> >; typename amgcl::math::scalar_of<typename IterativeSolver::backend_type::value_type>::type = double]'
~/Documents/Perso/minifem/src/fem/linear_solvers/./amgcl_linear_solvers.h:42:22:   required from 'std::tuple<int, double> AMGCLLinearSolver<Solver>::solve(const VectorXd&, Eigen::VectorXd&) const [with Solver = amgcl::make_solver<amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>, amgcl::solver::cg<amgcl::backend::builtin<double> > >; Eigen::VectorXd = Eigen::Matrix<double, -1, 1>]'
~/Documents/Perso/minifem/src/fem/linear_solvers/amgcl_linear_solvers.cc:22:16:   required from here
~/Documents/Perso/minifem/./external/amgcl/amgcl/backend/interface.hpp:188:55: error: no type named 'RESIDUAL_NOT_IMPLEMENTED' in class Eigen::SparseMatrix<double, 1>'
  188 |     typedef typename Matrix::RESIDUAL_NOT_IMPLEMENTED type;
      |                                                       ^~~~
~/Documents/Perso/minifem/./external/amgcl/amgcl/backend/interface.hpp: In instantiation of 'void amgcl::backend::residual(const Vector1&, const Matrix&, const Vector2&, Vector3&) [with Matrix = Eigen::SparseMatrix<double, 1>; Vector1 = Eigen::Matrix<double, -1, 1>; Vector2 = Eigen::Matrix<double, -1, 1>; Vector3 = amgcl::backend::numa_vector<double>]':
~/Documents/Perso/minifem/./external/amgcl/amgcl/solver/cg.hpp:156:30:   required from 'std::tuple<long unsigned int, typename amgcl::math::scalar_of<typename Backend::value_type, void>::type> amgcl::solver::cg<Backend, InnerProduct>::operator()(const Matrix&, const Precond&, const Vec1&, Vec2&&) const [with Matrix = Eigen::SparseMatrix<double, 1>; Precond = amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>; Vec1 = Eigen::Matrix<double, -1, 1>; Vec2 = Eigen::Matrix<double, -1, 1>&; Backend = amgcl::backend::builtin<double>; InnerProduct = amgcl::solver::detail::default_inner_product; typename amgcl::math::scalar_of<typename Backend::value_type>::type = double]'
~/Documents/Perso/minifem/./external/amgcl/amgcl/make_solver.hpp:132:34:   required from 'std::tuple<long unsigned int, typename amgcl::math::scalar_of<typename IterativeSolver::backend_type::value_type, void>::type> amgcl::make_solver<Precond, IterativeSolver>::operator()(const Matrix&, const Vec1&, Vec2&&) const [with Matrix = Eigen::SparseMatrix<double, 1>; Vec1 = Eigen::Matrix<double, -1, 1>; Vec2 = Eigen::Matrix<double, -1, 1>&; Precond = amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>; IterativeSolver = amgcl::solver::cg<amgcl::backend::builtin<double> >; typename amgcl::math::scalar_of<typename IterativeSolver::backend_type::value_type>::type = double]'
~/Documents/Perso/minifem/src/fem/linear_solvers/./amgcl_linear_solvers.h:42:22:   required from 'std::tuple<int, double> AMGCLLinearSolver<Solver>::solve(const VectorXd&, Eigen::VectorXd&) const [with Solver = amgcl::make_solver<amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>, amgcl::solver::cg<amgcl::backend::builtin<double> > >; Eigen::VectorXd = Eigen::Matrix<double, -1, 1>]'
~/Documents/Perso/minifem/src/fem/linear_solvers/amgcl_linear_solvers.cc:22:16:   required from here
~/Documents/Perso/minifem/./external/amgcl/amgcl/backend/interface.hpp:315:60: error: 'apply' is not a member of 'amgcl::backend::residual_impl<Eigen::SparseMatrix<double, 1>, Eigen::Matrix<double, -1, 1>, Eigen::Matrix<double, -1, 1>, amgcl::backend::numa_vector<double>, void>'
  315 |     residual_impl<Matrix, Vector1, Vector2, Vector3>::apply(rhs, A, x, r);
      |     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~
~/Documents/Perso/minifem/./external/amgcl/amgcl/backend/interface.hpp: In instantiation of 'struct amgcl::backend::spmv_impl<double, Eigen::SparseMatrix<double, 1>, amgcl::backend::numa_vector<double>, double, amgcl::backend::numa_vector<double>, void>':
~/Documents/Perso/minifem/./external/amgcl/amgcl/backend/interface.hpp:303:60:   required from 'void amgcl::backend::spmv(Alpha, const Matrix&, const Vector1&, Beta, Vector2&) [with Alpha = double; Matrix = Eigen::SparseMatrix<double, 1>; Vector1 = amgcl::backend::numa_vector<double>; Beta = double; Vector2 = amgcl::backend::numa_vector<double>]'
~/Documents/Perso/minifem/./external/amgcl/amgcl/solver/cg.hpp:171:30:   required from 'std::tuple<long unsigned int, typename amgcl::math::scalar_of<typename Backend::value_type, void>::type> amgcl::solver::cg<Backend, InnerProduct>::operator()(const Matrix&, const Precond&, const Vec1&, Vec2&&) const [with Matrix = Eigen::SparseMatrix<double, 1>; Precond = amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>; Vec1 = Eigen::Matrix<double, -1, 1>; Vec2 = Eigen::Matrix<double, -1, 1>&; Backend = amgcl::backend::builtin<double>; InnerProduct = amgcl::solver::detail::default_inner_product; typename amgcl::math::scalar_of<typename Backend::value_type>::type = double]'
~/Documents/Perso/minifem/./external/amgcl/amgcl/make_solver.hpp:132:34:   required from 'std::tuple<long unsigned int, typename amgcl::math::scalar_of<typename IterativeSolver::backend_type::value_type, void>::type> amgcl::make_solver<Precond, IterativeSolver>::operator()(const Matrix&, const Vec1&, Vec2&&) const [with Matrix = Eigen::SparseMatrix<double, 1>; Vec1 = Eigen::Matrix<double, -1, 1>; Vec2 = Eigen::Matrix<double, -1, 1>&; Precond = amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>; IterativeSolver = amgcl::solver::cg<amgcl::backend::builtin<double> >; typename amgcl::math::scalar_of<typename IterativeSolver::backend_type::value_type>::type = double]'
~/Documents/Perso/minifem/src/fem/linear_solvers/./amgcl_linear_solvers.h:42:22:   required from 'std::tuple<int, double> AMGCLLinearSolver<Solver>::solve(const VectorXd&, Eigen::VectorXd&) const [with Solver = amgcl::make_solver<amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>, amgcl::solver::cg<amgcl::backend::builtin<double> > >; Eigen::VectorXd = Eigen::Matrix<double, -1, 1>]'
~/Documents/Perso/minifem/src/fem/linear_solvers/amgcl_linear_solvers.cc:22:16:   required from here
~/Documents/Perso/minifem/./external/amgcl/amgcl/backend/interface.hpp:181:51: error: no type named 'SPMV_NOT_IMPLEMENTED' in 'class Eigen::SparseMatrix<double, 1>'
  181 |     typedef typename Matrix::SPMV_NOT_IMPLEMENTED type;
      |                                                   ^~~~
~/Documents/Perso/minifem/./external/amgcl/amgcl/backend/interface.hpp: In instantiation of 'void amgcl::backend::spmv(Alpha, const Matrix&, const Vector1&, Beta, Vector2&) [with Alpha = double; Matrix = Eigen::SparseMatrix<double, 1>; Vector1 = amgcl::backend::numa_vector<double>; Beta = double; Vector2 = amgcl::backend::numa_vector<double>]':
~/Documents/Perso/minifem/./external/amgcl/amgcl/solver/cg.hpp:171:30:   required from 'std::tuple<long unsigned int, typename amgcl::math::scalar_of<typename Backend::value_type, void>::type> amgcl::solver::cg<Backend, InnerProduct>::operator()(const Matrix&, const Precond&, const Vec1&, Vec2&&) const [with Matrix = Eigen::SparseMatrix<double, 1>; Precond = amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>; Vec1 = Eigen::Matrix<double, -1, 1>; Vec2 = Eigen::Matrix<double, -1, 1>&; Backend = amgcl::backend::builtin<double>; InnerProduct = amgcl::solver::detail::default_inner_product; typename amgcl::math::scalar_of<typename Backend::value_type>::type = double]'
~/Documents/Perso/minifem/./external/amgcl/amgcl/make_solver.hpp:132:34:   required from 'std::tuple<long unsigned int, typename amgcl::math::scalar_of<typename IterativeSolver::backend_type::value_type, void>::type> amgcl::make_solver<Precond, IterativeSolver>::operator()(const Matrix&, const Vec1&, Vec2&&) const [with Matrix = Eigen::SparseMatrix<double, 1>; Vec1 = Eigen::Matrix<double, -1, 1>; Vec2 = Eigen::Matrix<double, -1, 1>&; Precond = amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>; IterativeSolver = amgcl::solver::cg<amgcl::backend::builtin<double> >; typename amgcl::math::scalar_of<typename IterativeSolver::backend_type::value_type>::type = double]'
~/Documents/Perso/minifem/src/fem/linear_solvers/./amgcl_linear_solvers.h:42:22:   required from 'std::tuple<int, double> AMGCLLinearSolver<Solver>::solve(const VectorXd&, Eigen::VectorXd&) const [with Solver = amgcl::make_solver<amgcl::amg<amgcl::backend::builtin<double>, amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0>, amgcl::solver::cg<amgcl::backend::builtin<double> > >; Eigen::VectorXd = Eigen::Matrix<double, -1, 1>]'
~/Documents/Perso/minifem/src/fem/linear_solvers/amgcl_linear_solvers.cc:22:16:   required from here
~/Documents/Perso/minifem/./external/amgcl/amgcl/backend/interface.hpp:303:60: error: 'apply' is not a member of 'amgcl::backend::spmv_impl<double, Eigen::SparseMatrix<double, 1>, amgcl::backend::numa_vector<double>, double, amgcl::backend::numa_vector<double>, void>'
  303 |     spmv_impl<Alpha, Matrix, Vector1, Beta, Vector2>::apply(alpha, A, x, beta, y);
      |     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~
In file included from ~local/boost/1.72/include/boost/serialization/strong_typedef.hpp:27,
                 from ~local/boost/1.72/include/boost/serialization/serialization.hpp:14,
                 from ~/local/boost/1.72/include/boost/serialization/split_free.hpp:22,
                 from ~/local/boost/1.72/include/boost/serialization/nvp.hpp:35,
                 from ~/local/boost/1.72/include/boost/multi_index/detail/index_loader.hpp:22,
                 from ~/local/boost/1.72/include/boost/multi_index/detail/index_base.hpp:33,
                 from ~/local/boost/1.72/include/boost/multi_index/detail/base_type.hpp:21,
                 from ~/local/boost/1.72/include/boost/multi_index_container.hpp:36,
                 from ~/local/boost/1.72/include/boost/property_tree/ptree.hpp:21,
                 from ~/Documents/Perso/minifem/src/fem/linear_solvers/././linear_solver.h:4,
                 from ~/Documents/Perso/minifem/src/fem/linear_solvers/./amgcl_linear_solvers.h:2,
                 from ~/Documents/Perso/minifem/src/fem/linear_solvers/amgcl_linear_solvers.cc:1:
~/local/boost/1.72/include/boost/operators.hpp:162:18: warning: inline function 'bool boost::operators_impl::operator!=(const boost::multi_index::detail::bidir_node_iterator<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::null_augment_policy, boost::multi_index::detail::index_node_base<std::pair<const std::__cxx11::basic_string<char>, boost::property_tree::basic_ptree<std::__cxx11::basic_string<char>, std::__cxx11::basic_string<char> > >, std::allocator<std::pair<const std::__cxx11::basic_string<char>, boost::property_tree::basic_ptree<std::__cxx11::basic_string<char>, std::__cxx11::basic_string<char> > > > > > >&, const boost::multi_index::detail::bidir_node_iterator<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::null_augment_policy, boost::multi_index::detail::index_node_base<std::pair<const std::__cxx11::basic_string<char>, boost::property_tree::basic_ptree<std::__cxx11::basic_string<char>, std::__cxx11::basic_string<char> > >, std::allocator<std::pair<const std::__cxx11::basic_string<char>, boost::property_tree::basic_ptree<std::__cxx11::basic_string<char>, std::__cxx11::basic_string<char> > > > > > >&)' used but never defined
  162 |      friend bool operator!=(const T& x, const T& y) { return !static_cast<bool>(x == y); }
      |                  ^~~~~~~~
make[2]: *** [src/fem/CMakeFiles/minifem_fem.dir/linear_solvers/amgcl_linear_solvers.cc.o] Error 1
make[1]: *** [src/fem/CMakeFiles/minifem_fem.dir/all] Error 2
make: *** [all] Error 2

Do you know how I can use this function? Best regards

[ddemidov]

Would this workaround work for you? (Using crs_tuple adapter with pointers to eigen data)

[ddemidov]

Or you could just use eigen backend instead of the builtin one: https://gist.github.com/cc1d78bc82b21cb9d09f7426ea0e5e56

[xgarnaud]

Perfect! Thank you very much