7 #ifndef LINEARSOLVER_HPP
8 #define LINEARSOLVER_HPP
10 #include <Eigen/Dense>
11 #include <Eigen/Sparse>
14 #include <Eigen/UmfPackSupport>
18 #include <Eigen/PardisoSupport>
26 #include <boost/algorithm/string.hpp>
41 {
"eps_refinement",0},{
"eps_ctrl",1}
55 "Set the solver parameters; Usage [with v1, v2 as doubles]: \n"
56 "--solver-parameters param1=v1 param2=v2\n"
57 "--solver-parameters param1=v1 --solver-parameters param2=v2 \n"
59 "Possible parameters are [all optionals]:\n eps_refinement \n eps_ctrl\n"
61 "Set the solver parameters, not implemented with this solver, see src/Common/linearsolver.hpp for modifications\n"
66 template<
typename MatrixType>
68 template<
typename MatrixType>
71 template<
typename MatrixType>
78 template<
typename MatrixType>
121 template<
typename MatrixType>
122 using ListSolvers = std::tuple<std::unique_ptr<EigenLUType<MatrixType>>,
123 std::unique_ptr<EigenBiCGStabType<MatrixType>>,
124 std::unique_ptr<PardisoLUType<MatrixType>>,
125 std::unique_ptr<UMFPACKType<MatrixType>>,
126 std::unique_ptr<PastixLUType>,
127 std::unique_ptr<PastixLLTType>>;
131 template<
typename MatrixType>
139 : m_name(
map_solver[boost::algorithm::to_lower_copy(namesolver)])
145 std::cout <<
"[linearsolver] Pardiso not available, reverting to Eigen LU" << std::endl;
151 std::cout <<
"[linearsolver] UMFPACK not available, reverting to Eigen LU" << std::endl;
157 std::cout <<
"[linearsolver] PaStiX not available, reverting to Eigen LU" << std::endl;
176 unsigned nb_threads_hint = std::thread::hardware_concurrency();
178 mkl_set_num_threads(nb_threads_hint);
217 return m_info_factorize;
233 std::get<0>(m_list_solvers).get()->analyzePattern(
A);
238 std::get<1>(m_list_solvers).get()->analyzePattern(
A);
244 std::get<2>(m_list_solvers).get()->analyzePattern(
A);
252 std::get<3>(m_list_solvers).get()->analyzePattern(
A);
260 std::get<4>(m_list_solvers).get()->analyzePattern(
A);
266 std::get<5>(m_list_solvers).get()->analyzePattern(
A);
284 std::get<0>(m_list_solvers).get()->factorize(
A);
285 m_info_factorize = std::get<0>(m_list_solvers).get()->info();
290 std::get<1>(m_list_solvers).get()->factorize(
A);
291 m_info_factorize = std::get<1>(m_list_solvers).get()->info();
297 std::get<2>(m_list_solvers).get()->factorize(
A);
298 m_info_factorize = std::get<2>(m_list_solvers).get()->info();
306 std::get<3>(m_list_solvers).get()->factorize(
A);
307 m_info_factorize = std::get<3>(m_list_solvers).get()->info();
315 std::get<4>(m_list_solvers).get()->factorize(
A);
316 m_info_factorize = std::get<4>(m_list_solvers).get()->info();
322 std::get<5>(m_list_solvers).get()->factorize(
A);
323 m_info_factorize = std::get<5>(m_list_solvers).get()->info();
332 _check_info(m_info_factorize,
"factorize");
345 template<
typename VectorType>
354 x = std::get<0>(m_list_solvers).get()->solve(b);
355 m_info_solve = std::get<0>(m_list_solvers).get()->info();
360 x = std::get<1>(m_list_solvers).get()->solve(b);
361 m_info_solve = std::get<1>(m_list_solvers).get()->info();
367 x = std::get<2>(m_list_solvers).get()->solve(b);
368 m_info_solve = std::get<2>(m_list_solvers).get()->info();
376 x = std::get<3>(m_list_solvers).get()->solve(b);
377 m_info_solve = std::get<3>(m_list_solvers).get()->info();
385 x = std::get<4>(m_list_solvers).get()->solve(b);
386 m_info_solve = std::get<4>(m_list_solvers).get()->info();
392 x = std::get<5>(m_list_solvers).get()->solve(b);
393 m_info_solve = std::get<5>(m_list_solvers).get()->info();
402 _check_info(m_info_solve,
"solve");
408 template<
typename VectorType>
419 template <
typename VectorType>
420 inline double residual(MatrixType &
A, VectorType & b, VectorType & x)
421 {
return (
A*x - b).template lpNorm<Eigen::Infinity>() / (b.template lpNorm<Eigen::Infinity>()+1e-8); };
426 void _check_info(Eigen::ComputationInfo
const & info, std::string step)
const
428 if (info != Eigen::Success) {
429 std::cerr <<
"[linearsolver] ERROR in step: " << step << std::endl;
430 std::cerr <<
"[linearsolver] Info: " << info <<
" (see Eigen::ComputationInfo for the meaning)" << std::endl;
437 ListSolvers<MatrixType> m_list_solvers;
438 Eigen::ComputationInfo m_info_factorize;
439 Eigen::ComputationInfo m_info_solve;
Class to invoke the Pastix LLT solver.
Definition: PastixInterface.hpp:543
Class to invoke the Pastix LU solver.
Definition: PastixInterface.hpp:424
This structure is a wrapper to allow a given code to select various linear solvers.
Definition: linearsolver.hpp:133
bool PardisoLUType
Definition: linearsolver.hpp:75
VectorType compute_and_solve(MatrixType &A, VectorType &b)
Perform all operators to solve Ax=b.
Definition: linearsolver.hpp:409
SolverName
Enumeration of all available solvers.
Definition: linearsolver.hpp:36
const char * SolverParamHelper
Create a solver specific description in boost::desc.options.
Definition: linearsolver.hpp:53
LinearSolver(const std::string &namesolver)
Constructor.
Definition: linearsolver.hpp:138
std::map< std::string, SolverName > map_solver
Map to associate to each lowercase name a solver.
Definition: linearsolver.hpp:96
double residual(MatrixType &A, VectorType &b, VectorType &x)
Check relative infinity norm of residual: ||Ax-b||/||b||.
Definition: linearsolver.hpp:420
Eigen::ComputationInfo info_solve() const
Returns the information message after the "solve" step.
Definition: linearsolver.hpp:221
std::string name() const
Returns the name of the solver.
Definition: linearsolver.hpp:209
Eigen::ComputationInfo info_factorize() const
Returns the information message after the "factorize" step.
Definition: linearsolver.hpp:215
void analyzePattern(MatrixType &A)
Analyze the pattern of the matrix.
Definition: linearsolver.hpp:227
std::tuple< std::unique_ptr< EigenLUType< MatrixType > >, std::unique_ptr< EigenBiCGStabType< MatrixType > >, std::unique_ptr< PardisoLUType< MatrixType > >, std::unique_ptr< UMFPACKType< MatrixType > >, std::unique_ptr< PastixLUType >, std::unique_ptr< PastixLLTType > > ListSolvers
Definition: linearsolver.hpp:127
std::map< std::string, int > ParamName
To parse optional Solver parameter from command line.
Definition: linearsolver.hpp:39
std::vector< double > SolverParam
Create space and stores default values of optional solver parameters.
Definition: linearsolver.hpp:46
Eigen::BiCGSTAB< MatrixType, Eigen::IncompleteLUT< double > > EigenBiCGStabType
Definition: linearsolver.hpp:69
bool UMFPACKType
Definition: linearsolver.hpp:82
std::map< SolverName, std::string > map_realname
Map to associate to each solver its proper name.
Definition: linearsolver.hpp:104
Eigen::SparseLU< MatrixType > EigenLUType
Definition: linearsolver.hpp:67
void compute(MatrixType &A)
Analyze and factorize the matrix.
Definition: linearsolver.hpp:337
std::map< SolverName, size_t > map_id
Definition: linearsolver.hpp:111
bool PastixLUType
Definition: linearsolver.hpp:89
bool PastixLLTType
Definition: linearsolver.hpp:90
VectorType solve(VectorType &b)
Solve the system Ax=b using the selected solver (after analysis of pattern and computation of matrix)
Definition: linearsolver.hpp:346
void factorize(MatrixType &A)
Factorize the matrix.
Definition: linearsolver.hpp:278
@ PaStiXLU
Definition: linearsolver.hpp:36
@ UMFPACK
Definition: linearsolver.hpp:36
@ EigenLU
Definition: linearsolver.hpp:36
@ PaStiXLLT
Definition: linearsolver.hpp:36
@ PardisoLU
Definition: linearsolver.hpp:36
@ EigenBiCGStab
Definition: linearsolver.hpp:36
A
Definition: mmread.m:102
Definition: ddr-klplate.hpp:27