/* * Umfpack_solver.hpp * * Created on: Nov 27, 2015 * Author: i-bird */ #ifndef OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ #define OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ #include "Vector/Vector.hpp" #include "Eigen/UmfPackSupport" #include #define UMFPACK_NONE 0 template class umfpack_solver { public: template static Vector solve(const SparseMatrix & A, const Vector & b) { std::cerr << "Error Umfpack only suppor double precision" << "/n"; } }; #define SOLVER_NOOPTION 0 #define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1 #define SOLVER_PRINT_DETERMINANT 2 template<> class umfpack_solver { public: /*! \brief Here we invert the matrix and solve the system * * \warning umfpack is not a parallel solver, this function work only with one processor * * \note if you want to use umfpack in a NON parallel, but on a distributed data, use solve with triplet * * \tparam impl Implementation of the SparseMatrix * */ template static Vector solve(SparseMatrix & A, const Vector & b, size_t opt = UMFPACK_NONE) { Vcluster & vcl = *global_v_cluster; Vector x; // only master processor solve Eigen::UmfPackLU > solver; // Collect the matrix on master auto mat_ei = A.getMat(); Eigen::Matrix x_ei; // Collect the vector on master auto b_ei = b.getVec(); // Copy b into x, this also copy the information on how to scatter back the information on x x = b; if (vcl.getProcessUnitID() == 0) { solver.compute(mat_ei); if(solver.info()!=Eigen::Success) { // decomposition failed std::cout << __FILE__ << ":" << __LINE__ << " solver failed" << "\n"; return x; } x_ei = solver.solve(b_ei); if (opt & SOLVER_PRINT_RESIDUAL_NORM_INFINITY) { Eigen::Matrix res; res = mat_ei * x_ei - b_ei; std::cout << "Infinity norm: " << res.lpNorm() << "\n"; } if (opt & SOLVER_PRINT_DETERMINANT) { std::cout << " Determinant: " << solver.determinant() << "\n"; } x = x_ei; } // Vector is only on master, scatter back the information x.scatter(); return x; } }; #endif /* OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ */