/* * 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_ #define UMFPACK_NONE 0 #define SOLVER_NOOPTION 0 #define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1 #define SOLVER_PRINT_DETERMINANT 2 #ifdef HAVE_EIGEN /////// Compiled with EIGEN support #include "Vector/Vector.hpp" #include "Eigen/UmfPackSupport" #include template class umfpack_solver { public: template static Vector solve(const SparseMatrix & A, const Vector & b) { std::cerr << "Error Umfpack only support double precision, and int ad id type" << "/n"; } void best_solve() { std::cerr << "Error Umfpack only support double precision, and int ad id type" << "/n"; } }; template<> class umfpack_solver { public: /*! \brief No nothing * * */ void best_solve() {}; /*! \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 * */ static Vector solve(SparseMatrix & A, const Vector & b, size_t opt = UMFPACK_NONE) { Vcluster & vcl = create_vcluster(); 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"; x.scatter(); 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; } }; #else /////// Compiled without EIGEN support #include "Vector/Vector.hpp" template class umfpack_solver { public: template static Vector solve(const SparseMatrix & A, const Vector & b) { std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n"; } void best_solve() { std::cerr << __FILE__ << ":" << __LINE__ << " Error Umfpack only support double precision" << "/n"; } }; template<> class umfpack_solver { public: template static Vector solve(SparseMatrix & A, const Vector & b, size_t opt = UMFPACK_NONE) { std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n"; Vector x; return x; } void best_solve() { std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n"; } }; #endif #endif /* OPENFPM_NUMERICS_SRC_SOLVERS_UMFPACK_SOLVER_HPP_ */