Skip to content

Commit 06e4150

Browse files
jcarpentjcarpent
authored
Merge pull request #246 from jcarpent/devel
Expose MINRES solver
2 parents 8dd520b + aad265b commit 06e4150

File tree

5 files changed

+181
-2
lines changed

5 files changed

+181
-2
lines changed

‎CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -97,6 +97,7 @@ SET(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
9797
include/eigenpy/decompositions/LDLT.hpp
9898
include/eigenpy/decompositions/LLT.hpp
9999
include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
100+
include/eigenpy/decompositions/minres.hpp
100101
)
101102

102103
SET(${PROJECT_NAME}_HEADERS

‎include/eigenpy/decompositions/minres.hpp‎

Lines changed: 158 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,158 @@
1+
/*
2+
* Copyright 2021 INRIA
3+
*/
4+
5+
#ifndef __eigenpy_decomposition_minres_hpp__
6+
#define __eigenpy_decomposition_minres_hpp__
7+
8+
#include "eigenpy/eigenpy.hpp"
9+
10+
#include <iostream>
11+
#include <Eigen/Core>
12+
#include <unsupported/Eigen/IterativeSolvers>
13+
14+
#include "eigenpy/utils/scalar-name.hpp"
15+
16+
namespace eigenpy
17+
{
18+
template<typename _Solver>
19+
struct IterativeSolverBaseVisitor
20+
: public boost::python::def_visitor< IterativeSolverBaseVisitor<_Solver> >
21+
{
22+
typedef _Solver Solver;
23+
typedef typename Solver::MatrixType MatrixType;
24+
typedef typename Solver::Preconditioner Preconditioner;
25+
typedef typename Solver::Scalar Scalar;
26+
typedef typename Solver::RealScalar RealScalar;
27+
28+
typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic,MatrixType::Options> MatrixXs;
29+
30+
template<class PyClass>
31+
void visit(PyClass& cl) const
32+
{
33+
namespace bp = boost::python;
34+
cl
35+
.def("analyzePattern",
36+
(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix))&Solver::analyzePattern,
37+
bp::args("self","A"),
38+
"Initializes the iterative solver for the sparsity pattern of the matrix A for further solving Ax=b problems.",
39+
bp::return_self<>())
40+
.def("factorize",
41+
(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix))&Solver::factorize,
42+
bp::args("self","A"),
43+
"Initializes the iterative solver with the numerical values of the matrix A for further solving Ax=b problems.",
44+
bp::return_self<>())
45+
.def("compute",
46+
(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> & matrix))&Solver::compute,
47+
bp::args("self","A"),
48+
"Initializes the iterative solver with the matrix A for further solving Ax=b problems.",
49+
bp::return_self<>())
50+
51+
.def("rows",&Solver::rows,bp::arg("self"),"Returns the number of rows.")
52+
.def("cols",&Solver::cols,bp::arg("self"),"Returns the number of columns.")
53+
.def("tolerance",&Solver::tolerance,bp::arg("self"),
54+
"Returns the tolerance threshold used by the stopping criteria.")
55+
.def("setTolerance",&Solver::setTolerance,bp::args("self","tolerance"),
56+
"Sets the tolerance threshold used by the stopping criteria.\n"
57+
"This value is used as an upper bound to the relative residual error: |Ax-b|/|b|.\n"
58+
"The default value is the machine precision given by NumTraits<Scalar>::epsilon().",
59+
bp::return_self<>())
60+
.def("preconditioner",(Preconditioner & (Solver::*)())&Solver::preconditioner,bp::arg("self"),
61+
"Returns a read-write reference to the preconditioner for custom configuration.",
62+
bp::return_internal_reference<>())
63+
64+
.def("maxIterations",&Solver::maxIterations,bp::arg("self"),
65+
"Returns the max number of iterations.\n"
66+
"It is either the value setted by setMaxIterations or, by default, twice the number of columns of the matrix.")
67+
.def("setMaxIterations",&Solver::setMaxIterations,bp::args("self","max_iterations"),
68+
"Sets the max number of iterations.\n"
69+
"Default is twice the number of columns of the matrix.",
70+
bp::return_self<>())
71+
72+
.def("iterations",&Solver::iterations,bp::arg("self"),
73+
"Returns the number of iterations performed during the last solve.")
74+
.def("error",&Solver::error,bp::arg("self"),
75+
"Returns the tolerance error reached during the last solve.\n"
76+
"It is a close approximation of the true relative residual error |Ax-b|/|b|.")
77+
.def("info",&Solver::error,bp::arg("info"),
78+
"Returns Success if the iterations converged, and NoConvergence otherwise.")
79+
80+
.def("solveWithGuess",&solveWithGuess<MatrixXs,MatrixXs>,bp::args("self","b","x0"),
81+
"Returns the solution x of A x = b using the current decomposition of A and x0 as an initial solution.")
82+
83+
.def("solve",&solve<MatrixXs>,bp::args("self","b"),
84+
"Returns the solution x of A x = b using the current decomposition of A where b is a right hand side matrix or vector.")
85+
;
86+
}
87+
88+
private:
89+
90+
template<typename MatrixOrVector1, typename MatrixOrVector2>
91+
static MatrixOrVector1 solveWithGuess(const Solver & self, const MatrixOrVector1 & b, const MatrixOrVector2 & guess)
92+
{
93+
return self.solveWithGuess(b,guess);
94+
}
95+
96+
template<typename MatrixOrVector>
97+
static MatrixOrVector solve(const Solver & self, const MatrixOrVector & mat_or_vec)
98+
{
99+
MatrixOrVector res = self.solve(mat_or_vec);
100+
return res;
101+
}
102+
};
103+
104+
template<typename _MatrixType>
105+
struct MINRESSolverVisitor
106+
: public boost::python::def_visitor< MINRESSolverVisitor<_MatrixType> >
107+
{
108+
109+
typedef _MatrixType MatrixType;
110+
typedef typename MatrixType::Scalar Scalar;
111+
typedef typename MatrixType::RealScalar RealScalar;
112+
typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1,MatrixType::Options> VectorXs;
113+
typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic,MatrixType::Options> MatrixXs;
114+
typedef Eigen::MINRES<MatrixType> Solver;
115+
116+
template<class PyClass>
117+
void visit(PyClass& cl) const
118+
{
119+
namespace bp = boost::python;
120+
cl
121+
.def(bp::init<>("Default constructor"))
122+
.def(bp::init<MatrixType>(bp::arg("matrix"),
123+
"Initialize the solver with matrix A for further Ax=b solving.\n"
124+
"This constructor is a shortcut for the default constructor followed by a call to compute()."))
125+
126+
.def(IterativeSolverBaseVisitor<Solver>())
127+
;
128+
}
129+
130+
static void expose()
131+
{
132+
static const std::string classname = "MINRES" + scalar_name<Scalar>::shortname();
133+
expose(classname);
134+
}
135+
136+
static void expose(const std::string & name)
137+
{
138+
namespace bp = boost::python;
139+
bp::class_<Solver, boost::noncopyable>(name.c_str(),
140+
"A minimal residual solver for sparse symmetric problems.\n"
141+
"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"
142+
"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",
143+
bp::no_init)
144+
.def(MINRESSolverVisitor());
145+
}
146+
147+
private:
148+
149+
template<typename MatrixOrVector>
150+
static MatrixOrVector solve(const Solver & self, const MatrixOrVector & vec)
151+
{
152+
return self.solve(vec);
153+
}
154+
};
155+
156+
} // namespace eigenpy
157+
158+
#endif // ifndef __eigenpy_decomposition_minres_hpp__

