Commit 76d0407d authored by incardon's avatar incardon

Add ORB

parent d4e837f0
......@@ -65,7 +65,8 @@ private:
//! the set of all local sub-domain as vector
data_s<SpaceBox<dim,T>,device_l<SpaceBox<dim,T>>,Memory,openfpm::vector_grow_policy_default, openfpm::vect_isel<SpaceBox<dim,T>>::value > sub_domains;
//! base structure
//! Structure that contain for each sub-domain box the processor id
//! exist for efficient global communication
openfpm::vector<size_t> fine_s;
//! number of total sub-domain
......@@ -128,6 +129,10 @@ private:
met.decompose<nm_part_v::id>();
// fill the structure that store the processor id for each sub-domain
fine_s.resize(N_tot);
// Optimize the decomposition creating bigger spaces
// And reducing Ghost over-stress
......
......@@ -10,8 +10,11 @@ BOOST_AUTO_TEST_SUITE( CartDecomposition_test )
BOOST_AUTO_TEST_CASE( CartDecomposition_test_use)
{
// Virtual cluster
Vcluster vcl(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
// Vcluster
Vcluster & vcl = *global_v_cluster;
// Initialize the global VCluster
init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
CartDecomposition<3,float> dec(vcl);
......
/*
* ORB.hpp
*
* Created on: Mar 13, 2015
* Author: Pietro Incardona
*/
#ifndef ORB_HPP_
#define ORB_HPP_
#include "data_type/scalar.hpp"
/*! \brief this class is a functor for "for_each" algorithm
*
* This class is a functor for "for_each" algorithm. For each
* element of a range_c the operator() is called.
* Is mainly used to unroll the bisect cycle inside one ORB cycle
*
*/
template<unsigned int dim, typename ORB>
struct bisect_unroll
{
ORB & orb;
// start
size_t start = 0;
// level
size_t n = 0;
/*! \brief constructor
*
*/
bisect_unroll(ORB & orb, size_t start, size_t n)
:orb(orb),start(start),n(n)
{};
//! It call the copy function for each property
template<typename T>
void operator()(T& t)
{
start += orb.template bisect<T::value>(n,start);
n = (n+1)%dim;
}
};
/*! \brief This structure use SFINAE to avoid instantiation of invalid code
*
* birect_unroll cannot be performed when i > dim (basically you will never cut
* on a dimension that does not exist)
*
* What is happening here is that at compile-time enable_if< (i < dim) >::type
* is evaluated, if the condition is true enable_if<true>::type is a valid expression
* the the second specialization will be used, if is false enable_if<false>::type
* is not a valid expression, and the first specialization is chosen
*
*/
template<unsigned int dim, unsigned int i, typename ORB, class Enable=void>
struct do_when_dim_gr_i
{
static void bisect_loop(bisect_unroll<dim,ORB> & bu)
{
}
};
template<unsigned int dim, unsigned int i, typename ORB>
struct do_when_dim_gr_i<dim,i,ORB,typename boost::enable_if< boost::mpl::bool_<(i < dim)> >::type>
{
static void bisect_loop(bisect_unroll<dim,ORB> & bu, typename boost::enable_if< boost::mpl::bool_<(i < dim)> >::type* dummy = 0)
{
boost::mpl::for_each< boost::mpl::range_c<int,0,i> >(bu);
}
};
/*! \brief ORB node
* \tparam type of space float, double, complex ...
*
* it store the Center of mass of the local calculated particles
* on one direction (only one direction, mean that one scalar is enough)
*
*/
template<typename T> class ORB_node : public scalar<T>
{
public:
static const unsigned int CM = 0;
};
/*! \brief This class implement orthogonal recursive bisection
*
* \tparam dim Dimensionality of the ORB
* \tparam T type of the space
* \tparam loc_wg type of structure that store the local weight of particles
* \tparam loc_pos type of structure that store the local position of particles
* \tparam Box type of structure that contain the domain extension
* \tparam Tree type of structure that store the tree structure
*
*/
template<unsigned int dim, typename T, typename loc_wg=openfpm::vector<float>, typename loc_pos=openfpm::vector<space<dim,T>> , typename Box=Box<dim,T>, template<typename,typename> class Tree=Graph_CSR_s>
class ORB
{
// Virtual cluster
Vcluster & v_cl;
// particle coordinate accumulator
openfpm::vector<T> cm;
// number of elements summed into cm
openfpm::vector<size_t> cm_cnt;
// Local particle vector
loc_pos & lp;
// Label id of the particle, the label identify where the "local" particles
// are in the local decomposition
openfpm::vector<size_t> lp_lbl;
size_t dim_div;
// Structure that store the orthogonal bisection tree
Tree<ORB_node<T>,no_edge> grp;
/*! \brief Calculate the local center of mass on direction dir
*
* WARNING: with CM we mean the sum of the particle coordinate over one direction
*
* \tparam dir Direction witch to calculate the center of mass
*
* \param start from where the last leafs start
*/
template<unsigned int dir> void local_cm(size_t start)
{
typedef space<dim,T> s;
// reset the counters and accumulators
cm.fill(0);
cm_cnt.fill(0);
// Calculate the local CM
auto it = lp.getIterator();
while (it.isNext())
{
// vector key
auto key = it.get();
size_t lbl = lp_lbl.get(key) - start;
// add the particle coordinate to the CM accumulator
cm.get(lbl) += lp.template get<s::x>(key)[dir];
// increase the counter
cm_cnt.get(lbl)++;
++it;
}
}
/*! \brief It label the particles
*
* It identify where the particles are
*
* \param n level
*
*/
template<unsigned int dir> inline void local_label()
{
typedef space<dim,T> s;
// we just have to create the labels array with zero
if (dir == 0)
{lp_lbl.resize(lp.size());lp_lbl.fill(0); return;}
// we check where (the node) the particles live
auto p_it = lp.getIterator();
while(p_it.isNext())
{
auto key = p_it.get();
// Get the label
size_t lbl = lp_lbl.get(key);
// each node has two child's (FOR SURE)
// No leaf are processed here
size_t n1 = grp.getChild(lbl,0);
size_t n2 = grp.getChild(lbl,1);
// get the splitting center of mass
T cm = grp.template vertex_p<ORB_node<T>::CM>(lbl);
// if the particle n-coordinate is smaller than the CM is inside the child n1
// otherwise is inside the child n2
if (lp.template get<s::x>(key)[dir] < cm)
{lp_lbl.get(key) = n1;}
else
{lp_lbl.get(key) = n2;}
++p_it;
}
}
template<unsigned int dir> size_t bisect(size_t n, size_t start)
{
// first label the local particles
local_label<dir>();
// Index from where start the first leaf
size_t n_node = (n == 0)?1: 2 << (n-1);
// Calculate the local cm
local_cm<dir>(start - n_node);
// reduce the local cm and cm_cnt
v_cl.reduce(cm);
v_cl.reduce(cm_cnt);
v_cl.execute();
// set the CM for the previous leaf (they are not anymore leaf)
for (size_t i = 0 ; i < n_node ; i++)
{
grp.template vertex_p<ORB_node<T>::CM>(start-n_node+i) = cm.get(i) / cm_cnt.get(i);
}
// append on the 2^(n-1) previous end node to the 2^n leaf
for (size_t i = start - n_node, s = 0 ; i < start ; i++, s++)
{
// Add 2 leaf nodes and connect them with the node
grp.addVertex();
grp.addVertex();
grp.template addEdge(i,start+2*s);
grp.template addEdge(i,start+2*s+1);
}
return 2*n_node;
}
// define the friend class bisect_unroll
friend class bisect_unroll<dim,ORB<dim,T,loc_wg,loc_pos,Box,Tree>>;
public:
/*! \brief constructor
*
* \param dom Box domain
* \param n_sub number of sub-domain to create (it is approximated to the biggest 2^n number)
* \param lp Local position of the particles
*
*/
ORB(Box dom, size_t n_sub, loc_pos & lp)
:v_cl(*global_v_cluster),lp(lp)
{
typedef ORB<dim,T,loc_wg,loc_pos,Box,Tree> ORB_class;
dim_div = 0;
n_sub = round_big_2(n_sub);
size_t nsub = log2(n_sub);
// number of center or mass needed
cm.resize(2 << nsub);
cm_cnt.resize(2 << nsub);
// Every time we divide from 0 to dim-1, calculate how many cycle from 0 to dim-1 we have
size_t dim_cycle = nsub / dim;
// Create the root node in the graph
grp.addVertex();
// unroll bisection cycle
bisect_unroll<dim,ORB_class> bu(*this,1,0);
for (size_t i = 0 ; i < dim_cycle ; i++)
{
boost::mpl::for_each< boost::mpl::range_c<int,0,dim> >(bu);
}
// calculate and execute the remaining cycles
switch(nsub - dim_cycle * dim)
{
case 1:
do_when_dim_gr_i<dim,1,ORB_class>::bisect_loop(bu);
break;
case 2:
do_when_dim_gr_i<dim,2,ORB_class>::bisect_loop(bu);
break;
case 3:
do_when_dim_gr_i<dim,3,ORB_class>::bisect_loop(bu);
break;
case 4:
do_when_dim_gr_i<dim,4,ORB_class>::bisect_loop(bu);
break;
case 5:
do_when_dim_gr_i<dim,5,ORB_class>::bisect_loop(bu);
break;
case 6:
do_when_dim_gr_i<dim,6,ORB_class>::bisect_loop(bu);
break;
case 7:
do_when_dim_gr_i<dim,7,ORB_class>::bisect_loop(bu);
break;
default:
// To extend on higher dimension create other cases or create runtime
// version of local_cm and bisect
std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " ORB is not working for dimension bigger than 8";
}
}
};
#endif /* ORB_HPP_ */
/*
* ORB_unit_test.hpp
*
* Created on: Mar 16, 2015
* Author: i-bird
*/
#ifndef ORB_UNIT_TEST_HPP_
#define ORB_UNIT_TEST_HPP_
#include "ORB.hpp"
#include <random>
BOOST_AUTO_TEST_SUITE( ORB_test )
#define N_POINTS 1024
BOOST_AUTO_TEST_CASE( ORB_test_use)
{
// Initialize the global VCluster
init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
// Seed with a real random value, if available
// and create the random generator engine
std::srand(global_v_cluster->getProcessUnitID());
std::default_random_engine eg;
std::uniform_real_distribution<float> ud(0.0f, 1.0f);
typedef space<3,float> p;
// create a random local vector of particles
openfpm::vector<space<3,float>> vp(N_POINTS);
// fill the particles
auto vp_it = vp.getIterator();
while (vp_it.isNext())
{
auto key = vp_it.get();
vp.template get<p::x>(key)[0] = ud(eg);
vp.template get<p::x>(key)[1] = ud(eg);
vp.template get<p::x>(key)[2] = ud(eg);
++vp_it;
}
// Orthogonal Recursive Bisection
Box<3,float> dom({0.0,0.0,0.0},{1.0,1.0,1.0});
ORB<3,float> orb(dom,16,vp);
//
}
BOOST_AUTO_TEST_SUITE_END()
#endif /* ORB_UNIT_TEST_HPP_ */
......@@ -12,142 +12,16 @@ template<unsigned int dim, typename T> class space
{
public:
typedef boost::fusion::vector<> type;
typedef boost::fusion::vector<T[dim]> type;
typedef typename memory_traits_inte<type>::type memory_int;
typedef typename memory_traits_lin<type>::type memory_lin;
type data;
static const unsigned int max_prop = 0;
};
template<typename T> class space<1,T>
{
public:
typedef boost::fusion::vector<T> type;
typedef typename memory_traits_inte<type>::type memory_int;
typedef typename memory_traits_lin<type>::type memory_lin;
typedef T stype;
type data;
static const unsigned int x = 0;
static const unsigned int max_prop = 1;
};
template<typename T> class space<2,T>
{
public:
typedef boost::fusion::vector<T,T> type;
typedef typename memory_traits_inte<type>::type memory_int;
typedef typename memory_traits_lin<type>::type memory_lin;
typedef T stype;
type data;
static const unsigned int x = 0;
static const unsigned int y = 1;
static const unsigned int max_prop = 2;
};
template<typename T> class space<3,T>
{
public:
typedef boost::fusion::vector<T,T,T> type;
typedef typename memory_traits_inte<type>::type memory_int;
typedef typename memory_traits_lin<type>::type memory_lin;
typedef T stype;
type data;
static const unsigned int x = 0;
static const unsigned int y = 1;
static const unsigned int z = 2;
static const unsigned int max_prop = 3;
};
template<typename T> class space<4,T>
{
public:
typedef boost::fusion::vector<T,T,T,T> type;
typedef typename memory_traits_inte<type>::type memory_int;
typedef typename memory_traits_lin<type>::type memory_lin;
typedef T stype;
type data;
static const unsigned int x = 0;
static const unsigned int y = 1;
static const unsigned int z = 2;
static const unsigned int t = 3;
static const unsigned int max_prop = 4;
};
template<typename T> class space<5,T>
{
public:
typedef boost::fusion::vector<T,T,T,T,T> type;
typedef typename memory_traits_inte<type>::type memory_int;
typedef typename memory_traits_lin<type>::type memory_lin;
typedef T stype;
type data;
static const unsigned int max_prop = 5;
};
template<typename T> class space<6,T>
{
public:
typedef boost::fusion::vector<T,T,T,T,T,T> type;
typedef typename memory_traits_inte<type>::type memory_int;
typedef typename memory_traits_lin<type>::type memory_lin;
typedef T stype;
type data;
static const unsigned int max_prop = 6;
};
template<typename T> class space<7,T>
{
public:
typedef boost::fusion::vector<T,T,T,T,T,T,T> type;
typedef typename memory_traits_inte<type>::type memory_int;
typedef typename memory_traits_lin<type>::type memory_lin;
typedef T stype;
type data;
static const unsigned int max_prop = 7;
};
template<typename T> class space<8,T>
{
public:
typedef boost::fusion::vector<T,T,T,T,T,T,T,T> type;
typedef typename memory_traits_inte<type>::type memory_int;
typedef typename memory_traits_lin<type>::type memory_lin;
typedef T stype;
type data;
static const unsigned int max_prop = 8;
};
#endif /* SPACES_HPP_ */
......@@ -28,15 +28,13 @@ BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use)
{
auto key = it.get();
vd.template getPos<space<2,float>::x>(key) = (rand()/(double)(RAND_MAX));
vd.template getPos<space<2,float>::y>(key) = (rand()/(double)(RAND_MAX));
vd.template getPos<space<2,float>::x>(key)[0] = (rand()/(double)(RAND_MAX));
vd.template getPos<space<2,float>::x>(key)[1] = (rand()/(double)(RAND_MAX));
++it;
}
// collect weight
}
BOOST_AUTO_TEST_SUITE_END()
......
......@@ -15,6 +15,7 @@
#include "Space/Shape/Box.hpp"
#include "util.hpp"
#include "Decomposition/ORB_unit_test.hpp"
#include "Graph/CartesianGraphFactory_unit_test.hpp"
#include "metis_util_unit_test.hpp"
#include "dec_optimizer_unit_test.hpp"
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment