From 824a52b3b1ca1925b1f7b045f6171fe97cb930d3 Mon Sep 17 00:00:00 2001 From: jcarpent <jcarpent@laas.fr> Date: Wed, 18 Oct 2017 10:24:06 +0200 Subject: [PATCH] [Solvers] Expose iterative solvers of Eigen --- CMakeLists.txt | 8 ++- python/main.cpp | 2 + src/solvers/ConjugateGradient.hpp | 61 ++++++++++++++++ src/solvers/IterativeSolverBase.hpp | 103 ++++++++++++++++++++++++++++ src/solvers/SparseSolverBase.hpp | 58 ++++++++++++++++ src/solvers/solvers.cpp | 34 +++++++++ src/solvers/solvers.hpp | 27 ++++++++ 7 files changed, 292 insertions(+), 1 deletion(-) create mode 100644 src/solvers/ConjugateGradient.hpp create mode 100644 src/solvers/IterativeSolverBase.hpp create mode 100644 src/solvers/SparseSolverBase.hpp create mode 100644 src/solvers/solvers.cpp create mode 100644 src/solvers/solvers.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 627d0dcc..adff098c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ # -# Copyright (c) 2015-2016 LAAS-CNRS +# Copyright (c) 2015-2017 LAAS-CNRS # # This file is part of eigenpy. # eigenpy is free software: you can redistribute it and/or @@ -90,8 +90,13 @@ SET(HEADERS registration.hpp angle-axis.hpp quaternion.hpp + solvers/solvers.hpp + solvers/IterativeSolverBase.hpp + solvers/ConjugateGradient.hpp + solvers/SparseSolverBase.hpp ) MAKE_DIRECTORY("${${PROJECT_NAME}_BINARY_DIR}/include/eigenpy") +MAKE_DIRECTORY("${${PROJECT_NAME}_BINARY_DIR}/include/eigenpy/solvers") INCLUDE_DIRECTORIES(${${PROJECT_NAME}_BINARY_DIR}/include/eigenpy) #FOREACH(header ${HEADERS}) @@ -131,6 +136,7 @@ SET(${PROJECT_NAME}_SOURCES src/eigenpy.cpp src/angle-axis.cpp src/quaternion.cpp + src/solvers/solvers.cpp ) ADD_LIBRARY(${PROJECT_NAME} SHARED ${${PROJECT_NAME}_SOURCES} ${${PROJECT_NAME}_HEADERS}) diff --git a/python/main.cpp b/python/main.cpp index 95ec5c44..891cacc0 100644 --- a/python/main.cpp +++ b/python/main.cpp @@ -17,6 +17,7 @@ #include "eigenpy/eigenpy.hpp" #include "eigenpy/geometry.hpp" +#include "eigenpy/solvers/solvers.hpp" #include <iostream> @@ -26,5 +27,6 @@ BOOST_PYTHON_MODULE(eigenpy) eigenpy::enableEigenPy(); eigenpy::exposeAngleAxis(); eigenpy::exposeQuaternion(); + eigenpy::exposeSolvers(); } diff --git a/src/solvers/ConjugateGradient.hpp b/src/solvers/ConjugateGradient.hpp new file mode 100644 index 00000000..710d9f14 --- /dev/null +++ b/src/solvers/ConjugateGradient.hpp @@ -0,0 +1,61 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_conjugate_gradient_hpp__ +#define __eigenpy_conjugate_gradient_hpp__ + +#include <boost/python.hpp> +#include <Eigen/IterativeLinearSolvers> + +#include "eigenpy/solvers/IterativeSolverBase.hpp" + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename ConjugateGradient> + struct ConjugateGradientVisitor + : public boost::python::def_visitor< ConjugateGradientVisitor<ConjugateGradient> > + { + typedef typename ConjugateGradient::MatrixType MatrixType; + + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(bp::init<>("Default constructor")) + .def(bp::init<MatrixType>(bp::arg("A"),"Initialize the solver with matrix A for further Ax=b solving.\n" + "This constructor is a shortcut for the default constructor followed by a call to compute().")) + ; + + } + + static void expose() + { + bp::class_<ConjugateGradient,boost::noncopyable>("ConjugateGradient", + bp::no_init) + .def(IterativeSolverVisitor<ConjugateGradient>()) + .def(ConjugateGradientVisitor<ConjugateGradient>()) + ; + + } + + }; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_conjugate_gradient_hpp__ diff --git a/src/solvers/IterativeSolverBase.hpp b/src/solvers/IterativeSolverBase.hpp new file mode 100644 index 00000000..a341b082 --- /dev/null +++ b/src/solvers/IterativeSolverBase.hpp @@ -0,0 +1,103 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_iterative_solver_base_hpp__ +#define __eigenpy_iterative_solver_base_hpp__ + +#include <boost/python.hpp> +#include <Eigen/Core> + +#include "eigenpy/solvers/SparseSolverBase.hpp" + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename IterativeSolver> + struct IterativeSolverVisitor + : public boost::python::def_visitor< IterativeSolverVisitor<IterativeSolver> > + { + + typedef typename IterativeSolver::MatrixType MatrixType; + typedef Eigen::VectorXd VectorType; + + template<class PyClass> + void visit(PyClass& cl) const + { + typedef IterativeSolver IS; + + cl + .def(SparseSolverVisitor<IS>()) + .def("error",&IS::error,"Returns the tolerance error reached during the last solve.\n" + "It is a close approximation of the true relative residual error |Ax-b|/|b|.") + .def("info",&IS::info,"Returns success if the iterations converged, and NoConvergence otherwise.") + .def("iterations",&IS::iterations,"Returns the number of iterations performed during the last solve.") + .def("maxIterations",&IS::maxIterations,"Returns the max number of iterations.\n" + "It is either the value setted by setMaxIterations or, by default, twice the number of columns of the matrix.") + .def("setMaxIterations",&IS::setMaxIterations,"Sets the max number of iterations.\n" + "Default is twice the number of columns of the matrix.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("tolerance",&IS::tolerance,"Returns he tolerance threshold used by the stopping criteria.") + .def("setTolerance",&IS::setTolerance,"Sets the tolerance threshold used by the stopping criteria.\n" + "This value is used as an upper bound to the relative residual error: |Ax-b|/|b|. The default value is the machine precision.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("analyzePattern",&analyzePattern,bp::arg("A"),"Initializes the iterative solver for the sparsity pattern of the matrix A for further solving Ax=b problems.\n" + "Currently, this function mostly calls analyzePattern on the preconditioner.\n" + "In the future we might, for instance, implement column reordering for faster matrix vector products.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("factorize",&factorize,bp::arg("A"),"Initializes the iterative solver with the numerical values of the matrix A for further solving Ax=b problems.\n" + "Currently, this function mostly calls factorize on the preconditioner.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("compute",&compute,bp::arg("A"),"Initializes the iterative solver with the numerical values of the matrix A for further solving Ax=b problems.\n" + "Currently, this function mostly calls factorize on the preconditioner.\n" + "In the future we might, for instance, implement column reordering for faster matrix vector products.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("solveWithGuess",&solveWithGuess,bp::args("b","x0"), + "Returns the solution x of Ax = b using the current decomposition of A and x0 as an initial solution.") + ; + + } + + private: + + static IterativeSolver & factorize(IterativeSolver & self, const MatrixType & m) + { + return self.factorize(m); + } + + static IterativeSolver & compute(IterativeSolver & self, const MatrixType & m) + { + return self.compute(m); + } + + static IterativeSolver & analyzePattern(IterativeSolver & self, const MatrixType & m) + { + return self.analyzePattern(m); + } + + static VectorType solveWithGuess(IterativeSolver & self, const Eigen::VectorXd & b, const Eigen::VectorXd & x0) + { + return self.solveWithGuess(b,x0); + } + + }; + + + +} // namespace eigenpy + +#endif // ifndef __eigenpy_iterative_solver_base_hpp__ diff --git a/src/solvers/SparseSolverBase.hpp b/src/solvers/SparseSolverBase.hpp new file mode 100644 index 00000000..e96a9741 --- /dev/null +++ b/src/solvers/SparseSolverBase.hpp @@ -0,0 +1,58 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_sparse_solver_base_hpp__ +#define __eigenpy_sparse_solver_base_hpp__ + +#include <boost/python.hpp> +#include <Eigen/Core> + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename SparseSolver> + struct SparseSolverVisitor + : public bp::def_visitor< SparseSolverVisitor<SparseSolver> > + { + + typedef Eigen::VectorXd VectorType; + + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def("solve",&solve,bp::arg("b"), + "Returns the solution x of Ax = b using the current decomposition of A.") + ; + + } + + private: + + static VectorType solve(SparseSolver & self, const VectorType & b) + { + return self.solve(b); + } + + }; + + + +} // namespace eigenpy + +#endif // ifndef __eigenpy_sparse_solver_base_hpp__ diff --git a/src/solvers/solvers.cpp b/src/solvers/solvers.cpp new file mode 100644 index 00000000..ba594d0e --- /dev/null +++ b/src/solvers/solvers.cpp @@ -0,0 +1,34 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#include "eigenpy/solvers/solvers.hpp" +#include "eigenpy/solvers/ConjugateGradient.hpp" + +namespace eigenpy +{ + void exposeSolvers() + { + using namespace Eigen; + ConjugateGradientVisitor< ConjugateGradient<MatrixXd,Lower|Upper> >::expose(); + + boost::python::enum_<Eigen::ComputationInfo>("ComputationInfo") + .value("Success",Eigen::Success) + .value("NumericalIssue",Eigen::NumericalIssue) + .value("NoConvergence",Eigen::NoConvergence) + .value("InvalidInput",Eigen::InvalidInput) + ; + } +} // namespace eigenpy diff --git a/src/solvers/solvers.hpp b/src/solvers/solvers.hpp new file mode 100644 index 00000000..ff4bc7f1 --- /dev/null +++ b/src/solvers/solvers.hpp @@ -0,0 +1,27 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_solvers_hpp__ +#define __eigenpy_solvers_hpp__ + +namespace eigenpy +{ + + void exposeSolvers(); + +} // namespace eigenpy + +#endif // define __eigenpy_solvers_hpp__ -- GitLab