Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 #include "Zoltan2_InputTraits.hpp"
61 
62 namespace Zoltan2{
63 
64 
65 #define Z2_ABS(x) ((x) >= 0 ? (x) : -(x))
66 
71 
72 public:
73  typedef double coord_t;
75 
79  pID(pid),
80  dim(dim_),
81  lmins(0), lmaxs(0),
82  maxScalar (std::numeric_limits<coord_t>::max()),
83  _EPSILON(std::numeric_limits<coord_t>::epsilon()),
84  minHashIndices(0),
85  maxHashIndices(0),
86  gridIndices(0), neighbors()
87  {
88  lmins = new coord_t [dim];
89  lmaxs = new coord_t [dim];
90 
91  minHashIndices = new part_t [dim];
92  maxHashIndices = new part_t [dim];
93  gridIndices = new std::vector <part_t> ();
94  for (int i = 0; i < dim; ++i){
95  lmins[i] = -this->maxScalar;
96  lmaxs[i] = this->maxScalar;
97  }
98  }
102  template <typename scalar_t>
103  coordinateModelPartBox(part_t pid, int dim_, scalar_t *lmi, scalar_t *lma):
104  pID(pid),
105  dim(dim_),
106  lmins(0), lmaxs(0),
107  maxScalar (std::numeric_limits<coord_t>::max()),
108  _EPSILON(std::numeric_limits<coord_t>::epsilon()),
109  minHashIndices(0),
110  maxHashIndices(0),
111  gridIndices(0), neighbors()
112  {
113  lmins = new coord_t [dim];
114  lmaxs = new coord_t [dim];
115  minHashIndices = new part_t [dim];
116  maxHashIndices = new part_t [dim];
117  gridIndices = new std::vector <part_t> ();
118  for (int i = 0; i < dim; ++i){
119  lmins[i] = static_cast<coord_t>(lmi[i]);
120  lmaxs[i] = static_cast<coord_t>(lma[i]);
121  }
122  }
123 
124 
129  pID(other.getpId()),
130  dim(other.getDim()),
131  lmins(0), lmaxs(0),
132  maxScalar (std::numeric_limits<coord_t>::max()),
133  _EPSILON(std::numeric_limits<coord_t>::epsilon()),
134  minHashIndices(0),
135  maxHashIndices(0),
136  gridIndices(0), neighbors()
137  {
138 
139  lmins = new coord_t [dim];
140  lmaxs = new coord_t [dim];
141  minHashIndices = new part_t [dim];
142  maxHashIndices = new part_t [dim];
143  gridIndices = new std::vector <part_t> ();
144  coord_t *othermins = other.getlmins();
145  coord_t *othermaxs = other.getlmaxs();
146  for (int i = 0; i < dim; ++i){
147  lmins[i] = othermins[i];
148  lmaxs[i] = othermaxs[i];
149  }
150  }
154  delete []this->lmins;
155  delete [] this->lmaxs;
156  delete []this->minHashIndices;
157  delete [] this->maxHashIndices;
158  delete gridIndices;
159  }
160 
163  void setpId(part_t pid){
164  this->pID = pid;
165  }
168  part_t getpId() const{
169  return this->pID;
170  }
171 
172 
175  int getDim()const{
176  return this->dim;
177  }
180  coord_t * getlmins()const{
181  return this->lmins;
182  }
185  coord_t * getlmaxs()const{
186  return this->lmaxs;
187  }
190  void computeCentroid(coord_t *centroid)const {
191  for (int i = 0; i < this->dim; i++)
192  centroid[i] = 0.5 * (this->lmaxs[i] + this->lmins[i]);
193  }
194 
198  std::vector <part_t> * getGridIndices () {
199  return this->gridIndices;
200  }
201 
204  std::set<part_t> *getNeighbors() {
205  return &(this->neighbors);
206  }
207 
210  template <typename scalar_t>
211  bool pointInBox(int pointdim, scalar_t *point) const {
212  if (pointdim != this->dim)
213  throw std::logic_error("dim of point must match dim of box");
214  for (int i = 0; i < pointdim; i++) {
215  if (static_cast<coord_t>(point[i]) < this->lmins[i]) return false;
216  if (static_cast<coord_t>(point[i]) > this->lmaxs[i]) return false;
217  }
218  return true;
219  }
220 
223  template <typename scalar_t>
224  bool boxesOverlap(int cdim, scalar_t *lower, scalar_t *upper) const {
225  if (cdim != this->dim)
226  throw std::logic_error("dim of given box must match dim of box");
227 
228  // Check for at least partial overlap
229  bool found = true;
230  for (int i = 0; i < cdim; i++) {
231  if (!((static_cast<coord_t>(lower[i]) >= this->lmins[i] &&
232  static_cast<coord_t>(lower[i]) <= this->lmaxs[i])
233  // lower i-coordinate in the box
234  || (static_cast<coord_t>(upper[i]) >= this->lmins[i] &&
235  static_cast<coord_t>(upper[i]) <= this->lmaxs[i])
236  // upper i-coordinate in the box
237  || (static_cast<coord_t>(lower[i]) < this->lmins[i] &&
238  static_cast<coord_t>(upper[i]) > this->lmaxs[i]))) {
239  // i-coordinates straddle the box
240  found = false;
241  break;
242  }
243  }
244  return found;
245  }
246 
250  const coordinateModelPartBox &other) const{
251 
252  coord_t *omins = other.getlmins();
253  coord_t *omaxs = other.getlmaxs();
254 
255  int equality = 0;
256  for (int i = 0; i < dim; ++i){
257 
258  if (omins[i] - this->lmaxs[i] > _EPSILON ||
259  this->lmins[i] - omaxs[i] > _EPSILON ) {
260  return false;
261  }
262  else if (Z2_ABS(omins[i] - this->lmaxs[i]) < _EPSILON ||
263  Z2_ABS(this->lmins[i] - omaxs[i]) < _EPSILON ){
264  if (++equality > 1){
265  return false;
266  }
267  }
268  }
269  if (equality == 1) {
270  return true;
271  }
272  else {
273  std::cout << "something is wrong: equality:"
274  << equality << std::endl;
275  return false;
276  }
277  }
278 
279 
282  void addNeighbor(part_t nIndex){
283  neighbors.insert(nIndex);
284  }
287  bool isAlreadyNeighbor(part_t nIndex){
288 
289  if (neighbors.end() != neighbors.find(nIndex)){
290  return true;
291  }
292  return false;
293 
294  }
295 
296 
300  coord_t *minMaxBoundaries,
301  coord_t *sliceSizes,
302  part_t numSlicePerDim
303  ){
304  for (int j = 0; j < dim; ++j){
305  coord_t distance = (lmins[j] - minMaxBoundaries[j]);
306  part_t minInd = 0;
307  if (distance > _EPSILON && sliceSizes[j] > _EPSILON){
308  minInd = static_cast<part_t>(floor((lmins[j] - minMaxBoundaries[j])/ sliceSizes[j]));
309  }
310 
311  if(minInd >= numSlicePerDim){
312  minInd = numSlicePerDim - 1;
313  }
314  part_t maxInd = 0;
315  distance = (lmaxs[j] - minMaxBoundaries[j]);
316  if (distance > _EPSILON && sliceSizes[j] > _EPSILON){
317  maxInd = static_cast<part_t>(ceil((lmaxs[j] - minMaxBoundaries[j])/ sliceSizes[j]));
318  }
319  if(maxInd >= numSlicePerDim){
320  maxInd = numSlicePerDim - 1;
321  }
322 
323  //cout << "j:" << j << " lmins:" << lmins[j] << " lmaxs:" << lmaxs[j] << endl;
324  //cout << "j:" << j << " min:" << minInd << " max:" << maxInd << endl;
325  minHashIndices[j] = minInd;
326  maxHashIndices[j] = maxInd;
327  }
328  std::vector <part_t> *in = new std::vector <part_t> ();
329  in->push_back(0);
330  std::vector <part_t> *out = new std::vector <part_t> ();
331 
332  for (int j = 0; j < dim; ++j){
333 
334  part_t minInd = minHashIndices[j];
335  part_t maxInd = maxHashIndices[j];
336 
337 
338  part_t pScale = part_t(pow (float(numSlicePerDim), int(dim - j -1)));
339 
340  part_t inSize = in->size();
341 
342  for (part_t k = minInd; k <= maxInd; ++k){
343  for (part_t i = 0; i < inSize; ++i){
344  out->push_back((*in)[i] + k * pScale);
345  }
346  }
347  in->clear();
348  std::vector <part_t> *tmp = in;
349  in= out;
350  out= tmp;
351  }
352 
353  std::vector <part_t> *tmp = in;
354  in = gridIndices;
355  gridIndices = tmp;
356 
357 
358  delete in;
359  delete out;
360  }
361 
364  void print(){
365  for(int i = 0; i < this->dim; ++i){
366  std::cout << "\tbox:" << this->pID << " dim:" << i << " min:" << lmins[i] << " max:" << lmaxs[i] << std::endl;
367  }
368  }
369 
372  template <typename scalar_t>
373  void updateMinMax (scalar_t newBoundary, int isMax, int dimInd){
374  if (isMax){
375  lmaxs[dimInd] = static_cast<coord_t>(newBoundary);
376  }
377  else {
378  lmins[dimInd] = static_cast<coord_t>(newBoundary);
379  }
380  }
381 
384  void writeGnuPlot(std::ofstream &file,std::ofstream &mm){
385  int numCorners = (int(1)<<dim);
386  coord_t *corner1 = new coord_t [dim];
387  coord_t *corner2 = new coord_t [dim];
388 
389  for (int i = 0; i < dim; ++i){
390  /*
391  if (-maxScalar == lmins[i]){
392  if (lmaxs[i] > 0){
393  lmins[i] = lmaxs[i] / 2;
394  }
395  else{
396  lmins[i] = lmaxs[i] * 2;
397  }
398  }
399  */
400  //std::cout << lmins[i] << " ";
401  mm << lmins[i] << " ";
402  }
403  //std::cout << std::endl;
404  mm << std::endl;
405  for (int i = 0; i < dim; ++i){
406 
407  /*
408  if (maxScalar == lmaxs[i]){
409  if (lmins[i] < 0){
410  lmaxs[i] = lmins[i] / 2;
411  }
412  else{
413  lmaxs[i] = lmins[i] * 2;
414  }
415  }
416  */
417 
418  //std::cout << lmaxs[i] << " ";
419  mm << lmaxs[i] << " ";
420  }
421  //std::cout << std::endl;
422  mm << std::endl;
423 
424  for (int j = 0; j < numCorners; ++j){
425  std::vector <int> neighborCorners;
426  for (int i = 0; i < dim; ++i){
427  if(int(j & (int(1)<<i)) == 0){
428  corner1[i] = lmins[i];
429  }
430  else {
431  corner1[i] = lmaxs[i];
432  }
433  if (j % (int(1)<<(i + 1)) >= (int(1)<<i)){
434  int c1 = j - (int(1)<<i);
435 
436  if (c1 > 0) {
437  neighborCorners.push_back(c1);
438  }
439  }
440  else {
441 
442  int c1 = j + (int(1)<<i);
443  if (c1 < (int(1) << dim)) {
444  neighborCorners.push_back(c1);
445  }
446  }
447  }
448  //std::cout << "me:" << j << " nc:" << int (neighborCorners.size()) << std::endl;
449  for (int m = 0; m < int (neighborCorners.size()); ++m){
450 
451  int n = neighborCorners[m];
452  //std::cout << "me:" << j << " n:" << n << std::endl;
453  for (int i = 0; i < dim; ++i){
454  if(int(n & (int(1)<<i)) == 0){
455  corner2[i] = lmins[i];
456  }
457  else {
458  corner2[i] = lmaxs[i];
459  }
460  }
461 
462  std::string arrowline = "set arrow from ";
463  for (int i = 0; i < dim - 1; ++i){
464  arrowline +=
465  Teuchos::toString<coord_t>(corner1[i]) + ",";
466  }
467  arrowline +=
468  Teuchos::toString<coord_t>(corner1[dim -1]) + " to ";
469 
470  for (int i = 0; i < dim - 1; ++i){
471  arrowline +=
472  Teuchos::toString<coord_t>(corner2[i]) + ",";
473  }
474  arrowline +=
475  Teuchos::toString<coord_t>(corner2[dim -1]) +
476  " nohead\n";
477 
478  file << arrowline;
479  }
480  }
481  delete []corner1;
482  delete []corner2;
483  }
484 
485 private:
486  part_t pID; //part Id
487  int dim; //dimension of the box
488  coord_t *lmins; //minimum boundaries of the box along all dimensions.
489  coord_t *lmaxs; //maximum boundaries of the box along all dimensions.
490  coord_t maxScalar;
491  coord_t _EPSILON;
492 
493  //to calculate the neighbors of the box and avoid the p^2 comparisons,
494  //we use hashing. A box can be put into multiple hash buckets.
495  //the following 2 variable holds the minimum and maximum of the
496  //hash values along all dimensions.
497  part_t *minHashIndices;
498  part_t *maxHashIndices;
499 
500  //result hash bucket indices.
501  std::vector <part_t> *gridIndices;
502  //neighbors of the box.
503  std::set <part_t> neighbors;
504 };
505 
506 
510 class GridHash{
511 private:
512 
513  typedef typename Zoltan2::coordinateModelPartBox::coord_t coord_t;
514  typedef typename Zoltan2::coordinateModelPartBox::part_t part_t;
515 
516  const RCP<std::vector<Zoltan2::coordinateModelPartBox> > pBoxes;
517 
518  //minimum of the maximum box boundaries
519  coord_t *minMaxBoundaries;
520  //maximum of the minimum box boundaries
521  coord_t *maxMinBoundaries;
522  //the size of each slice along dimensions
523  coord_t *sliceSizes;
524  part_t nTasks;
525  int dim;
526  //the number of slices per dimension
527  part_t numSlicePerDim;
528  //the number of grids - buckets
529  part_t numGrids;
530  //hash vector
531  std::vector <std::vector <part_t> > grids;
532  //result communication graph.
533  ArrayRCP <part_t> comXAdj;
534  ArrayRCP <part_t> comAdj;
535 public:
536 
540  GridHash(const RCP<std::vector<Zoltan2::coordinateModelPartBox> > &pBoxes_,
541  part_t ntasks_, int dim_):
542  pBoxes(pBoxes_),
543  minMaxBoundaries(0),
544  maxMinBoundaries(0), sliceSizes(0),
545  nTasks(ntasks_),
546  dim(dim_),
547  numSlicePerDim(part_t(pow(double(ntasks_), 1.0 / dim))),
548  numGrids(0),
549  grids(),
550  comXAdj(), comAdj()
551  {
552 
553  minMaxBoundaries = new coord_t[dim];
554  maxMinBoundaries = new coord_t[dim];
555  sliceSizes = new coord_t[dim];
556  //calculate the number of slices in each dimension.
557  numSlicePerDim /= 2;
558  if (numSlicePerDim == 0) numSlicePerDim = 1;
559 
560  numGrids = part_t(pow(float(numSlicePerDim), int(dim)));
561 
562  //allocate memory for buckets.
563  std::vector <std::vector <part_t> > grids_ (numGrids);
564  this->grids = grids_;
565  //get the boundaries of buckets.
566  this->getMinMaxBoundaries();
567  //insert boxes to buckets
568  this->insertToHash();
569  //calculate the neighbors for each bucket.
570  part_t nCount = this->calculateNeighbors();
571 
572  //allocate memory for communication graph
573  ArrayRCP <part_t> tmpComXadj(ntasks_+1);
574  ArrayRCP <part_t> tmpComAdj(nCount);
575  comXAdj = tmpComXadj;
576  comAdj = tmpComAdj;
577  //fill communication graph
578  this->fillAdjArrays();
579  }
580 
581 
586  delete []minMaxBoundaries;
587  delete []maxMinBoundaries;
588  delete []sliceSizes;
589  }
590 
595 
596  part_t adjIndex = 0;
597 
598  comXAdj[0] = 0;
599  for(part_t i = 0; i < this->nTasks; ++i){
600  std::set<part_t> *neigbors = (*pBoxes)[i].getNeighbors();
601 
602  part_t s = neigbors->size();
603 
604  comXAdj[i+1] = comXAdj[i] + s;
605  typedef typename std::set<part_t> mySet;
606  typedef typename mySet::iterator myIT;
607  myIT it;
608  for (it=neigbors->begin(); it!=neigbors->end(); ++it)
609 
610  comAdj[adjIndex++] = *it;
611  //TODO not needed anymore.
612  neigbors->clear();
613  }
614  }
615 
616 
617 
622  ArrayRCP <part_t> &comXAdj_,
623  ArrayRCP <part_t> &comAdj_){
624  comXAdj_ = this->comXAdj;
625  comAdj_ = this->comAdj;
626  }
627 
632  part_t nCount = 0;
633  for(part_t i = 0; i < this->nTasks; ++i){
634  std::vector <part_t> *gridIndices =(*pBoxes)[i].getGridIndices();
635  part_t gridCount = gridIndices->size();
636 
637  for (part_t j = 0; j < gridCount; ++j){
638  part_t grid = (*gridIndices)[j];
639  part_t boxCount = grids[grid].size();
640  for (part_t k = 0; k < boxCount; ++k){
641  part_t boxIndex = grids[grid][k];
642  if (boxIndex > i){
643  if((!(*pBoxes)[i].isAlreadyNeighbor(boxIndex))&& (*pBoxes)[i].isNeighborWith((*pBoxes)[boxIndex])){
644  //cout << "i:" << i << " n:" << boxIndex << " are neighbors."<< endl;
645  (*pBoxes)[i].addNeighbor(boxIndex);
646  (*pBoxes)[boxIndex].addNeighbor(i);
647  nCount += 2;
648  }
649  }
650  }
651  }
652  }
653 
654  return nCount;
655  }
656 
660  void insertToHash(){
661 
662  //cout << "ntasks:" << this->nTasks << endl;
663  for(part_t i = 0; i < this->nTasks; ++i){
664  (*pBoxes)[i].setMinMaxHashIndices(minMaxBoundaries, sliceSizes, numSlicePerDim);
665 
666  std::vector <part_t> *gridIndices =(*pBoxes)[i].getGridIndices();
667 
668  part_t gridCount = gridIndices->size();
669  //cout << "i:" << i << " gridsize:" << gridCount << endl;
670  for (part_t j = 0; j < gridCount; ++j){
671  part_t grid = (*gridIndices)[j];
672 
673  //cout << "i:" << i << " is being inserted to:" << grid << endl;
674  (grids)[grid].push_back(i);
675  }
676  }
677 
678 
679 /*
680  for(part_t i = 0; i < grids.size(); ++i){
681  cout << "grid:" << i << " gridsuze:" << (grids)[i].size() << " elements:";
682  for(part_t j = 0; j < (grids)[i].size(); ++j){
683  cout <<(grids)[i][j] << " ";
684  }
685  cout << endl;
686 
687  }
688 */
689  }
690 
695  coord_t *mins = (*pBoxes)[0].getlmins();
696  coord_t *maxs = (*pBoxes)[0].getlmaxs();
697 
698  for (int j = 0; j < dim; ++j){
699  minMaxBoundaries[j] = maxs[j];
700  maxMinBoundaries[j] = mins[j];
701  }
702 
703  for (part_t i = 1; i < nTasks; ++i){
704 
705  mins = (*pBoxes)[i].getlmins();
706  maxs = (*pBoxes)[i].getlmaxs();
707 
708  for (int j = 0; j < dim; ++j){
709 
710  if (minMaxBoundaries[j] > maxs[j]){
711  minMaxBoundaries[j] = maxs[j];
712  }
713  if (maxMinBoundaries[j] < mins[j]){
714  maxMinBoundaries[j] = mins[j];
715  }
716  }
717  }
718 
719 
720  for (int j = 0; j < dim; ++j){
721  sliceSizes[j] = (maxMinBoundaries[j] - minMaxBoundaries[j]) / numSlicePerDim;
722  if (sliceSizes[j] < 0) sliceSizes[j] = 0;
723  /*
724  cout << "dim:" << j <<
725  " minMax:" << minMaxBoundaries[j] <<
726  " maxMin:" << maxMinBoundaries[j] <<
727  " sliceSizes:" << sliceSizes[j] << endl;
728  */
729  }
730  }
731 };
732 /*
733 template <typename coord_t,typename part_t>
734 class coordinatePartBox{
735 public:
736  part_t pID;
737  int dim;
738  int numCorners;
739  coord_t **corners;
740  coord_t *lmins, *gmins;
741  coord_t *lmaxs, *gmaxs;
742  coord_t maxScalar;
743  std::vector <part_t> hash_indices;
744  coordinatePartBox(part_t pid, int dim_, coord_t *lMins, coord_t *gMins,
745  coord_t *lMaxs, coord_t *gMaxs):
746  pID(pid),
747  dim(dim_),
748  numCorners(int(pow(2, dim_))),
749  corners(0),
750  lmins(lMins), gmins(gMins), lmaxs(lMaxs), gmaxs(gMaxs),
751  maxScalar (std::numeric_limits<coord_t>::max()){
752  this->corners = new coord_t *[dim];
753  for (int i = 0; i < dim; ++i){
754  this->corners[i] = new coord_t[this->numCorners];
755  lmins[i] = this->maxScalar;
756  lmaxs[i] = -this->maxScalar;
757  }
758 
759 
760  for (int j = 0; j < this->numCorners; ++j){
761  for (int i = 0; i < dim; ++i){
762  std::cout << "j:" << j << " i:" << i << " 2^i:" << pow(2,i) << " and:" << int(j & int(pow(2,i))) << std::endl;
763  if(int(j & int(pow(2,i))) == 0){
764  corners[i][j] = gmins[i];
765  }
766  else {
767  corners[i][j] = gmaxs[i];
768  }
769 
770  }
771  }
772  }
773 
774 };
775 
776 template <typename Adapter, typename part_t>
777 class CoordinateCommGraph{
778 private:
779 
780  typedef typename Adapter::lno_t lno_t;
781  typedef typename Adapter::gno_t gno_t;
782  typedef typename Adapter::coord_t coord_t;
783 
784  const Environment *env;
785  const Teuchos::Comm<int> *comm;
786  const Zoltan2::CoordinateModel<typename Adapter::base_adapter_t> *coords;
787  const Zoltan2::PartitioningSolution<Adapter> *soln;
788  std::vector<coordinatePartBox, part_t> cpb;
789  int coordDim;
790  part_t numParts;
791 
792 
793 public:
794 
795  CoordinateCommGraph(
796  const Environment *env_,
797  const Teuchos::Comm<int> *comm_,
798  const Zoltan2::CoordinateModel<typename Adapter::base_adapter_t> *coords_,
799  const Zoltan2::PartitioningSolution<Adapter> *soln_
800  ):
801  env(env_),
802  comm(comm_),
803  coords(coords_),
804  soln(soln_),
805  coordDim (coords_->getCoordinateDim()),
806  numParts (this->soln->getActualGlobalNumberOfParts())
807  {
808  this->create_part_boxes();
809  this->hash_part_boxes();
810  this->find_neighbors();
811  }
812 
813  void create_part_boxes(){
814 
815 
816  size_t allocSize = numParts * coordDim;
817  coord_t *lmins = new coord_t [allocSize];
818  coord_t *gmins = new coord_t [allocSize];
819  coord_t *lmaxs = new coord_t [allocSize];
820  coord_t *gmaxs = new coord_t [allocSize];
821 
822  for(part_t i = 0; i < numParts; ++i){
823  coordinatePartBox tmp(
824  i,
825  this->coordDim,
826  lmins + i * coordDim,
827  gmins + i * coordDim,
828  lmaxs + i * coordDim,
829  gmaxs + i * coordDim
830  );
831  cpb.push_back(tmp);
832  }
833 
834  typedef StridedData<lno_t, coord_t> input_t;
835  Teuchos::ArrayView<const gno_t> gnos;
836  Teuchos::ArrayView<input_t> xyz;
837  Teuchos::ArrayView<input_t> wgts;
838  coords->getCoordinates(gnos, xyz, wgts);
839 
840  //local and global num coordinates.
841  lno_t numLocalCoords = coords->getLocalNumCoordinates();
842 
843  coord_t **pqJagged_coordinates = new coord_t *[coordDim];
844 
845  for (int dim=0; dim < coordDim; dim++){
846  Teuchos::ArrayRCP<const coord_t> ar;
847  xyz[dim].getInputArray(ar);
848  //pqJagged coordinate values assignment
849  pqJagged_coordinates[dim] = (coord_t *)ar.getRawPtr();
850  }
851 
852  part_t *sol_part = soln->getPartList();
853  for(lno_t i = 0; i < numLocalCoords; ++i){
854  part_t p = sol_part[i];
855  cpb[p].updateMinMax(pqJagged_coordinates, i);
856  }
857  delete []pqJagged_coordinates;
858 
859 
860  reduceAll<int, gno_t>(*comm, Teuchos::REDUCE_MIN,
861  dim * numParts, lmins, gmins
862  );
863  reduceAll<int, gno_t>(*comm, Teuchos::REDUCE_MAX,
864  dim * numParts, lmaxs, gmaxs
865  );
866  }
867 
868  void hash_part_boxes (){
869  part_t pSingleDim = pow(double(numParts), double(1.0 / coordDim));
870  if (pSingleDim == 0) pSingleDim = 1;
871  std::vector < std::vector <part_t> > hash
872  (
873  part_t ( pow ( part_t (pSingleDim),
874  part_t(coordDim)
875  )
876  )
877  );
878 
879  //calculate the corners of the dataset.
880  coord_t *allMins = new coord_t [coordDim];
881  coord_t *allMaxs = new coord_t [coordDim];
882  part_t *hash_scales= new coord_t [coordDim];
883 
884  for (int j = 0; j < coordDim; ++j){
885  allMins[j] = cpb[0].gmins[j];
886  allMaxs[j] = cpb[0].gmaxs[j];
887  hash_scales[j] = part_t ( pow ( part_t (pSingleDim), part_t(coordDim - j - 1)));
888  }
889 
890  for (part_t i = 1; i < numParts; ++i){
891  for (int j = 0; j < coordDim; ++j){
892  coord_t minC = cpb[i].gmins[i];
893  coord_t maxC = cpb[i].gmaxs[i];
894  if (minC < allMins[j]) allMins[j] = minC;
895  if (maxC > allMaxs[j]) allMaxs[j] = maxC;
896  }
897  }
898 
899  //get size of each hash for each dimension
900  coord_t *hash_slices_size = new coord_t [coordDim];
901  for (int j = 0; j < coordDim; ++j){
902  hash_slices_size[j] = (allMaxs[j] - allMins[j]) / pSingleDim;
903 
904  }
905 
906  delete []allMaxs;
907  delete []allMins;
908 
909 
910 
911  std::vector <part_t> *hashIndices = new std::vector <part_t>();
912  std::vector <part_t> *resultHashIndices = new std::vector <part_t>();
913  std::vector <part_t> *tmp_swap;
914  for (part_t i = 0; i < numParts; ++i){
915  hashIndices->clear();
916  resultHashIndices->clear();
917 
918  hashIndices->push_back(0);
919 
920  for (int j = 0; j < coordDim; ++j){
921 
922  coord_t minC = cpb[i].gmins[i];
923  coord_t maxC = cpb[i].gmaxs[i];
924  part_t minHashIndex = part_t ((minC - allMins[j]) / hash_slices_size[j]);
925  part_t maxHashIndex = part_t ((maxC - allMins[j]) / hash_slices_size[j]);
926 
927  part_t hashIndexSize = hashIndices->size();
928 
929  for (part_t k = minHashIndex; k <= maxHashIndex; ++k ){
930 
931  for (part_t i = 0; i < hashIndexSize; ++i){
932  resultHashIndices->push_back(hashIndices[i] + k * hash_scales[j]);
933  }
934  }
935  tmp_swap = hashIndices;
936  hashIndices = resultHashIndices;
937  resultHashIndices = tmp_swap;
938  }
939 
940  part_t hashIndexSize = hashIndices->size();
941  for (part_t j = 0; j < hashIndexSize; ++j){
942  hash[(*hashIndices)[j]].push_back(i);
943  }
944  cpb[i].hash_indices = (*hashIndices);
945  }
946  delete hashIndices;
947  delete resultHashIndices;
948  }
949 
950  void find_neighbors(){
951 
952  }
953 
954 
955 };
956 
957 */
958 } // namespace Zoltan2
959 
960 #endif
~GridHash()
GridHash Class, Destructor.
GridHash Class, Hashing Class for part boxes.
std::vector< part_t > * getGridIndices()
function to get the indices of the buckets that the part is inserted to
coord_t * getlmaxs() const
function to get maximum 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.
std::set< part_t > * getNeighbors()
function to get the indices of the neighboring parts.
part_t calculateNeighbors()
GridHash Class, For each box compares the adjacency against the boxes that are in the same buckets...
void getAdjArrays(ArrayRCP< part_t > &comXAdj_, ArrayRCP< part_t > &comAdj_)
GridHash Class, returns the adj arrays.
void updateMinMax(scalar_t newBoundary, int isMax, int dimInd)
function to update the boundary of the box.
coordinateModelPartBox(part_t pid, int dim_)
Constructor.
GridHash(const RCP< std::vector< Zoltan2::coordinateModelPartBox > > &pBoxes_, part_t ntasks_, int dim_)
GridHash Class, Constructor.
void insertToHash()
GridHash Class, For each box calculates the buckets which it should be inserted to.
void computeCentroid(coord_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...
SparseMatrixAdapter_t::part_t part_t
bool isNeighborWith(const coordinateModelPartBox &other) const
function to check if two boxes are neighbors.
void getMinMaxBoundaries()
GridHash Class, calculates the minimum of maximum box boundaries, and maxium of minimum box boundarie...
void setMinMaxHashIndices(coord_t *minMaxBoundaries, coord_t *sliceSizes, part_t numSlicePerDim)
function to obtain the min and max hash values along all dimensions.
void setpId(part_t pid)
function to set the part id
coordinateModelPartBox(const coordinateModelPartBox &other)
Copy Constructor deep copy of the maximum and minimum boundaries.
Traits for application input objects.
#define epsilon
coord_t * getlmins() const
function to get minimum values along all dimensions
void addNeighbor(part_t nIndex)
function to add a new neighbor to the neighbor list.
bool isAlreadyNeighbor(part_t nIndex)
function to check if a given part is already in the neighbor list.
void writeGnuPlot(std::ofstream &file, std::ofstream &mm)
function for visualization.
bool boxesOverlap(int cdim, scalar_t *lower, scalar_t *upper) const
function to test whether this box overlaps a given box
bool pointInBox(int pointdim, scalar_t *point) const
function to test whether a point is in the box
part_t getpId() const
function to get the part id
void print()
function to print the boundaries.
int getDim() const
function to set the dimension
void fillAdjArrays()
GridHash Class, Function to fill adj arrays.