Zoltan2
Zoltan2_CoordinatePartitioningGraph.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Zoltan2: A package of combinatorial algorithms for scientific computing
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Karen Devine (kddevin@sandia.gov)
39 // Erik Boman (egboman@sandia.gov)
40 // Siva Rajamanickam (srajama@sandia.gov)
41 //
42 // ***********************************************************************
43 //
44 // @HEADER
45 
46 #ifndef _ZOLTAN2_COORDCOMMGRAPH_HPP_
47 #define _ZOLTAN2_COORDCOMMGRAPH_HPP_
48 
49 
50 #include <cmath>
51 #include <limits>
52 #include <iostream>
53 #include <vector>
54 #include <set>
55 #include <fstream>
56 #include "Teuchos_CommHelpers.hpp"
57 #include "Teuchos_Comm.hpp"
58 #include "Teuchos_ArrayViewDecl.hpp"
59 #include "Teuchos_RCPDecl.hpp"
60 
61 namespace Zoltan2{
62 
63 
64 #define Z2_ABS(x) ((x) >= 0 ? (x) : -(x))
65 
69 template <typename scalar_t,typename part_t>
71 
72  part_t pID; //part Id
73  int dim; //dimension of the box
74  scalar_t *lmins; //minimum boundaries of the box along all dimensions.
75  scalar_t *lmaxs; //maximum boundaries of the box along all dimensions.
76  scalar_t maxScalar;
77  scalar_t _EPSILON;
78 
79  //to calculate the neighbors of the box and avoid the p^2 comparisons,
80  //we use hashing. A box can be put into multiple hash buckets.
81  //the following 2 variable holds the minimum and maximum of the
82  //hash values along all dimensions.
83  part_t *minHashIndices;
84  part_t *maxHashIndices;
85 
86  //result hash bucket indices.
87  std::vector <part_t> *gridIndices;
88  //neighbors of the box.
89  std::set <part_t> neighbors;
90 public:
93  coordinateModelPartBox(part_t pid, int dim_):
94  pID(pid),
95  dim(dim_),
96  lmins(0), lmaxs(0),
97  maxScalar (std::numeric_limits<scalar_t>::max()),
98  _EPSILON(std::numeric_limits<scalar_t>::epsilon()),
99  minHashIndices(0),
100  maxHashIndices(0),
101  gridIndices(0), neighbors()
102  {
103  lmins = new scalar_t [dim];
104  lmaxs = new scalar_t [dim];
105 
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;
112  }
113  }
117  coordinateModelPartBox(part_t pid, int dim_, scalar_t *lmi, scalar_t *lma):
118  pID(pid),
119  dim(dim_),
120  lmins(0), lmaxs(0),
121  maxScalar (std::numeric_limits<scalar_t>::max()),
122  _EPSILON(std::numeric_limits<scalar_t>::epsilon()),
123  minHashIndices(0),
124  maxHashIndices(0),
125  gridIndices(0), neighbors()
126  {
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){
133  lmins[i] = lmi[i];
134  lmaxs[i] = lma[i];
135  }
136  }
137 
138 
143  pID(other.getpId()),
144  dim(other.getDim()),
145  lmins(0), lmaxs(0),
146  maxScalar (std::numeric_limits<scalar_t>::max()),
147  _EPSILON(std::numeric_limits<scalar_t>::epsilon()),
148  minHashIndices(0),
149  maxHashIndices(0),
150  gridIndices(0), neighbors()
151  {
152 
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];
163  }
164  }
168  delete []this->lmins;
169  delete [] this->lmaxs;
170  delete []this->minHashIndices;
171  delete [] this->maxHashIndices;
172  delete gridIndices;
173  }
174 
177  void setpId(part_t pid){
178  this->pID = pid;
179  }
182  part_t getpId() const{
183  return this->pID;
184  }
185 
186 
189  int getDim()const{
190  return this->dim;
191  }
194  scalar_t * getlmins()const{
195  return this->lmins;
196  }
199  scalar_t * getlmaxs()const{
200  return this->lmaxs;
201  }
204  void computeCentroid(scalar_t *centroid)const {
205  for (int i = 0; i < this->dim; i++)
206  centroid[i] = 0.5 * (this->lmaxs[i] + this->lmins[i]);
207  }
208 
212  std::vector <part_t> * getGridIndices () {
213  return this->gridIndices;
214  }
215 
218  std::set<part_t> *getNeighbors() {
219  return &(this->neighbors);
220  }
221 
224  bool pointInBox(int pointdim, scalar_t *point) const {
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;
230  }
231  return true;
232  }
233 
236  bool boxesOverlap(int cdim, scalar_t *lower, scalar_t *upper) const {
237  if (cdim != this->dim)
238  throw std::logic_error("dim of given box must match dim of box");
239 
240  // Check for at least partial overlap
241  bool found = true;
242  for (int i = 0; i < cdim; i++) {
243  if (!((lower[i] >= this->lmins[i] && lower[i] <= this->lmaxs[i])
244  // lower i-coordinate in the box
245  || (upper[i] >= this->lmins[i] && upper[i] <= this->lmaxs[i])
246  // upper i-coordinate in the box
247  || (lower[i] < this->lmins[i] && upper[i] > this->lmaxs[i]))) {
248  // i-coordinates straddle the box
249  found = false;
250  break;
251  }
252  }
253  return found;
254  }
255 
259  const coordinateModelPartBox <scalar_t, part_t> &other) const{
260 
261 
262  scalar_t *omins = other.getlmins();
263  scalar_t *omaxs = other.getlmaxs();
264 
265  int equality = 0;
266  for (int i = 0; i < dim; ++i){
267 
268  if (omins[i] - this->lmaxs[i] > _EPSILON ||
269  this->lmins[i] - omaxs[i] > _EPSILON ) {
270  return false;
271  }
272  else if (Z2_ABS(omins[i] - this->lmaxs[i]) < _EPSILON ||
273  Z2_ABS(this->lmins[i] - omaxs[i]) < _EPSILON ){
274  if (++equality > 1){
275  return false;
276  }
277  }
278  }
279  if (equality == 1) {
280  return true;
281  }
282  else {
283  std::cout << "something is wrong: equality:"
284  << equality << std::endl;
285  return false;
286  }
287  }
288 
289 
292  void addNeighbor(part_t nIndex){
293  neighbors.insert(nIndex);
294  }
297  bool isAlreadyNeighbor(part_t nIndex){
298 
299  if (neighbors.end() != neighbors.find(nIndex)){
300  return true;
301  }
302  return false;
303 
304  }
305 
306 
310  scalar_t *minMaxBoundaries,
311  scalar_t *sliceSizes,
312  part_t numSlicePerDim
313  ){
314  for (int j = 0; j < dim; ++j){
315  scalar_t distance = (lmins[j] - minMaxBoundaries[j]);
316  part_t minInd = 0;
317  if (distance > _EPSILON && sliceSizes[j] > _EPSILON){
318  minInd = static_cast<part_t>(floor((lmins[j] - minMaxBoundaries[j])/ sliceSizes[j]));
319  }
320 
321  if(minInd >= numSlicePerDim){
322  minInd = numSlicePerDim - 1;
323  }
324  part_t maxInd = 0;
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]));
328  }
329  if(maxInd >= numSlicePerDim){
330  maxInd = numSlicePerDim - 1;
331  }
332 
333  //cout << "j:" << j << " lmins:" << lmins[j] << " lmaxs:" << lmaxs[j] << endl;
334  //cout << "j:" << j << " min:" << minInd << " max:" << maxInd << endl;
335  minHashIndices[j] = minInd;
336  maxHashIndices[j] = maxInd;
337  }
338  std::vector <part_t> *in = new std::vector <part_t> ();
339  in->push_back(0);
340  std::vector <part_t> *out = new std::vector <part_t> ();
341 
342  for (int j = 0; j < dim; ++j){
343 
344  part_t minInd = minHashIndices[j];
345  part_t maxInd = maxHashIndices[j];
346 
347 
348  part_t pScale = part_t(pow (float(numSlicePerDim), int(dim - j -1)));
349 
350  part_t inSize = in->size();
351 
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);
355  }
356  }
357  in->clear();
358  std::vector <part_t> *tmp = in;
359  in= out;
360  out= tmp;
361  }
362 
363  std::vector <part_t> *tmp = in;
364  in = gridIndices;
365  gridIndices = tmp;
366 
367 
368  delete in;
369  delete out;
370  }
371 
374  void print(){
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;
377  }
378  }
379 
382  void updateMinMax (scalar_t newBoundary, int isMax, int dimInd){
383  if (isMax){
384  lmaxs[dimInd] = newBoundary;
385  }
386  else {
387  lmins[dimInd] = newBoundary;
388  }
389  }
390 
391  template <typename tt>
392  std::string toString(tt obj){
393  std::stringstream ss (std::stringstream::in |std::stringstream::out);
394  ss << obj;
395  std::string tmp = "";
396  ss >> tmp;
397  return tmp;
398  }
399 
400 
403  void writeGnuPlot(std::ofstream &file,std::ofstream &mm){
404  int numCorners = (int(1)<<dim);
405  scalar_t *corner1 = new scalar_t [dim];
406  scalar_t *corner2 = new scalar_t [dim];
407 
408  for (int i = 0; i < dim; ++i){
409  /*
410  if (-maxScalar == lmins[i]){
411  if (lmaxs[i] > 0){
412  lmins[i] = lmaxs[i] / 2;
413  }
414  else{
415  lmins[i] = lmaxs[i] * 2;
416  }
417  }
418  */
419  //std::cout << lmins[i] << " ";
420  mm << lmins[i] << " ";
421  }
422  //std::cout << std::endl;
423  mm << std::endl;
424  for (int i = 0; i < dim; ++i){
425 
426  /*
427  if (maxScalar == lmaxs[i]){
428  if (lmins[i] < 0){
429  lmaxs[i] = lmins[i] / 2;
430  }
431  else{
432  lmaxs[i] = lmins[i] * 2;
433  }
434  }
435  */
436 
437  //std::cout << lmaxs[i] << " ";
438  mm << lmaxs[i] << " ";
439  }
440  //std::cout << std::endl;
441  mm << std::endl;
442 
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];
448  }
449  else {
450  corner1[i] = lmaxs[i];
451  }
452  if (j % (int(1)<<(i + 1)) >= (int(1)<<i)){
453  int c1 = j - (int(1)<<i);
454 
455  if (c1 > 0) {
456  neighborCorners.push_back(c1);
457  }
458  }
459  else {
460 
461  int c1 = j + (int(1)<<i);
462  if (c1 < (int(1) << dim)) {
463  neighborCorners.push_back(c1);
464  }
465  }
466  }
467  //std::cout << "me:" << j << " nc:" << int (neighborCorners.size()) << std::endl;
468  for (int m = 0; m < int (neighborCorners.size()); ++m){
469 
470  int n = neighborCorners[m];
471  //std::cout << "me:" << j << " n:" << n << std::endl;
472  for (int i = 0; i < dim; ++i){
473  if(int(n & (int(1)<<i)) == 0){
474  corner2[i] = lmins[i];
475  }
476  else {
477  corner2[i] = lmaxs[i];
478  }
479  }
480 
481  std::string arrowline = "set arrow from ";
482  for (int i = 0; i < dim - 1; ++i){
483  arrowline += toString<scalar_t>(corner1[i]) + ",";
484  }
485  arrowline += toString<scalar_t>(corner1[dim -1]) + " to ";
486 
487  for (int i = 0; i < dim - 1; ++i){
488  arrowline += toString<scalar_t>(corner2[i]) + ",";
489  }
490  arrowline += toString<scalar_t>(corner2[dim -1]) + " nohead\n";
491 
492  file << arrowline;
493  }
494  }
495  delete []corner1;
496  delete []corner2;
497  }
498 
499 
500 };
501 
502 
506 template <typename scalar_t, typename part_t>
507 class GridHash{
508 private:
509 
510  const RCP < std::vector <Zoltan2::coordinateModelPartBox <scalar_t, part_t> > > pBoxes;
511 
512  //minimum of the maximum box boundaries
513  scalar_t *minMaxBoundaries;
514  //maximum of the minimum box boundaries
515  scalar_t *maxMinBoundaries;
516  //the size of each slice along dimensions
517  scalar_t *sliceSizes;
518  part_t nTasks;
519  int dim;
520  //the number of slices per dimension
521  part_t numSlicePerDim;
522  //the number of grids - buckets
523  part_t numGrids;
524  //hash vector
525  std::vector <std::vector <part_t> > grids;
526  //result communication graph.
527  ArrayRCP <part_t> comXAdj;
528  ArrayRCP <part_t> comAdj;
529 public:
530 
534  GridHash(const RCP < std::vector <Zoltan2::coordinateModelPartBox <scalar_t, part_t> > > &pBoxes_,
535  part_t ntasks_, int dim_):
536  pBoxes(pBoxes_),
537  minMaxBoundaries(0),
538  maxMinBoundaries(0), sliceSizes(0),
539  nTasks(ntasks_),
540  dim(dim_),
541  numSlicePerDim(part_t(pow(double(ntasks_), 1.0 / dim))),
542  numGrids(0),
543  grids(),
544  comXAdj(), comAdj()
545  {
546 
547  minMaxBoundaries = new scalar_t[dim];
548  maxMinBoundaries = new scalar_t[dim];
549  sliceSizes = new scalar_t[dim];
550  //calculate the number of slices in each dimension.
551  numSlicePerDim /= 2;
552  if (numSlicePerDim == 0) numSlicePerDim = 1;
553 
554  numGrids = part_t(pow(float(numSlicePerDim), int(dim)));
555 
556  //allocate memory for buckets.
557  std::vector <std::vector <part_t> > grids_ (numGrids);
558  this->grids = grids_;
559  //get the boundaries of buckets.
560  this->getMinMaxBoundaries();
561  //insert boxes to buckets
562  this->insertToHash();
563  //calculate the neighbors for each bucket.
564  part_t nCount = this->calculateNeighbors();
565 
566  //allocate memory for communication graph
567  ArrayRCP <part_t> tmpComXadj(ntasks_+1);
568  ArrayRCP <part_t> tmpComAdj(nCount);
569  comXAdj = tmpComXadj;
570  comAdj = tmpComAdj;
571  //fill communication graph
572  this->fillAdjArrays();
573  }
574 
575 
580  delete []minMaxBoundaries;
581  delete []maxMinBoundaries;
582  delete []sliceSizes;
583  }
584 
589 
590  part_t adjIndex = 0;
591 
592  comXAdj[0] = 0;
593  for(part_t i = 0; i < this->nTasks; ++i){
594  std::set<part_t> *neigbors = (*pBoxes)[i].getNeighbors();
595 
596  part_t s = neigbors->size();
597 
598  comXAdj[i+1] = comXAdj[i] + s;
599  typedef typename std::set<part_t> mySet;
600  typedef typename mySet::iterator myIT;
601  myIT it;
602  for (it=neigbors->begin(); it!=neigbors->end(); ++it)
603 
604  comAdj[adjIndex++] = *it;
605  //TODO not needed anymore.
606  neigbors->clear();
607  }
608  }
609 
610 
611 
616  ArrayRCP <part_t> &comXAdj_,
617  ArrayRCP <part_t> &comAdj_){
618  comXAdj_ = this->comXAdj;
619  comAdj_ = this->comAdj;
620  }
621 
626  part_t nCount = 0;
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();
630 
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];
636  if (boxIndex > i){
637  if((!(*pBoxes)[i].isAlreadyNeighbor(boxIndex))&& (*pBoxes)[i].isNeighborWith((*pBoxes)[boxIndex])){
638  //cout << "i:" << i << " n:" << boxIndex << " are neighbors."<< endl;
639  (*pBoxes)[i].addNeighbor(boxIndex);
640  (*pBoxes)[boxIndex].addNeighbor(i);
641  nCount += 2;
642  }
643  }
644  }
645  }
646  }
647 
648  return nCount;
649  }
650 
654  void insertToHash(){
655 
656  //cout << "ntasks:" << this->nTasks << endl;
657  for(part_t i = 0; i < this->nTasks; ++i){
658  (*pBoxes)[i].setMinMaxHashIndices(minMaxBoundaries, sliceSizes, numSlicePerDim);
659 
660  std::vector <part_t> *gridIndices =(*pBoxes)[i].getGridIndices();
661 
662  part_t gridCount = gridIndices->size();
663  //cout << "i:" << i << " gridsize:" << gridCount << endl;
664  for (part_t j = 0; j < gridCount; ++j){
665  part_t grid = (*gridIndices)[j];
666 
667  //cout << "i:" << i << " is being inserted to:" << grid << endl;
668  (grids)[grid].push_back(i);
669  }
670  }
671 
672 
673 /*
674  for(part_t i = 0; i < grids.size(); ++i){
675  cout << "grid:" << i << " gridsuze:" << (grids)[i].size() << " elements:";
676  for(part_t j = 0; j < (grids)[i].size(); ++j){
677  cout <<(grids)[i][j] << " ";
678  }
679  cout << endl;
680 
681  }
682 */
683  }
684 
689  scalar_t *mins = (*pBoxes)[0].getlmins();
690  scalar_t *maxs = (*pBoxes)[0].getlmaxs();
691 
692  for (int j = 0; j < dim; ++j){
693  minMaxBoundaries[j] = maxs[j];
694  maxMinBoundaries[j] = mins[j];
695  }
696 
697  for (part_t i = 1; i < nTasks; ++i){
698 
699  mins = (*pBoxes)[i].getlmins();
700  maxs = (*pBoxes)[i].getlmaxs();
701 
702  for (int j = 0; j < dim; ++j){
703 
704  if (minMaxBoundaries[j] > maxs[j]){
705  minMaxBoundaries[j] = maxs[j];
706  }
707  if (maxMinBoundaries[j] < mins[j]){
708  maxMinBoundaries[j] = mins[j];
709  }
710  }
711  }
712 
713 
714  for (int j = 0; j < dim; ++j){
715  sliceSizes[j] = (maxMinBoundaries[j] - minMaxBoundaries[j]) / numSlicePerDim;
716  if (sliceSizes[j] < 0) sliceSizes[j] = 0;
717  /*
718  cout << "dim:" << j <<
719  " minMax:" << minMaxBoundaries[j] <<
720  " maxMin:" << maxMinBoundaries[j] <<
721  " sliceSizes:" << sliceSizes[j] << endl;
722  */
723  }
724  }
725 };
726 /*
727 template <typename scalar_t,typename part_t>
728 class coordinatePartBox{
729 public:
730  part_t pID;
731  int dim;
732  int numCorners;
733  scalar_t **corners;
734  scalar_t *lmins, *gmins;
735  scalar_t *lmaxs, *gmaxs;
736  scalar_t maxScalar;
737  std::vector <part_t> hash_indices;
738  coordinatePartBox(part_t pid, int dim_, scalar_t *lMins, scalar_t *gMins,
739  scalar_t *lMaxs, scalar_t *gMaxs):
740  pID(pid),
741  dim(dim_),
742  numCorners(int(pow(2, dim_))),
743  corners(0),
744  lmins(lMins), gmins(gMins), lmaxs(lMaxs), gmaxs(gMaxs),
745  maxScalar (std::numeric_limits<scalar_t>::max()){
746  this->corners = new scalar_t *[dim];
747  for (int i = 0; i < dim; ++i){
748  this->corners[i] = new scalar_t[this->numCorners];
749  lmins[i] = this->maxScalar;
750  lmaxs[i] = -this->maxScalar;
751  }
752 
753 
754  for (int j = 0; j < this->numCorners; ++j){
755  for (int i = 0; i < dim; ++i){
756  std::cout << "j:" << j << " i:" << i << " 2^i:" << pow(2,i) << " and:" << int(j & int(pow(2,i))) << std::endl;
757  if(int(j & int(pow(2,i))) == 0){
758  corners[i][j] = gmins[i];
759  }
760  else {
761  corners[i][j] = gmaxs[i];
762  }
763 
764  }
765  }
766  }
767 
768 };
769 
770 template <typename Adapter, typename part_t>
771 class CoordinateCommGraph{
772 private:
773 
774  typedef typename Adapter::lno_t lno_t;
775  typedef typename Adapter::gno_t gno_t;
776  typedef typename Adapter::scalar_t scalar_t;
777 
778  const Environment *env;
779  const Teuchos::Comm<int> *comm;
780  const Zoltan2::CoordinateModel<typename Adapter::base_adapter_t> *coords;
781  const Zoltan2::PartitioningSolution<Adapter> *soln;
782  std::vector<coordinatePartBox, part_t> cpb;
783  int coordDim;
784  part_t numParts;
785 
786 
787 public:
788 
789  CoordinateCommGraph(
790  const Environment *env_,
791  const Teuchos::Comm<int> *comm_,
792  const Zoltan2::CoordinateModel<typename Adapter::base_adapter_t> *coords_,
793  const Zoltan2::PartitioningSolution<Adapter> *soln_
794  ):
795  env(env_),
796  comm(comm_),
797  coords(coords_),
798  soln(soln_),
799  coordDim (coords_->getCoordinateDim()),
800  numParts (this->soln->getActualGlobalNumberOfParts())
801  {
802  this->create_part_boxes();
803  this->hash_part_boxes();
804  this->find_neighbors();
805  }
806 
807  void create_part_boxes(){
808 
809 
810  size_t allocSize = numParts * coordDim;
811  scalar_t *lmins = new scalar_t [allocSize];
812  scalar_t *gmins = new scalar_t [allocSize];
813  scalar_t *lmaxs = new scalar_t [allocSize];
814  scalar_t *gmaxs = new scalar_t [allocSize];
815 
816  for(part_t i = 0; i < numParts; ++i){
817  coordinatePartBox tmp(
818  i,
819  this->coordDim,
820  lmins + i * coordDim,
821  gmins + i * coordDim,
822  lmaxs + i * coordDim,
823  gmaxs + i * coordDim
824  );
825  cpb.push_back(tmp);
826  }
827 
828  typedef StridedData<lno_t, scalar_t> input_t;
829  Teuchos::ArrayView<const gno_t> gnos;
830  Teuchos::ArrayView<input_t> xyz;
831  Teuchos::ArrayView<input_t> wgts;
832  coords->getCoordinates(gnos, xyz, wgts);
833 
834  //local and global num coordinates.
835  lno_t numLocalCoords = coords->getLocalNumCoordinates();
836 
837  scalar_t **pqJagged_coordinates = new scalar_t *[coordDim];
838 
839  for (int dim=0; dim < coordDim; dim++){
840  Teuchos::ArrayRCP<const scalar_t> ar;
841  xyz[dim].getInputArray(ar);
842  //pqJagged coordinate values assignment
843  pqJagged_coordinates[dim] = (scalar_t *)ar.getRawPtr();
844  }
845 
846  part_t *sol_part = soln->getPartList();
847  for(lno_t i = 0; i < numLocalCoords; ++i){
848  part_t p = sol_part[i];
849  cpb[p].updateMinMax(pqJagged_coordinates, i);
850  }
851  delete []pqJagged_coordinates;
852 
853 
854  reduceAll<int, gno_t>(*comm, Teuchos::REDUCE_MIN,
855  dim * numParts, lmins, gmins
856  );
857  reduceAll<int, gno_t>(*comm, Teuchos::REDUCE_MAX,
858  dim * numParts, lmaxs, gmaxs
859  );
860  }
861 
862  void hash_part_boxes (){
863  part_t pSingleDim = pow(double(numParts), double(1.0 / coordDim));
864  if (pSingleDim == 0) pSingleDim = 1;
865  std::vector < std::vector <part_t> > hash
866  (
867  part_t ( pow ( part_t (pSingleDim),
868  part_t(coordDim)
869  )
870  )
871  );
872 
873  //calculate the corners of the dataset.
874  scalar_t *allMins = new scalar_t [coordDim];
875  scalar_t *allMaxs = new scalar_t [coordDim];
876  part_t *hash_scales= new scalar_t [coordDim];
877 
878  for (int j = 0; j < coordDim; ++j){
879  allMins[j] = cpb[0].gmins[j];
880  allMaxs[j] = cpb[0].gmaxs[j];
881  hash_scales[j] = part_t ( pow ( part_t (pSingleDim), part_t(coordDim - j - 1)));
882  }
883 
884  for (part_t i = 1; i < numParts; ++i){
885  for (int j = 0; j < coordDim; ++j){
886  scalar_t minC = cpb[i].gmins[i];
887  scalar_t maxC = cpb[i].gmaxs[i];
888  if (minC < allMins[j]) allMins[j] = minC;
889  if (maxC > allMaxs[j]) allMaxs[j] = maxC;
890  }
891  }
892 
893  //get size of each hash for each dimension
894  scalar_t *hash_slices_size = new scalar_t [coordDim];
895  for (int j = 0; j < coordDim; ++j){
896  hash_slices_size[j] = (allMaxs[j] - allMins[j]) / pSingleDim;
897 
898  }
899 
900  delete []allMaxs;
901  delete []allMins;
902 
903 
904 
905  std::vector <part_t> *hashIndices = new std::vector <part_t>();
906  std::vector <part_t> *resultHashIndices = new std::vector <part_t>();
907  std::vector <part_t> *tmp_swap;
908  for (part_t i = 0; i < numParts; ++i){
909  hashIndices->clear();
910  resultHashIndices->clear();
911 
912  hashIndices->push_back(0);
913 
914  for (int j = 0; j < coordDim; ++j){
915 
916  scalar_t minC = cpb[i].gmins[i];
917  scalar_t maxC = cpb[i].gmaxs[i];
918  part_t minHashIndex = part_t ((minC - allMins[j]) / hash_slices_size[j]);
919  part_t maxHashIndex = part_t ((maxC - allMins[j]) / hash_slices_size[j]);
920 
921  part_t hashIndexSize = hashIndices->size();
922 
923  for (part_t k = minHashIndex; k <= maxHashIndex; ++k ){
924 
925  for (part_t i = 0; i < hashIndexSize; ++i){
926  resultHashIndices->push_back(hashIndices[i] + k * hash_scales[j]);
927  }
928  }
929  tmp_swap = hashIndices;
930  hashIndices = resultHashIndices;
931  resultHashIndices = tmp_swap;
932  }
933 
934  part_t hashIndexSize = hashIndices->size();
935  for (part_t j = 0; j < hashIndexSize; ++j){
936  hash[(*hashIndices)[j]].push_back(i);
937  }
938  cpb[i].hash_indices = (*hashIndices);
939  }
940  delete hashIndices;
941  delete resultHashIndices;
942  }
943 
944  void find_neighbors(){
945 
946  }
947 
948 
949 };
950 
951 */
952 } // namespace Zoltan2
953 
954 #endif
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
~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...
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.
#define epsilon
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.