diff --git a/CMakeLists.txt b/CMakeLists.txt index c66588eb756cd6a13e8b8749e1e4b96703114176..1b57e649e6e0c3a11552e238774b5e91aef0232e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -82,10 +82,13 @@ INCLUDE_DIRECTORIES(${NUMPY_INCLUDE_DIRS}) # ---------------------------------------------------- SET(${PROJECT_NAME}_SOLVERS_HEADERS solvers/solvers.hpp + solvers/preconditioners.hpp solvers/IterativeSolverBase.hpp solvers/LeastSquaresConjugateGradient.hpp solvers/ConjugateGradient.hpp solvers/SparseSolverBase.hpp + solvers/BasicPreconditioners.hpp + solvers/BFGSPreconditioners.hpp ) SET(HEADERS @@ -139,6 +142,7 @@ ENDFOREACH(header) # --- TARGETS ---------------------------------------- # ---------------------------------------------------- SET(${PROJECT_NAME}_SOLVERS_SOURCES + src/solvers/preconditioners.cpp src/solvers/solvers.cpp ) SET(${PROJECT_NAME}_SOURCES diff --git a/python/main.cpp b/python/main.cpp index 891cacc0fb8537fbb068a35545ef4d0ae281ea52..2f57c4e444240fb22357f45ec4f453990e5948f1 100644 --- a/python/main.cpp +++ b/python/main.cpp @@ -18,6 +18,7 @@ #include "eigenpy/eigenpy.hpp" #include "eigenpy/geometry.hpp" #include "eigenpy/solvers/solvers.hpp" +#include "eigenpy/solvers/preconditioners.hpp" #include <iostream> @@ -28,5 +29,6 @@ BOOST_PYTHON_MODULE(eigenpy) eigenpy::exposeAngleAxis(); eigenpy::exposeQuaternion(); eigenpy::exposeSolvers(); + eigenpy::exposePreconditioners(); } diff --git a/src/solvers/BFGSPreconditioners.hpp b/src/solvers/BFGSPreconditioners.hpp new file mode 100644 index 0000000000000000000000000000000000000000..e3cbb22a1dc2a6d2a24ca34f718375cefb3a0cd0 --- /dev/null +++ b/src/solvers/BFGSPreconditioners.hpp @@ -0,0 +1,100 @@ +/* + * 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_bfgs_preconditioners_hpp__ +#define __eigenpy_bfgs_preconditioners_hpp__ + +#include <boost/python.hpp> +#include <Eigen/Core> +#include <Eigen/IterativeLinearSolvers> + +#include "eigenpy/solvers/BasicPreconditioners.hpp" + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename Preconditioner> + struct BFGSPreconditionerBaseVisitor + : public bp::def_visitor< BFGSPreconditionerBaseVisitor<Preconditioner> > + { + + typedef Eigen::VectorXd VectorType; + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(PreconditionerBaseVisitor<Preconditioner>()) + .def("rows",&Preconditioner::rows,"Returns the number of rows in the preconditioner.") + .def("cols",&Preconditioner::cols,"Returns the number of cols in the preconditioner.") + .def("dim",&Preconditioner::dim,"Returns the dimension of the BFGS preconditioner") + .def("update",(const Preconditioner& (Preconditioner::*)(const VectorType &, const VectorType &) const)&Preconditioner::update,bp::args("s","y"),"Update the BFGS estimate of the matrix A.",bp::return_value_policy<bp::reference_existing_object>()) + .def("reset",&Preconditioner::reset,"Reset the BFGS estimate.") + ; + + } + + static void expose(const std::string & name) + { + bp::class_<Preconditioner>(name, + bp::no_init) + .def(BFGSPreconditionerBaseVisitor<Preconditioner>()) + ; + + } + + + + }; + + template<typename Preconditioner> + struct LimitedBFGSPreconditionerBaseVisitor + : public bp::def_visitor< LimitedBFGSPreconditionerBaseVisitor<Preconditioner> > + { + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(PreconditionerBaseVisitor<Preconditioner>()) + .def(BFGSPreconditionerBaseVisitor<Preconditioner>()) + .def("resize",&Preconditioner::resize,bp::arg("dim"),"Resizes the preconditionner with size dim.",bp::return_value_policy<bp::reference_existing_object>()) + ; + + } + + static void expose(const std::string & name) + { + bp::class_<Preconditioner>(name.c_str(), + bp::no_init) + .def(LimitedBFGSPreconditionerBaseVisitor<Preconditioner>()) + ; + + } + + + + }; + + + + + + + +} // namespace eigenpy + +#endif // ifndef __eigenpy_bfgs_preconditioners_hpp__ diff --git a/src/solvers/BasicPreconditioners.hpp b/src/solvers/BasicPreconditioners.hpp new file mode 100644 index 0000000000000000000000000000000000000000..b9a085ef5628fc2273ae96353f23d90a53edd073 --- /dev/null +++ b/src/solvers/BasicPreconditioners.hpp @@ -0,0 +1,138 @@ +/* + * 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_basic_preconditioners_hpp__ +#define __eigenpy_basic_preconditioners_hpp__ + +#include <boost/python.hpp> +#include <Eigen/Core> +#include <Eigen/IterativeLinearSolvers> + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename Preconditioner> + struct PreconditionerBaseVisitor + : public bp::def_visitor< PreconditionerBaseVisitor<Preconditioner> > + { + typedef Eigen::MatrixXd MatrixType; + typedef Eigen::VectorXd VectorType; + + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(bp::init<>("Default constructor")) + .def(bp::init<MatrixType>(bp::arg("A"),"Initialize the preconditioner with matrix A for further Az=b solving.")) + .def("info",&Preconditioner::info, + "Returns success if the Preconditioner has been well initialized.") + .def("solve",&solve,bp::arg("b"), + "Returns the solution A * z = b where the preconditioner is an estimate of A^-1.") + + .def("compute",&Preconditioner::template compute<MatrixType>,bp::arg("mat"), + "Initialize the preconditioner from the matrix value.", + bp::return_value_policy<bp::reference_existing_object>()) + ; + + } + + private: + + static VectorType solve(Preconditioner & self, const VectorType & b) + { + return self.solve(b); + } + + }; + + template<typename Scalar> + struct DiagonalPreconditionerVisitor : PreconditionerBaseVisitor<Eigen::DiagonalPreconditioner<Scalar> > + { + typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> MatrixType; + typedef Eigen::DiagonalPreconditioner<Scalar> Preconditioner; + + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(PreconditionerBaseVisitor<Preconditioner>()) + .def("rows",&Preconditioner::rows,"Returns the number of rows in the preconditioner.") + .def("cols",&Preconditioner::cols,"Returns the number of cols in the preconditioner.") + ; + + } + + static void expose() + { + bp::class_<Preconditioner>("DiagonalPreconditioner", + "A preconditioner based on the digonal entrie.\n" + "This class allows to approximately solve for A.x = b problems assuming A is a diagonal matrix.", + bp::no_init) + ; + + } + }; + + template<typename Scalar> + struct LeastSquareDiagonalPreconditionerVisitor : PreconditionerBaseVisitor<Eigen::LeastSquareDiagonalPreconditioner<Scalar> > + { + typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> MatrixType; + typedef Eigen::LeastSquareDiagonalPreconditioner<Scalar> Preconditioner; + + template<class PyClass> + void visit(PyClass&) const + { + } + + static void expose() + { + bp::class_<Preconditioner>("LeastSquareDiagonalPreconditioner", + "Jacobi preconditioner for LeastSquaresConjugateGradient.\n" + "his class allows to approximately solve for A' A x = A' b problems assuming A' A is a diagonal matrix.", + bp::no_init) + .def(DiagonalPreconditionerVisitor<Scalar>()) + ; + + } + }; + + struct IdentityPreconditionerVisitor : PreconditionerBaseVisitor<Eigen::IdentityPreconditioner > + { + typedef Eigen::IdentityPreconditioner Preconditioner; + + template<class PyClass> + void visit(PyClass&) const + { + } + + static void expose() + { + bp::class_<Preconditioner>("IdentityPreconditioner", + bp::no_init) + .def(PreconditionerBaseVisitor<Preconditioner>()) + ; + + } + }; + + + + +} // namespace eigenpy + +#endif // ifndef __eigenpy_basic_preconditioners_hpp__ diff --git a/src/solvers/IterativeSolverBase.hpp b/src/solvers/IterativeSolverBase.hpp index a341b082b979267a4c10c3689cd2ec15cd9e6a67..ffa827de77b89d28efdb7463dc750022beb9366c 100644 --- a/src/solvers/IterativeSolverBase.hpp +++ b/src/solvers/IterativeSolverBase.hpp @@ -33,6 +33,7 @@ namespace eigenpy { typedef typename IterativeSolver::MatrixType MatrixType; + typedef typename IterativeSolver::Preconditioner Preconditioner; typedef Eigen::VectorXd VectorType; template<class PyClass> @@ -68,6 +69,7 @@ namespace eigenpy 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.") + .def("preconditioner",(Preconditioner & (IS::*)(void))&IS::preconditioner,"Returns a read-write reference to the preconditioner for custom configuration.",bp::return_internal_reference<>()) ; } diff --git a/src/solvers/preconditioners.cpp b/src/solvers/preconditioners.cpp new file mode 100644 index 0000000000000000000000000000000000000000..e0d4501db9204987172ea7ab4ec0d0f9fcced9b8 --- /dev/null +++ b/src/solvers/preconditioners.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/preconditioners.hpp" +#include "eigenpy/solvers/BasicPreconditioners.hpp" +#include "eigenpy/solvers/BFGSPreconditioners.hpp" + +namespace eigenpy +{ + void exposePreconditioners() + { + using namespace Eigen; + + DiagonalPreconditionerVisitor<double>::expose(); + LeastSquareDiagonalPreconditionerVisitor<double>::expose(); + IdentityPreconditionerVisitor::expose(); + LimitedBFGSPreconditionerBaseVisitor< LimitedBFGSPreconditioner<double,Eigen::Dynamic,Eigen::Upper|Eigen::Lower> >::expose("LimitedBFGSPreconditioner"); + + + } +} // namespace eigenpy diff --git a/src/solvers/preconditioners.hpp b/src/solvers/preconditioners.hpp new file mode 100644 index 0000000000000000000000000000000000000000..111a9564058b2efced35fd333a2b7a029f041635 --- /dev/null +++ b/src/solvers/preconditioners.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_preconditioners_hpp__ +#define __eigenpy_preconditioners_hpp__ + +namespace eigenpy +{ + + void exposePreconditioners(); + +} // namespace eigenpy + +#endif // define __eigenpy_preconditioners_hpp__