46 #ifndef _ZOLTAN2_COORDCOMMGRAPH_HPP_ 47 #define _ZOLTAN2_COORDCOMMGRAPH_HPP_ 56 #include "Teuchos_CommHelpers.hpp" 57 #include "Teuchos_Comm.hpp" 58 #include "Teuchos_ArrayViewDecl.hpp" 59 #include "Teuchos_RCPDecl.hpp" 64 #define Z2_ABS(x) ((x) >= 0 ? (x) : -(x)) 69 template <
typename scalar_t,
typename part_t>
83 part_t *minHashIndices;
84 part_t *maxHashIndices;
87 std::vector <part_t> *gridIndices;
89 std::set <part_t> neighbors;
97 maxScalar (
std::numeric_limits<scalar_t>::max()),
98 _EPSILON(
std::numeric_limits<scalar_t>::
epsilon()),
101 gridIndices(0), neighbors()
103 lmins =
new scalar_t [dim];
104 lmaxs =
new scalar_t [dim];
106 minHashIndices =
new part_t [dim];
107 maxHashIndices =
new part_t [dim];
108 gridIndices =
new std::vector <part_t> ();
109 for (
int i = 0; i < dim; ++i){
110 lmins[i] = -this->maxScalar;
111 lmaxs[i] = this->maxScalar;
121 maxScalar (
std::numeric_limits<scalar_t>::max()),
122 _EPSILON(
std::numeric_limits<scalar_t>::
epsilon()),
125 gridIndices(0), neighbors()
127 lmins =
new scalar_t [dim];
128 lmaxs =
new scalar_t [dim];
129 minHashIndices =
new part_t [dim];
130 maxHashIndices =
new part_t [dim];
131 gridIndices =
new std::vector <part_t> ();
132 for (
int i = 0; i < dim; ++i){
146 maxScalar (
std::numeric_limits<scalar_t>::max()),
147 _EPSILON(
std::numeric_limits<scalar_t>::
epsilon()),
150 gridIndices(0), neighbors()
153 lmins =
new scalar_t [dim];
154 lmaxs =
new scalar_t [dim];
155 minHashIndices =
new part_t [dim];
156 maxHashIndices =
new part_t [dim];
157 gridIndices =
new std::vector <part_t> ();
158 scalar_t *othermins = other.
getlmins();
159 scalar_t *othermaxs = other.
getlmaxs();
160 for (
int i = 0; i < dim; ++i){
161 lmins[i] = othermins[i];
162 lmaxs[i] = othermaxs[i];
168 delete []this->lmins;
169 delete [] this->lmaxs;
170 delete []this->minHashIndices;
171 delete [] this->maxHashIndices;
205 for (
int i = 0; i < this->dim; i++)
206 centroid[i] = 0.5 * (this->lmaxs[i] + this->lmins[i]);
213 return this->gridIndices;
219 return &(this->neighbors);
225 if (pointdim != this->dim)
226 throw std::logic_error(
"dim of point must match dim of box");
227 for (
int i = 0; i < pointdim; i++) {
228 if (point[i] < this->lmins[i])
return false;
229 if (point[i] > this->lmaxs[i])
return false;
237 if (cdim != this->dim)
238 throw std::logic_error(
"dim of given box must match dim of box");
242 for (
int i = 0; i < cdim; i++) {
243 if (!((lower[i] >= this->lmins[i] && lower[i] <= this->lmaxs[i])
245 || (upper[i] >= this->lmins[i] && upper[i] <= this->lmaxs[i])
247 || (lower[i] < this->lmins[i] && upper[i] > this->lmaxs[i]))) {
266 for (
int i = 0; i < dim; ++i){
268 if (omins[i] - this->lmaxs[i] > _EPSILON ||
269 this->lmins[i] - omaxs[i] > _EPSILON ) {
272 else if (
Z2_ABS(omins[i] - this->lmaxs[i]) < _EPSILON ||
273 Z2_ABS(this->lmins[i] - omaxs[i]) < _EPSILON ){
283 std::cout <<
"something is wrong: equality:" 284 << equality << std::endl;
293 neighbors.insert(nIndex);
299 if (neighbors.end() != neighbors.find(nIndex)){
310 scalar_t *minMaxBoundaries,
311 scalar_t *sliceSizes,
312 part_t numSlicePerDim
314 for (
int j = 0; j < dim; ++j){
315 scalar_t distance = (lmins[j] - minMaxBoundaries[j]);
317 if (distance > _EPSILON && sliceSizes[j] > _EPSILON){
318 minInd =
static_cast<part_t
>(floor((lmins[j] - minMaxBoundaries[j])/ sliceSizes[j]));
321 if(minInd >= numSlicePerDim){
322 minInd = numSlicePerDim - 1;
325 distance = (lmaxs[j] - minMaxBoundaries[j]);
326 if (distance > _EPSILON && sliceSizes[j] > _EPSILON){
327 maxInd =
static_cast<part_t
>(ceil((lmaxs[j] - minMaxBoundaries[j])/ sliceSizes[j]));
329 if(maxInd >= numSlicePerDim){
330 maxInd = numSlicePerDim - 1;
335 minHashIndices[j] = minInd;
336 maxHashIndices[j] = maxInd;
338 std::vector <part_t> *in =
new std::vector <part_t> ();
340 std::vector <part_t> *out =
new std::vector <part_t> ();
342 for (
int j = 0; j < dim; ++j){
344 part_t minInd = minHashIndices[j];
345 part_t maxInd = maxHashIndices[j];
348 part_t pScale = part_t(pow (
float(numSlicePerDim),
int(dim - j -1)));
350 part_t inSize = in->size();
352 for (part_t k = minInd; k <= maxInd; ++k){
353 for (part_t i = 0; i < inSize; ++i){
354 out->push_back((*in)[i] + k * pScale);
358 std::vector <part_t> *tmp = in;
363 std::vector <part_t> *tmp = in;
375 for(
int i = 0; i < this->dim; ++i){
376 std::cout <<
"\tbox:" << this->pID <<
" dim:" << i <<
" min:" << lmins[i] <<
" max:" << lmaxs[i] << std::endl;
384 lmaxs[dimInd] = newBoundary;
387 lmins[dimInd] = newBoundary;
391 template <
typename tt>
393 std::stringstream ss (std::stringstream::in |std::stringstream::out);
395 std::string tmp =
"";
404 int numCorners = (int(1)<<dim);
405 scalar_t *corner1 =
new scalar_t [dim];
406 scalar_t *corner2 =
new scalar_t [dim];
408 for (
int i = 0; i < dim; ++i){
420 mm << lmins[i] <<
" ";
424 for (
int i = 0; i < dim; ++i){
438 mm << lmaxs[i] <<
" ";
443 for (
int j = 0; j < numCorners; ++j){
444 std::vector <int> neighborCorners;
445 for (
int i = 0; i < dim; ++i){
446 if(
int(j & (
int(1)<<i)) == 0){
447 corner1[i] = lmins[i];
450 corner1[i] = lmaxs[i];
452 if (j % (
int(1)<<(i + 1)) >= (int(1)<<i)){
453 int c1 = j - (int(1)<<i);
456 neighborCorners.push_back(c1);
461 int c1 = j + (int(1)<<i);
462 if (c1 < (
int(1) << dim)) {
463 neighborCorners.push_back(c1);
468 for (
int m = 0; m < int (neighborCorners.size()); ++m){
470 int n = neighborCorners[m];
472 for (
int i = 0; i < dim; ++i){
473 if(
int(n & (
int(1)<<i)) == 0){
474 corner2[i] = lmins[i];
477 corner2[i] = lmaxs[i];
481 std::string arrowline =
"set arrow from ";
482 for (
int i = 0; i < dim - 1; ++i){
483 arrowline += toString<scalar_t>(corner1[i]) +
",";
485 arrowline += toString<scalar_t>(corner1[dim -1]) +
" to ";
487 for (
int i = 0; i < dim - 1; ++i){
488 arrowline += toString<scalar_t>(corner2[i]) +
",";
490 arrowline += toString<scalar_t>(corner2[dim -1]) +
" nohead\n";
506 template <
typename scalar_t,
typename part_t>
510 const RCP < std::vector <Zoltan2::coordinateModelPartBox <scalar_t, part_t> > > pBoxes;
513 scalar_t *minMaxBoundaries;
515 scalar_t *maxMinBoundaries;
517 scalar_t *sliceSizes;
521 part_t numSlicePerDim;
525 std::vector <std::vector <part_t> > grids;
527 ArrayRCP <part_t> comXAdj;
528 ArrayRCP <part_t> comAdj;
535 part_t ntasks_,
int dim_):
538 maxMinBoundaries(0), sliceSizes(0),
541 numSlicePerDim(part_t(pow(double(ntasks_), 1.0 / dim))),
547 minMaxBoundaries =
new scalar_t[dim];
548 maxMinBoundaries =
new scalar_t[dim];
549 sliceSizes =
new scalar_t[dim];
552 if (numSlicePerDim == 0) numSlicePerDim = 1;
554 numGrids = part_t(pow(
float(numSlicePerDim),
int(dim)));
557 std::vector <std::vector <part_t> > grids_ (numGrids);
558 this->grids = grids_;
560 this->getMinMaxBoundaries();
562 this->insertToHash();
564 part_t nCount = this->calculateNeighbors();
567 ArrayRCP <part_t> tmpComXadj(ntasks_+1);
568 ArrayRCP <part_t> tmpComAdj(nCount);
569 comXAdj = tmpComXadj;
572 this->fillAdjArrays();
580 delete []minMaxBoundaries;
581 delete []maxMinBoundaries;
593 for(part_t i = 0; i < this->nTasks; ++i){
594 std::set<part_t> *neigbors = (*pBoxes)[i].getNeighbors();
596 part_t s = neigbors->size();
598 comXAdj[i+1] = comXAdj[i] + s;
599 typedef typename std::set<part_t> mySet;
600 typedef typename mySet::iterator myIT;
602 for (it=neigbors->begin(); it!=neigbors->end(); ++it)
604 comAdj[adjIndex++] = *it;
616 ArrayRCP <part_t> &comXAdj_,
617 ArrayRCP <part_t> &comAdj_){
618 comXAdj_ = this->comXAdj;
619 comAdj_ = this->comAdj;
627 for(part_t i = 0; i < this->nTasks; ++i){
628 std::vector <part_t> *gridIndices =(*pBoxes)[i].getGridIndices();
629 part_t gridCount = gridIndices->size();
631 for (part_t j = 0; j < gridCount; ++j){
632 part_t grid = (*gridIndices)[j];
633 part_t boxCount = grids[grid].size();
634 for (part_t k = 0; k < boxCount; ++k){
635 part_t boxIndex = grids[grid][k];
639 (*pBoxes)[i].addNeighbor(boxIndex);
640 (*pBoxes)[boxIndex].addNeighbor(i);
657 for(part_t i = 0; i < this->nTasks; ++i){
658 (*pBoxes)[i].setMinMaxHashIndices(minMaxBoundaries, sliceSizes, numSlicePerDim);
660 std::vector <part_t> *gridIndices =(*pBoxes)[i].getGridIndices();
662 part_t gridCount = gridIndices->size();
664 for (part_t j = 0; j < gridCount; ++j){
665 part_t grid = (*gridIndices)[j];
668 (grids)[grid].push_back(i);
689 scalar_t *mins = (*pBoxes)[0].getlmins();
690 scalar_t *maxs = (*pBoxes)[0].getlmaxs();
692 for (
int j = 0; j < dim; ++j){
693 minMaxBoundaries[j] = maxs[j];
694 maxMinBoundaries[j] = mins[j];
697 for (part_t i = 1; i < nTasks; ++i){
699 mins = (*pBoxes)[i].getlmins();
700 maxs = (*pBoxes)[i].getlmaxs();
702 for (
int j = 0; j < dim; ++j){
704 if (minMaxBoundaries[j] > maxs[j]){
705 minMaxBoundaries[j] = maxs[j];
707 if (maxMinBoundaries[j] < mins[j]){
708 maxMinBoundaries[j] = mins[j];
714 for (
int j = 0; j < dim; ++j){
715 sliceSizes[j] = (maxMinBoundaries[j] - minMaxBoundaries[j]) / numSlicePerDim;
716 if (sliceSizes[j] < 0) sliceSizes[j] = 0;
GridHash Class, Hashing Class for part boxes.
std::set< part_t > * getNeighbors()
function to get the indices of the neighboring parts.
coordinateModelPartBox(part_t pid, int dim_)
Constructor.
part_t getpId() const
function to get the part id
std::string toString(tt obj)
~GridHash()
GridHash Class, Destructor.
void setpId(part_t pid)
function to set the part id
void getAdjArrays(ArrayRCP< part_t > &comXAdj_, ArrayRCP< part_t > &comAdj_)
GridHash Class, returns the adj arrays.
scalar_t * getlmins() const
function to get minimum values along all dimensions
bool pointInBox(int pointdim, scalar_t *point) const
function to test whether a point is in the box
void fillAdjArrays()
GridHash Class, Function to fill adj arrays.
std::vector< part_t > * getGridIndices()
function to get the indices of the buckets that the part is inserted to
void computeCentroid(scalar_t *centroid) const
compute the centroid of the box
coordinateModelPartBox Class, represents the boundaries of the box which is a result of a geometric p...
void writeGnuPlot(std::ofstream &file, std::ofstream &mm)
function for visualization.
void print()
function to print the boundaries.
void updateMinMax(scalar_t newBoundary, int isMax, int dimInd)
function to update the boundary of the box.
scalar_t * getlmaxs() const
function to get maximum values along all dimensions
GridHash(const RCP< std::vector< Zoltan2::coordinateModelPartBox< scalar_t, part_t > > > &pBoxes_, part_t ntasks_, int dim_)
GridHash Class, Constructor.
part_t calculateNeighbors()
GridHash Class, For each box compares the adjacency against the boxes that are in the same buckets...
int getDim() const
function to set the dimension
void getMinMaxBoundaries()
GridHash Class, calculates the minimum of maximum box boundaries, and maxium of minimum box boundarie...
~coordinateModelPartBox()
Destructor.
bool isNeighborWith(const coordinateModelPartBox< scalar_t, part_t > &other) const
function to check if two boxes are neighbors.
void setMinMaxHashIndices(scalar_t *minMaxBoundaries, scalar_t *sliceSizes, part_t numSlicePerDim)
function to obtain the min and max hash values along all dimensions.
bool isAlreadyNeighbor(part_t nIndex)
function to check if a given part is already in the neighbor list.
void insertToHash()
GridHash Class, For each box calculates the buckets which it should be inserted to.
bool boxesOverlap(int cdim, scalar_t *lower, scalar_t *upper) const
function to test whether this box overlaps a given box
coordinateModelPartBox(const coordinateModelPartBox< scalar_t, part_t > &other)
Copy Constructor deep copy of the maximum and minimum boundaries.
void addNeighbor(part_t nIndex)
function to add a new neighbor to the neighbor list.
coordinateModelPartBox(part_t pid, int dim_, scalar_t *lmi, scalar_t *lma)
Constructor deep copy of the maximum and minimum boundaries.