7 #ifndef NIHU_FMM_X2X_PRECOMPUTE_HPP_INCLUDED
8 #define NIHU_FMM_X2X_PRECOMPUTE_HPP_INCLUDED
16 #ifdef NIHU_FMM_PARALLEL
20 #include "Eigen/SparseCore"
23 #include <type_traits>
38 template <
class Result,
class ClusterDerived,
class FmmTag>
43 typedef ClusterDerived cluster_t;
45 typedef interaction_lists::list_t list_t;
46 typedef Result result_t;
48 typedef typename std::conditional<
50 >::type container_elem_t;
52 typedef typename std::conditional<
56 typedef unsigned int index_t;
57 typedef Eigen::SparseMatrix<index_t> indices_t;
59 template <
class Operator>
61 : m_tree(op.get_tree())
62 , m_pindices(
new indices_t(m_tree.get_n_clusters(), m_tree.get_n_clusters()))
63 , m_pcontainer(
new container_t(m_tree.get_n_levels()))
65 typedef typename std::decay<Operator>::type operator_t;
66 typedef Eigen::Triplet<index_t> triplet_t;
67 std::vector<triplet_t> triplets, to_compute;
69 std::vector<std::vector<bool> > ready(m_tree.
get_n_levels() + 1);
71 auto tstart = std::chrono::steady_clock::now();
73 for (
size_t to = 0; to < list.size(); ++to)
75 for (
auto from : list[to])
77 cluster_t
const &cto = m_tree[to];
78 cluster_t
const &cfrom = m_tree[from];
79 size_t level = cto.get_level();
80 size_t idx = operator_t::operator_t::unique_idx(cto, cfrom);
81 if (idx >= (*m_pcontainer)[level].size() )
83 (*m_pcontainer)[level].resize(idx + 1);
84 ready[level].resize(idx + 1,
false);
86 if (!ready[level][idx])
88 to_compute.push_back(triplet_t(
int(to),
int(from), idx));
89 ready[level][idx] =
true;
91 triplets.push_back(triplet_t(
int(to),
int(from), idx));
95 #ifdef NIHU_FMM_PARALLEL
96 #pragma omp parallel for schedule(dynamic)
98 for (
int i = 0; i < to_compute.size(); ++i)
100 size_t to = to_compute[i].row();
101 size_t from = to_compute[i].col();
102 size_t idx = to_compute[i].value();
103 size_t level = m_tree[to].get_level();
104 (*m_pcontainer)[level][idx] = op(to, from);
106 #ifdef NIHU_FMM_PARALLEL
110 m_pindices->setFromTriplets(triplets.begin(), triplets.end());
111 auto tend = std::chrono::steady_clock::now();
112 m_assembly_time = std::chrono::duration_cast<std::chrono::microseconds>(tend - tstart).count();
115 result_t
const &operator()(
size_t to,
size_t from)
const
117 size_t level = m_tree[to].get_level();
118 size_t idx = m_pindices->coeff(to, from);
119 return (*m_pcontainer)[level][idx];
122 size_t get_assembly_time()
const
124 return m_assembly_time;
129 std::shared_ptr<indices_t> m_pindices;
130 std::shared_ptr<container_t> m_pcontainer;
131 size_t m_assembly_time;
135 template <
class Operator>
137 create_x2x_precompute(Operator
const &op,
typename interaction_lists::list_t
const &list)
139 typedef typename std::decay<Operator>::type operator_t;
141 typename operator_t::result_t,
142 typename operator_t::cluster_t,
143 typename operator_t::fmm_tag