1 #ifndef AMGCL_MAKE_BLOCK_SOLVER_HPP
2 #define AMGCL_MAKE_BLOCK_SOLVER_HPP
3 
4 #include <amgcl/backend/interface.hpp>
5 #include <amgcl/adapter/block_matrix.hpp>
6 #include <amgcl/value_type/static_matrix.hpp>
7 #include <amgcl/make_solver.hpp>
8 #include <amgcl/util.hpp>
9 
10 namespace amgcl {
11 
12 namespace backend {
13 
14 } // namespace backend
15 
16 /* Creates solver that operates in non-scalar domain but may take scalar inputs
17  * for the system matrix and the rhs/solution vectors.
18  */
19 template <class Precond, class IterativeSolver>
20 class make_block_solver {
21     public:
22         typedef typename Precond::backend_type             backend_type;
23         typedef typename backend_type::value_type          value_type;
24         typedef typename backend_type::params              backend_params;
25         typedef typename backend_type::vector              vector;
26         typedef typename math::scalar_of<value_type>::type scalar_type;
27         typedef typename math::rhs_of<value_type>::type    rhs_type;
28 
29         typedef typename make_solver<Precond, IterativeSolver>::params params;
30 
31         template <class Matrix>
make_block_solver(const Matrix & A,const params & prm=params (),const backend_params & bprm=backend_params ())32         make_block_solver(
33                 const Matrix &A,
34                 const params &prm = params(),
35                 const backend_params &bprm = backend_params()
36                 )
37         {
38             S = std::make_shared<Solver>(adapter::block_matrix<value_type>(A), prm, bprm);
39         }
40 
41         template <class Matrix, class Vec1, class Vec2>
operator ()(const Matrix & A,const Vec1 & rhs,Vec2 && x) const42         std::tuple<size_t, scalar_type> operator()(
43                 const Matrix &A, const Vec1 &rhs, Vec2 &&x) const
44         {
45             typedef typename math::scalar_of<typename backend::value_type<typename std::decay<Vec1>::type>::type>::type fs;
46             typedef typename math::scalar_of<typename backend::value_type<typename std::decay<Vec2>::type>::type>::type xs;
47 
48             typedef typename math::replace_scalar<rhs_type, fs>::type f_type;
49             typedef typename math::replace_scalar<rhs_type, xs>::type x_type;
50 
51             auto F = backend::reinterpret<const f_type>(rhs);
52             auto X = backend::reinterpret<x_type>(x);
53 
54             return (*S)(A, F, X);
55         }
56 
57         template <class Vec1, class Vec2>
58         std::tuple<size_t, scalar_type>
operator ()(const Vec1 & rhs,Vec2 && x) const59         operator()(const Vec1 &rhs, Vec2 &&x) const {
60             typedef typename math::scalar_of<typename backend::value_type<typename std::decay<Vec1>::type>::type>::type fs;
61             typedef typename math::scalar_of<typename backend::value_type<typename std::decay<Vec2>::type>::type>::type xs;
62 
63             typedef typename math::replace_scalar<rhs_type, fs>::type f_type;
64             typedef typename math::replace_scalar<rhs_type, xs>::type x_type;
65 
66             auto F = backend::reinterpret<const f_type>(rhs);
67             auto X = backend::reinterpret<x_type>(x);
68 
69             return (*S)(F, X);
70         }
71 
system_matrix_ptr() const72         std::shared_ptr<typename Precond::matrix> system_matrix_ptr() const {
73             return S->system_matrix_ptr();
74         }
75 
system_matrix() const76         typename Precond::matrix const& system_matrix() const {
77             return S->system_matrix();
78         }
79 
operator <<(std::ostream & os,const make_block_solver & p)80         friend std::ostream& operator<<(std::ostream &os, const make_block_solver &p) {
81             return os << *p.S << std::endl;
82         }
83 
bytes() const84         size_t bytes() const {
85             return backend::bytes(*S);
86         }
87     private:
88         typedef make_solver<Precond, IterativeSolver> Solver;
89         std::shared_ptr<Solver> S;
90 };
91 
92 } // namespace amgcl
93 
94 #endif
95