diff --git a/CMakeLists.txt b/CMakeLists.txt
index 689dc3dbad31bc1155646ae34097d909259cb7cf..97b46d60e7648d1b1e767a667df4c6f44c7006b4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -97,6 +97,7 @@ SET(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
include/eigenpy/decompositions/LDLT.hpp
include/eigenpy/decompositions/LLT.hpp
include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
+ include/eigenpy/decompositions/minres.hpp
)
SET(${PROJECT_NAME}_HEADERS
diff --git a/include/eigenpy/decompositions/minres.hpp b/include/eigenpy/decompositions/minres.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..45123ce2e188208357b84a246d6584f45400f1ee
--- /dev/null
+++ b/include/eigenpy/decompositions/minres.hpp
@@ -0,0 +1,158 @@
+/*
+ * Copyright 2021 INRIA
+ */
+
+#ifndef __eigenpy_decomposition_minres_hpp__
+#define __eigenpy_decomposition_minres_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+
+#include <iostream>
+#include <Eigen/Core>
+#include <unsupported/Eigen/IterativeSolvers>
+
+#include "eigenpy/utils/scalar-name.hpp"
+
+namespace eigenpy
+{
+ template<typename _Solver>
+ struct IterativeSolverBaseVisitor
+ : public boost::python::def_visitor< IterativeSolverBaseVisitor<_Solver> >
+ {
+ typedef _Solver Solver;
+ typedef typename Solver::MatrixType MatrixType;
+ typedef typename Solver::Preconditioner Preconditioner;
+ typedef typename Solver::Scalar Scalar;
+ typedef typename Solver::RealScalar RealScalar;
+
+ typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic,MatrixType::Options> MatrixXs;
+
+ template<class PyClass>
+ void visit(PyClass& cl) const
+ {
+ namespace bp = boost::python;
+ cl
+ .def("analyzePattern",
+ (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix))&Solver::analyzePattern,
+ bp::args("self","A"),
+ "Initializes the iterative solver for the sparsity pattern of the matrix A for further solving Ax=b problems.",
+ bp::return_self<>())
+ .def("factorize",
+ (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix))&Solver::factorize,
+ bp::args("self","A"),
+ "Initializes the iterative solver with the numerical values of the matrix A for further solving Ax=b problems.",
+ bp::return_self<>())
+ .def("compute",
+ (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix))&Solver::compute,
+ bp::args("self","A"),
+ "Initializes the iterative solver with the matrix A for further solving Ax=b problems.",
+ bp::return_self<>())
+
+ .def("rows",&Solver::rows,bp::arg("self"),"Returns the number of rows.")
+ .def("cols",&Solver::cols,bp::arg("self"),"Returns the number of columns.")
+ .def("tolerance",&Solver::tolerance,bp::arg("self"),
+ "Returns the tolerance threshold used by the stopping criteria.")
+ .def("setTolerance",&Solver::setTolerance,bp::args("self","tolerance"),
+ "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|.\n"
+ "The default value is the machine precision given by NumTraits<Scalar>::epsilon().",
+ bp::return_self<>())
+ .def("preconditioner",(Preconditioner & (Solver::*)())&Solver::preconditioner,bp::arg("self"),
+ "Returns a read-write reference to the preconditioner for custom configuration.",
+ bp::return_internal_reference<>())
+
+ .def("maxIterations",&Solver::maxIterations,bp::arg("self"),
+ "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",&Solver::setMaxIterations,bp::args("self","max_iterations"),
+ "Sets the max number of iterations.\n"
+ "Default is twice the number of columns of the matrix.",
+ bp::return_self<>())
+
+ .def("iterations",&Solver::iterations,bp::arg("self"),
+ "Returns the number of iterations performed during the last solve.")
+ .def("error",&Solver::error,bp::arg("self"),
+ "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",&Solver::error,bp::arg("info"),
+ "Returns Success if the iterations converged, and NoConvergence otherwise.")
+
+ .def("solveWithGuess",&solveWithGuess<MatrixXs,MatrixXs>,bp::args("self","b","x0"),
+ "Returns the solution x of A x = b using the current decomposition of A and x0 as an initial solution.")
+
+ .def("solve",&solve<MatrixXs>,bp::args("self","b"),
+ "Returns the solution x of A x = b using the current decomposition of A where b is a right hand side matrix or vector.")
+ ;
+ }
+
+ private:
+
+ template<typename MatrixOrVector1, typename MatrixOrVector2>
+ static MatrixOrVector1 solveWithGuess(const Solver & self, const MatrixOrVector1 & b, const MatrixOrVector2 & guess)
+ {
+ return self.solveWithGuess(b,guess);
+ }
+
+ template<typename MatrixOrVector>
+ static MatrixOrVector solve(const Solver & self, const MatrixOrVector & mat_or_vec)
+ {
+ MatrixOrVector res = self.solve(mat_or_vec);
+ return res;
+ }
+ };
+
+ template<typename _MatrixType>
+ struct MINRESSolverVisitor
+ : public boost::python::def_visitor< MINRESSolverVisitor<_MatrixType> >
+ {
+
+ typedef _MatrixType MatrixType;
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1,MatrixType::Options> VectorXs;
+ typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic,MatrixType::Options> MatrixXs;
+ typedef Eigen::MINRES<MatrixType> Solver;
+
+ template<class PyClass>
+ void visit(PyClass& cl) const
+ {
+ namespace bp = boost::python;
+ cl
+ .def(bp::init<>("Default constructor"))
+ .def(bp::init<MatrixType>(bp::arg("matrix"),
+ "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()."))
+
+ .def(IterativeSolverBaseVisitor<Solver>())
+ ;
+ }
+
+ static void expose()
+ {
+ static const std::string classname = "MINRES" + scalar_name<Scalar>::shortname();
+ expose(classname);
+ }
+
+ static void expose(const std::string & name)
+ {
+ namespace bp = boost::python;
+ bp::class_<Solver, boost::noncopyable>(name.c_str(),
+ "A minimal residual solver for sparse symmetric problems.\n"
+ "This class allows to solve for A.x = b sparse linear problems using the MINRES algorithm of Paige and Saunders (1975). The sparse matrix A must be symmetric (possibly indefinite). The vectors x and b can be either dense or sparse.\n"
+ "The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations and NumTraits<Scalar>::epsilon() for the tolerance.\n",
+ bp::no_init)
+ .def(MINRESSolverVisitor());
+ }
+
+ private:
+
+ template<typename MatrixOrVector>
+ static MatrixOrVector solve(const Solver & self, const MatrixOrVector & vec)
+ {
+ return self.solve(vec);
+ }
+ };
+
+} // namespace eigenpy
+
+#endif // ifndef __eigenpy_decomposition_minres_hpp__
diff --git a/src/decompositions/decompositions.cpp b/src/decompositions/decompositions.cpp
index a3bbd36e8427767e321bba9b71d1c6c4db8cd990..93ed609d6a215620a5f7bff37b2a2b9f871f2714 100644
--- a/src/decompositions/decompositions.cpp
+++ b/src/decompositions/decompositions.cpp
@@ -1,5 +1,5 @@
/*
- * Copyright 2020 INRIA
+ * Copyright 2020-2021 INRIA
*/
#include "eigenpy/fwd.hpp"
@@ -10,6 +10,7 @@
#include "eigenpy/decompositions/SelfAdjointEigenSolver.hpp"
#include "eigenpy/decompositions/LLT.hpp"
#include "eigenpy/decompositions/LDLT.hpp"
+#include "eigenpy/decompositions/minres.hpp"
namespace eigenpy
{
@@ -22,6 +23,8 @@ namespace eigenpy
SelfAdjointEigenSolverVisitor<MatrixXd>::expose("SelfAdjointEigenSolver");
LLTSolverVisitor<MatrixXd>::expose("LLT");
LDLTSolverVisitor<MatrixXd>::expose("LDLT");
+
+ MINRESSolverVisitor<MatrixXd>::expose("MINRES");
{
bp::enum_<DecompositionOptions>("DecompositionOptions")
diff --git a/unittest/CMakeLists.txt b/unittest/CMakeLists.txt
index ba18c9c60dc2beac50415faba1c8a34422e4d22c..d8730c43cbfe7fb3959434f7820fc0247e67d8ae 100644
--- a/unittest/CMakeLists.txt
+++ b/unittest/CMakeLists.txt
@@ -1,6 +1,6 @@
#
# Copyright (c) 2014-2019 CNRS
-# Copyright (c) 2018-2020 INRIA
+# Copyright (c) 2018-2021 INRIA
#
MACRO(ADD_LIB_UNIT_TEST test)
@@ -62,3 +62,6 @@ SET_TESTS_PROPERTIES("py-LLT" PROPERTIES DEPENDS ${PYWRAP})
ADD_PYTHON_UNIT_TEST("py-LDLT" "unittest/python/test_LDLT.py" "python/eigenpy;unittest")
SET_TESTS_PROPERTIES("py-LDLT" PROPERTIES DEPENDS ${PYWRAP})
+
+ADD_PYTHON_UNIT_TEST("py-MINRES" "unittest/python/test_MINRES.py" "python/eigenpy;unittest")
+SET_TESTS_PROPERTIES("py-MINRES" PROPERTIES DEPENDS ${PYWRAP})
diff --git a/unittest/python/test_MINRES.py b/unittest/python/test_MINRES.py
new file mode 100644
index 0000000000000000000000000000000000000000..67cbfac96376ed56e3841dc99868467326318956
--- /dev/null
+++ b/unittest/python/test_MINRES.py
@@ -0,0 +1,14 @@
+import eigenpy
+
+import numpy as np
+import numpy.linalg as la
+
+dim = 100
+A = np.eye(dim)
+
+minres = eigenpy.MINRES(A)
+
+X = np.random.rand(dim,20)
+B = A.dot(X)
+X_est = minres.solve(B)
+assert eigenpy.is_approx(A.dot(X_est),B,1e-6)