‎src/decompositions/decompositions.cpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
/*
2-
* Copyright 2020 INRIA
2+
* Copyright 2020-2021 INRIA
33
*/
44

55
#include "eigenpy/fwd.hpp"
@@ -10,6 +10,7 @@
1010
#include "eigenpy/decompositions/SelfAdjointEigenSolver.hpp"
1111
#include "eigenpy/decompositions/LLT.hpp"
1212
#include "eigenpy/decompositions/LDLT.hpp"
13+
#include "eigenpy/decompositions/minres.hpp"
1314

1415
namespace eigenpy
1516
{
@@ -22,6 +23,8 @@ namespace eigenpy
2223
SelfAdjointEigenSolverVisitor<MatrixXd>::expose("SelfAdjointEigenSolver");
2324
LLTSolverVisitor<MatrixXd>::expose("LLT");
2425
LDLTSolverVisitor<MatrixXd>::expose("LDLT");
26+
27+
MINRESSolverVisitor<MatrixXd>::expose("MINRES");
2528

2629
{
2730
bp::enum_<DecompositionOptions>("DecompositionOptions")

‎unittest/CMakeLists.txt‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
#
22
# Copyright (c) 2014-2019 CNRS
3-
# Copyright (c) 2018-2020 INRIA
3+
# Copyright (c) 2018-2021 INRIA
44
#
55

66
MACRO(ADD_LIB_UNIT_TEST test)
@@ -62,3 +62,6 @@ SET_TESTS_PROPERTIES("py-LLT" PROPERTIES DEPENDS ${PYWRAP})
6262

6363
ADD_PYTHON_UNIT_TEST("py-LDLT" "unittest/python/test_LDLT.py" "python/eigenpy;unittest")
6464
SET_TESTS_PROPERTIES("py-LDLT" PROPERTIES DEPENDS ${PYWRAP})
65+
66+
ADD_PYTHON_UNIT_TEST("py-MINRES" "unittest/python/test_MINRES.py" "python/eigenpy;unittest")
67+
SET_TESTS_PROPERTIES("py-MINRES" PROPERTIES DEPENDS ${PYWRAP})

‎unittest/python/test_MINRES.py‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
import eigenpy
2+
3+
import numpy as np
4+
import numpy.linalg as la
5+
6+
dim = 100
7+
A = np.eye(dim)
8+
9+
minres = eigenpy.MINRES(A)
10+
11+
X = np.random.rand(dim,20)
12+
B = A.dot(X)
13+
X_est = minres.solve(B)
14+
assert eigenpy.is_approx(A.dot(X_est),B,1e-6)

0 commit comments

Comments
 (0)