1 #ifndef DIY_ALGORITHMS_HPP 2 #define DIY_ALGORITHMS_HPP 3 4 #include <vector> 5 6 #include "master.hpp" 7 #include "assigner.hpp" 8 #include "reduce.hpp" 9 #include "reduce-operations.hpp" 10 #include "partners/swap.hpp" 11 12 #include "detail/algorithms/sort.hpp" 13 #include "detail/algorithms/kdtree.hpp" 14 #include "detail/algorithms/kdtree-sampling.hpp" 15 16 #include "log.hpp" 17 18 namespace diy 19 { 20 /** 21 * \ingroup Algorithms 22 * \brief sample sort `values` of each block, store the boundaries between blocks in `samples` 23 */ 24 template<class Block, class T, class Cmp> sort(Master & master,const Assigner & assigner,std::vector<T> Block::* values,std::vector<T> Block::* samples,size_t num_samples,const Cmp & cmp,int k=2,bool samples_only=false)25 void sort(Master& master, //!< master object 26 const Assigner& assigner, //!< assigner object 27 std::vector<T> Block::* values, //!< all values to sort 28 std::vector<T> Block::* samples, //!< (output) boundaries of blocks 29 size_t num_samples, //!< desired number of samples 30 const Cmp& cmp, //!< comparison function 31 int k = 2, //!< k-ary reduction will be used 32 bool samples_only = false) //!< false: results will be all_to_all exchanged; true: only sort but don't exchange results 33 { 34 bool immediate = master.immediate(); 35 master.set_immediate(false); 36 37 // NB: although sorter will go out of scope, its member functions sample() 38 // and exchange() will return functors whose copies get saved inside reduce 39 detail::SampleSort<Block,T,Cmp> sorter(values, samples, cmp, num_samples); 40 41 // swap-reduce to all-gather samples 42 RegularDecomposer<DiscreteBounds> decomposer(1, interval(0,assigner.nblocks()), assigner.nblocks()); 43 RegularSwapPartners partners(decomposer, k); 44 reduce(master, assigner, partners, sorter.sample(), detail::SkipIntermediate(partners.rounds())); 45 46 // all_to_all to exchange the values 47 if (!samples_only) 48 all_to_all(master, assigner, sorter.exchange(), k); 49 50 master.set_immediate(immediate); 51 } 52 53 54 /** 55 * \ingroup Algorithms 56 * \brief sample sort `values` of each block, store the boundaries between blocks in `samples` 57 * shorter version of above sort algorithm with the default less-than comparator used for T 58 * and all_to_all exchange included 59 */ 60 template<class Block, class T> sort(Master & master,const Assigner & assigner,std::vector<T> Block::* values,std::vector<T> Block::* samples,size_t num_samples,int k=2)61 void sort(Master& master, //!< master object 62 const Assigner& assigner, //!< assigner object 63 std::vector<T> Block::* values, //!< all values to sort 64 std::vector<T> Block::* samples, //!< (output) boundaries of blocks 65 size_t num_samples, //!< desired number of samples 66 int k = 2) //!< k-ary reduction will be used 67 { 68 sort(master, assigner, values, samples, num_samples, std::less<T>(), k); 69 } 70 71 /** 72 * \ingroup Algorithms 73 * \brief build a kd-tree and sort a set of points into it (use histograms to determine split values) 74 */ 75 template<class Block, class Point> kdtree(Master & master,const Assigner & assigner,int dim,const ContinuousBounds & domain,std::vector<Point> Block::* points,size_t bins,bool wrap=false)76 void kdtree(Master& master, //!< master object 77 const Assigner& assigner, //!< assigner object 78 int dim, //!< dimensionality 79 const ContinuousBounds& domain, //!< global data extents 80 std::vector<Point> Block::* points, //!< input points to sort into kd-tree 81 size_t bins, //!< number of histogram bins for splitting a dimension 82 bool wrap = false)//!< periodic boundaries in all dimensions 83 { 84 if (assigner.nblocks() & (assigner.nblocks() - 1)) 85 throw std::runtime_error(fmt::format("KD-tree requires a number of blocks that's a power of 2, got {}", assigner.nblocks())); 86 87 typedef diy::RegularContinuousLink RCLink; 88 89 for (size_t i = 0; i < master.size(); ++i) 90 { 91 RCLink* link = static_cast<RCLink*>(master.link(i)); 92 *link = RCLink(dim, domain, domain); 93 94 if (wrap) // set up the links to self 95 { 96 diy::BlockID self = { master.gid(i), master.communicator().rank() }; 97 for (int j = 0; j < dim; ++j) 98 { 99 diy::Direction dir(dim,0), wrap_dir(dim,0); 100 101 // left 102 dir[j] = -1; wrap_dir[j] = -1; 103 link->add_neighbor(self); 104 link->add_bounds(domain); 105 link->add_direction(dir); 106 link->add_wrap(wrap_dir); 107 108 // right 109 dir[j] = 1; wrap_dir[j] = 1; 110 link->add_neighbor(self); 111 link->add_bounds(domain); 112 link->add_direction(dir); 113 link->add_wrap(wrap_dir); 114 } 115 } 116 } 117 118 detail::KDTreePartition<Block,Point> kdtree_partition(dim, points, bins); 119 120 detail::KDTreePartners partners(dim, assigner.nblocks(), wrap, domain); 121 reduce(master, assigner, partners, kdtree_partition); 122 123 // update master.expected to match the links 124 int expected = 0; 125 for (size_t i = 0; i < master.size(); ++i) 126 expected += master.link(i)->size_unique(); 127 master.set_expected(expected); 128 } 129 130 /** 131 * \ingroup Algorithms 132 * \brief build a kd-tree and sort a set of points into it (use sampling to determine split values) 133 */ 134 template<class Block, class Point> kdtree_sampling(Master & master,const Assigner & assigner,int dim,const ContinuousBounds & domain,std::vector<Point> Block::* points,size_t samples,bool wrap=false)135 void kdtree_sampling 136 (Master& master, //!< master object 137 const Assigner& assigner, //!< assigner object 138 int dim, //!< dimensionality 139 const ContinuousBounds& domain, //!< global data extents 140 std::vector<Point> Block::* points, //!< input points to sort into kd-tree 141 size_t samples, //!< number of samples to take in each block 142 bool wrap = false)//!< periodic boundaries in all dimensions 143 { 144 if (assigner.nblocks() & (assigner.nblocks() - 1)) 145 throw std::runtime_error(fmt::format("KD-tree requires a number of blocks that's a power of 2, got {}", assigner.nblocks())); 146 147 typedef diy::RegularContinuousLink RCLink; 148 149 for (size_t i = 0; i < master.size(); ++i) 150 { 151 RCLink* link = static_cast<RCLink*>(master.link(i)); 152 *link = RCLink(dim, domain, domain); 153 154 if (wrap) // set up the links to self 155 { 156 diy::BlockID self = { master.gid(i), master.communicator().rank() }; 157 for (int j = 0; j < dim; ++j) 158 { 159 diy::Direction dir(dim,0), wrap_dir(dim,0); 160 161 // left 162 dir[j] = -1; wrap_dir[j] = -1; 163 link->add_neighbor(self); 164 link->add_bounds(domain); 165 link->add_direction(dir); 166 link->add_wrap(wrap_dir); 167 168 // right 169 dir[j] = 1; wrap_dir[j] = 1; 170 link->add_neighbor(self); 171 link->add_bounds(domain); 172 link->add_direction(dir); 173 link->add_wrap(wrap_dir); 174 } 175 } 176 } 177 178 detail::KDTreeSamplingPartition<Block,Point> kdtree_partition(dim, points, samples); 179 180 detail::KDTreePartners partners(dim, assigner.nblocks(), wrap, domain); 181 reduce(master, assigner, partners, kdtree_partition); 182 183 // update master.expected to match the links 184 int expected = 0; 185 for (size_t i = 0; i < master.size(); ++i) 186 expected += master.link(i)->size_unique(); 187 master.set_expected(expected); 188 } 189 } 190 191 #endif 192