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 
393  void writeGnuPlot(std::ofstream &file,std::ofstream &mm){
394  int numCorners = (int(1)<<dim);
395  scalar_t *corner1 = new scalar_t [dim];
396  scalar_t *corner2 = new scalar_t [dim];
397 
398  for (int i = 0; i < dim; ++i){
399  /*
400  if (-maxScalar == lmins[i]){
401  if (lmaxs[i] > 0){
402  lmins[i] = lmaxs[i] / 2;
403  }
404  else{
405  lmins[i] = lmaxs[i] * 2;
406  }
407  }
408  */
409  //std::cout << lmins[i] << " ";
410  mm << lmins[i] << " ";
411  }
412  //std::cout << std::endl;
413  mm << std::endl;
414  for (int i = 0; i < dim; ++i){
415 
416  /*
417  if (maxScalar == lmaxs[i]){
418  if (lmins[i] < 0){
419  lmaxs[i] = lmins[i] / 2;
420  }
421  else{
422  lmaxs[i] = lmins[i] * 2;
423  }
424  }
425  */
426 
427  //std::cout << lmaxs[i] << " ";
428  mm << lmaxs[i] << " ";
429  }
430  //std::cout << std::endl;
431  mm << std::endl;
432 
433  for (int j = 0; j < numCorners; ++j){
434  std::vector <int> neighborCorners;
435  for (int i = 0; i < dim; ++i){
436  if(int(j & (int(1)<<i)) == 0){
437  corner1[i] = lmins[i];
438  }
439  else {
440  corner1[i] = lmaxs[i];
441  }
442  if (j % (int(1)<<(i + 1)) >= (int(1)<<i)){
443  int c1 = j - (int(1)<<i);
444 
445  if (c1 > 0) {
446  neighborCorners.push_back(c1);
447  }
448  }
449  else {
450 
451  int c1 = j + (int(1)<<i);
452  if (c1 < (int(1) << dim)) {
453  neighborCorners.push_back(c1);
454  }
455  }
456  }
457  //std::cout << "me:" << j << " nc:" << int (neighborCorners.size()) << std::endl;
458  for (int m = 0; m < int (neighborCorners.size()); ++m){
459 
460  int n = neighborCorners[m];
461  //std::cout << "me:" << j << " n:" << n << std::endl;
462  for (int i = 0; i < dim; ++i){
463  if(int(n & (int(1)<<i)) == 0){
464  corner2[i] = lmins[i];
465  }
466  else {
467  corner2[i] = lmaxs[i];
468  }
469  }
470 
471  std::string arrowline = "set arrow from ";
472  for (int i = 0; i < dim - 1; ++i){
473  arrowline +=
474  Teuchos::toString<scalar_t>(corner1[i]) + ",";
475  }
476  arrowline +=
477  Teuchos::toString<scalar_t>(corner1[dim -1]) + " to ";
478 
479  for (int i = 0; i < dim - 1; ++i){
480  arrowline +=
481  Teuchos::toString<scalar_t>(corner2[i]) + ",";
482  }
483  arrowline +=
484  Teuchos::toString<scalar_t>(corner2[dim -1]) +
485  " nohead\n";
486 
487  file << arrowline;
488  }
489  }
490  delete []corner1;
491  delete []corner2;
492  }
493 
494 
495 };
496 
497 
501 template <typename scalar_t, typename part_t>
502 class GridHash{
503 private:
504 
505  const RCP < std::vector <Zoltan2::coordinateModelPartBox <scalar_t, part_t> > > pBoxes;
506 
507  //minimum of the maximum box boundaries
508  scalar_t *minMaxBoundaries;
509  //maximum of the minimum box boundaries
510  scalar_t *maxMinBoundaries;
511  //the size of each slice along dimensions
512  scalar_t *sliceSizes;
513  part_t nTasks;
514  int dim;
515  //the number of slices per dimension
516  part_t numSlicePerDim;
517  //the number of grids - buckets
518  part_t numGrids;
519  //hash vector
520  std::vector <std::vector <part_t> > grids;
521  //result communication graph.
522  ArrayRCP <part_t> comXAdj;
523  ArrayRCP <part_t> comAdj;
524 public:
525 
529  GridHash(const RCP < std::vector <Zoltan2::coordinateModelPartBox <scalar_t, part_t> > > &pBoxes_,
530  part_t ntasks_, int dim_):
531  pBoxes(pBoxes_),
532  minMaxBoundaries(0),
533  maxMinBoundaries(0), sliceSizes(0),
534  nTasks(ntasks_),
535  dim(dim_),
536  numSlicePerDim(part_t(pow(double(ntasks_), 1.0 / dim))),
537  numGrids(0),
538  grids(),
539  comXAdj(), comAdj()
540  {
541 
542  minMaxBoundaries = new scalar_t[dim];
543  maxMinBoundaries = new scalar_t[dim];
544  sliceSizes = new scalar_t[dim];
545  //calculate the number of slices in each dimension.
546  numSlicePerDim /= 2;
547  if (numSlicePerDim == 0) numSlicePerDim = 1;
548 
549  numGrids = part_t(pow(float(numSlicePerDim), int(dim)));
550 
551  //allocate memory for buckets.
552  std::vector <std::vector <part_t> > grids_ (numGrids);
553  this->grids = grids_;
554  //get the boundaries of buckets.
555  this->getMinMaxBoundaries();
556  //insert boxes to buckets
557  this->insertToHash();
558  //calculate the neighbors for each bucket.
559  part_t nCount = this->calculateNeighbors();
560 
561  //allocate memory for communication graph
562  ArrayRCP <part_t> tmpComXadj(ntasks_+1);
563  ArrayRCP <part_t> tmpComAdj(nCount);
564  comXAdj = tmpComXadj;
565  comAdj = tmpComAdj;
566  //fill communication graph
567  this->fillAdjArrays();
568  }
569 
570 
575  delete []minMaxBoundaries;
576  delete []maxMinBoundaries;
577  delete []sliceSizes;
578  }
579 
584 
585  part_t adjIndex = 0;
586 
587  comXAdj[0] = 0;
588  for(part_t i = 0; i < this->nTasks; ++i){
589  std::set<part_t> *neigbors = (*pBoxes)[i].getNeighbors();
590 
591  part_t s = neigbors->size();
592 
593  comXAdj[i+1] = comXAdj[i] + s;
594  typedef typename std::set<part_t> mySet;
595  typedef typename mySet::iterator myIT;
596  myIT it;
597  for (it=neigbors->begin(); it!=neigbors->end(); ++it)
598 
599  comAdj[adjIndex++] = *it;
600  //TODO not needed anymore.
601  neigbors->clear();
602  }
603  }
604 
605 
606 
611  ArrayRCP <part_t> &comXAdj_,
612  ArrayRCP <part_t> &comAdj_){
613  comXAdj_ = this->comXAdj;
614  comAdj_ = this->comAdj;
615  }
616 
621  part_t nCount = 0;
622  for(part_t i = 0; i < this->nTasks; ++i){
623  std::vector <part_t> *gridIndices =(*pBoxes)[i].getGridIndices();
624  part_t gridCount = gridIndices->size();
625 
626  for (part_t j = 0; j < gridCount; ++j){
627  part_t grid = (*gridIndices)[j];
628  part_t boxCount = grids[grid].size();
629  for (part_t k = 0; k < boxCount; ++k){
630  part_t boxIndex = grids[grid][k];
631  if (boxIndex > i){
632  if((!(*pBoxes)[i].isAlreadyNeighbor(boxIndex))&& (*pBoxes)[i].isNeighborWith((*pBoxes)[boxIndex])){
633  //cout << "i:" << i << " n:" << boxIndex << " are neighbors."<< endl;
634  (*pBoxes)[i].addNeighbor(boxIndex);
635  (*pBoxes)[boxIndex].addNeighbor(i);
636  nCount += 2;
637  }
638  }
639  }
640  }
641  }
642 
643  return nCount;
644  }
645 
649  void insertToHash(){
650 
651  //cout << "ntasks:" << this->nTasks << endl;
652  for(part_t i = 0; i < this->nTasks; ++i){
653  (*pBoxes)[i].setMinMaxHashIndices(minMaxBoundaries, sliceSizes, numSlicePerDim);
654 
655  std::vector <part_t> *gridIndices =(*pBoxes)[i].getGridIndices();
656 
657  part_t gridCount = gridIndices->size();
658  //cout << "i:" << i << " gridsize:" << gridCount << endl;
659  for (part_t j = 0; j < gridCount; ++j){
660  part_t grid = (*gridIndices)[j];
661 
662  //cout << "i:" << i << " is being inserted to:" << grid << endl;
663  (grids)[grid].push_back(i);
664  }
665  }
666 
667 
668 /*
669  for(part_t i = 0; i < grids.size(); ++i){
670  cout << "grid:" << i << " gridsuze:" << (grids)[i].size() << " elements:";
671  for(part_t j = 0; j < (grids)[i].size(); ++j){
672  cout <<(grids)[i][j] << " ";
673  }
674  cout << endl;
675 
676  }
677 */
678  }
679 
684  scalar_t *mins = (*pBoxes)[0].getlmins();
685  scalar_t *maxs = (*pBoxes)[0].getlmaxs();
686 
687  for (int j = 0; j < dim; ++j){
688  minMaxBoundaries[j] = maxs[j];
689  maxMinBoundaries[j] = mins[j];
690  }
691 
692  for (part_t i = 1; i < nTasks; ++i){
693 
694  mins = (*pBoxes)[i].getlmins();
695  maxs = (*pBoxes)[i].getlmaxs();
696 
697  for (int j = 0; j < dim; ++j){
698 
699  if (minMaxBoundaries[j] > maxs[j]){
700  minMaxBoundaries[j] = maxs[j];
701  }
702  if (maxMinBoundaries[j] < mins[j]){
703  maxMinBoundaries[j] = mins[j];
704  }
705  }
706  }
707 
708 
709  for (int j = 0; j < dim; ++j){
710  sliceSizes[j] = (maxMinBoundaries[j] - minMaxBoundaries[j]) / numSlicePerDim;
711  if (sliceSizes[j] < 0) sliceSizes[j] = 0;
712  /*
713  cout << "dim:" << j <<
714  " minMax:" << minMaxBoundaries[j] <<
715  " maxMin:" << maxMinBoundaries[j] <<
716  " sliceSizes:" << sliceSizes[j] << endl;
717  */
718  }
719  }
720 };
721 /*
722 template <typename scalar_t,typename part_t>
723 class coordinatePartBox{
724 public:
725  part_t pID;
726  int dim;
727  int numCorners;
728  scalar_t **corners;
729  scalar_t *lmins, *gmins;
730  scalar_t *lmaxs, *gmaxs;
731  scalar_t maxScalar;
732  std::vector <part_t> hash_indices;
733  coordinatePartBox(part_t pid, int dim_, scalar_t *lMins, scalar_t *gMins,
734  scalar_t *lMaxs, scalar_t *gMaxs):
735  pID(pid),
736  dim(dim_),
737  numCorners(int(pow(2, dim_))),
738  corners(0),
739  lmins(lMins), gmins(gMins), lmaxs(lMaxs), gmaxs(gMaxs),
740  maxScalar (std::numeric_limits<scalar_t>::max()){
741  this->corners = new scalar_t *[dim];
742  for (int i = 0; i < dim; ++i){
743  this->corners[i] = new scalar_t[this->numCorners];
744  lmins[i] = this->maxScalar;
745  lmaxs[i] = -this->maxScalar;
746  }
747 
748 
749  for (int j = 0; j < this->numCorners; ++j){
750  for (int i = 0; i < dim; ++i){
751  std::cout << "j:" << j << " i:" << i << " 2^i:" << pow(2,i) << " and:" << int(j & int(pow(2,i))) << std::endl;
752  if(int(j & int(pow(2,i))) == 0){
753  corners[i][j] = gmins[i];
754  }
755  else {
756  corners[i][j] = gmaxs[i];
757  }
758 
759  }
760  }
761  }
762 
763 };
764 
765 template <typename Adapter, typename part_t>
766 class CoordinateCommGraph{
767 private:
768 
769  typedef typename Adapter::lno_t lno_t;
770  typedef typename Adapter::gno_t gno_t;
771  typedef typename Adapter::scalar_t scalar_t;
772 
773  const Environment *env;
774  const Teuchos::Comm<int> *comm;
775  const Zoltan2::CoordinateModel<typename Adapter::base_adapter_t> *coords;
776  const Zoltan2::PartitioningSolution<Adapter> *soln;
777  std::vector<coordinatePartBox, part_t> cpb;
778  int coordDim;
779  part_t numParts;
780 
781 
782 public:
783 
784  CoordinateCommGraph(
785  const Environment *env_,
786  const Teuchos::Comm<int> *comm_,
787  const Zoltan2::CoordinateModel<typename Adapter::base_adapter_t> *coords_,
788  const Zoltan2::PartitioningSolution<Adapter> *soln_
789  ):
790  env(env_),
791  comm(comm_),
792  coords(coords_),
793  soln(soln_),
794  coordDim (coords_->getCoordinateDim()),
795  numParts (this->soln->getActualGlobalNumberOfParts())
796  {
797  this->create_part_boxes();
798  this->hash_part_boxes();
799  this->find_neighbors();
800  }
801 
802  void create_part_boxes(){
803 
804 
805  size_t allocSize = numParts * coordDim;
806  scalar_t *lmins = new scalar_t [allocSize];
807  scalar_t *gmins = new scalar_t [allocSize];
808  scalar_t *lmaxs = new scalar_t [allocSize];
809  scalar_t *gmaxs = new scalar_t [allocSize];
810 
811  for(part_t i = 0; i < numParts; ++i){
812  coordinatePartBox tmp(
813  i,
814  this->coordDim,
815  lmins + i * coordDim,
816  gmins + i * coordDim,
817  lmaxs + i * coordDim,
818  gmaxs + i * coordDim
819  );
820  cpb.push_back(tmp);
821  }
822 
823  typedef StridedData<lno_t, scalar_t> input_t;
824  Teuchos::ArrayView<const gno_t> gnos;
825  Teuchos::ArrayView<input_t> xyz;
826  Teuchos::ArrayView<input_t> wgts;
827  coords->getCoordinates(gnos, xyz, wgts);
828 
829  //local and global num coordinates.
830  lno_t numLocalCoords = coords->getLocalNumCoordinates();
831 
832  scalar_t **pqJagged_coordinates = new scalar_t *[coordDim];
833 
834  for (int dim=0; dim < coordDim; dim++){
835  Teuchos::ArrayRCP<const scalar_t> ar;
836  xyz[dim].getInputArray(ar);
837  //pqJagged coordinate values assignment
838  pqJagged_coordinates[dim] = (scalar_t *)ar.getRawPtr();
839  }
840 
841  part_t *sol_part = soln->getPartList();
842  for(lno_t i = 0; i < numLocalCoords; ++i){
843  part_t p = sol_part[i];
844  cpb[p].updateMinMax(pqJagged_coordinates, i);
845  }
846  delete []pqJagged_coordinates;
847 
848 
849  reduceAll<int, gno_t>(*comm, Teuchos::REDUCE_MIN,
850  dim * numParts, lmins, gmins
851  );
852  reduceAll<int, gno_t>(*comm, Teuchos::REDUCE_MAX,
853  dim * numParts, lmaxs, gmaxs
854  );
855  }
856 
857  void hash_part_boxes (){
858  part_t pSingleDim = pow(double(numParts), double(1.0 / coordDim));
859  if (pSingleDim == 0) pSingleDim = 1;
860  std::vector < std::vector <part_t> > hash
861  (
862  part_t ( pow ( part_t (pSingleDim),
863  part_t(coordDim)
864  )
865  )
866  );
867 
868  //calculate the corners of the dataset.
869  scalar_t *allMins = new scalar_t [coordDim];
870  scalar_t *allMaxs = new scalar_t [coordDim];
871  part_t *hash_scales= new scalar_t [coordDim];
872 
873  for (int j = 0; j < coordDim; ++j){
874  allMins[j] = cpb[0].gmins[j];
875  allMaxs[j] = cpb[0].gmaxs[j];
876  hash_scales[j] = part_t ( pow ( part_t (pSingleDim), part_t(coordDim - j - 1)));
877  }
878 
879  for (part_t i = 1; i < numParts; ++i){
880  for (int j = 0; j < coordDim; ++j){
881  scalar_t minC = cpb[i].gmins[i];
882  scalar_t maxC = cpb[i].gmaxs[i];
883  if (minC < allMins[j]) allMins[j] = minC;
884  if (maxC > allMaxs[j]) allMaxs[j] = maxC;
885  }
886  }
887 
888  //get size of each hash for each dimension
889  scalar_t *hash_slices_size = new scalar_t [coordDim];
890  for (int j = 0; j < coordDim; ++j){
891  hash_slices_size[j] = (allMaxs[j] - allMins[j]) / pSingleDim;
892 
893  }
894 
895  delete []allMaxs;
896  delete []allMins;
897 
898 
899 
900  std::vector <part_t> *hashIndices = new std::vector <part_t>();
901  std::vector <part_t> *resultHashIndices = new std::vector <part_t>();
902  std::vector <part_t> *tmp_swap;
903  for (part_t i = 0; i < numParts; ++i){
904  hashIndices->clear();
905  resultHashIndices->clear();
906 
907  hashIndices->push_back(0);
908 
909  for (int j = 0; j < coordDim; ++j){
910 
911  scalar_t minC = cpb[i].gmins[i];
912  scalar_t maxC = cpb[i].gmaxs[i];
913  part_t minHashIndex = part_t ((minC - allMins[j]) / hash_slices_size[j]);
914  part_t maxHashIndex = part_t ((maxC - allMins[j]) / hash_slices_size[j]);
915 
916  part_t hashIndexSize = hashIndices->size();
917 
918  for (part_t k = minHashIndex; k <= maxHashIndex; ++k ){
919 
920  for (part_t i = 0; i < hashIndexSize; ++i){
921  resultHashIndices->push_back(hashIndices[i] + k * hash_scales[j]);
922  }
923  }
924  tmp_swap = hashIndices;
925  hashIndices = resultHashIndices;
926  resultHashIndices = tmp_swap;
927  }
928 
929  part_t hashIndexSize = hashIndices->size();
930  for (part_t j = 0; j < hashIndexSize; ++j){
931  hash[(*hashIndices)[j]].push_back(i);
932  }
933  cpb[i].hash_indices = (*hashIndices);
934  }
935  delete hashIndices;
936  delete resultHashIndices;
937  }
938 
939  void find_neighbors(){
940 
941  }
942 
943 
944 };
945 
946 */
947 } // namespace Zoltan2
948 
949 #endif
int getDim() const
function to set the dimension
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.
~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.
void fillAdjArrays()
GridHash Class, Function to fill adj arrays.
bool boxesOverlap(int cdim, scalar_t *lower, scalar_t *upper) const
function to test whether this box overlaps a given box
std::vector< part_t > * getGridIndices()
function to get the indices of the buckets that the part is inserted to
coordinateModelPartBox Class, represents the boundaries of the box which is a result of a geometric p...
scalar_t * getlmaxs() const
function to get maximum values along all dimensions
bool pointInBox(int pointdim, scalar_t *point) const
function to test whether a point is in the box
void computeCentroid(scalar_t *centroid) const
compute the centroid of the box
bool isNeighborWith(const coordinateModelPartBox< scalar_t, part_t > &other) const
function to check if two boxes are neighbors.
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.
GridHash(const RCP< std::vector< Zoltan2::coordinateModelPartBox< scalar_t, part_t > > > &pBoxes_, part_t ntasks_, int dim_)
GridHash Class, Constructor.
part_t getpId() const
function to get the part id
part_t calculateNeighbors()
GridHash Class, For each box compares the adjacency against the boxes that are in the same buckets...
void getMinMaxBoundaries()
GridHash Class, calculates the minimum of maximum box boundaries, and maxium of minimum box boundarie...
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.
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.
scalar_t * getlmins() const
function to get minimum values along all dimensions
coordinateModelPartBox(part_t pid, int dim_, scalar_t *lmi, scalar_t *lma)
Constructor deep copy of the maximum and minimum boundaries.