Zoltan2
Loading...
Searching...
No Matches
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"
61
62namespace Zoltan2{
63
64
65#define Z2_ABS(x) ((x) >= 0 ? (x) : -(x))
66
71
72public:
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 }
181 return this->lmins;
182 }
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 }
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
485private:
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
511private:
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;
535public:
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
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/*
733template <typename coord_t,typename part_t>
734class coordinatePartBox{
735public:
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
776template <typename Adapter, typename part_t>
777class CoordinateCommGraph{
778private:
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
793public:
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
Traits for application input objects.
GridHash Class, Hashing Class for part boxes.
void fillAdjArrays()
GridHash Class, Function to fill adj arrays.
part_t calculateNeighbors()
GridHash Class, For each box compares the adjacency against the boxes that are in the same buckets.
~GridHash()
GridHash Class, Destructor.
GridHash(const RCP< std::vector< Zoltan2::coordinateModelPartBox > > &pBoxes_, part_t ntasks_, int dim_)
GridHash Class, Constructor.
void getMinMaxBoundaries()
GridHash Class, calculates the minimum of maximum box boundaries, and maxium of minimum box boundarie...
void insertToHash()
GridHash Class, For each box calculates the buckets which it should be inserted to.
void getAdjArrays(ArrayRCP< part_t > &comXAdj_, ArrayRCP< part_t > &comAdj_)
GridHash Class, returns the adj arrays.
coordinateModelPartBox Class, represents the boundaries of the box which is a result of a geometric p...
coordinateModelPartBox(part_t pid, int dim_, scalar_t *lmi, scalar_t *lma)
Constructor deep copy of the maximum and minimum boundaries.
void setMinMaxHashIndices(coord_t *minMaxBoundaries, coord_t *sliceSizes, part_t numSlicePerDim)
function to obtain the min and max hash values along all dimensions.
bool boxesOverlap(int cdim, scalar_t *lower, scalar_t *upper) const
function to test whether this box overlaps a given box
bool isNeighborWith(const coordinateModelPartBox &other) const
function to check if two boxes are neighbors.
bool isAlreadyNeighbor(part_t nIndex)
function to check if a given part is already in the neighbor list.
coord_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
std::vector< part_t > * getGridIndices()
function to get the indices of the buckets that the part is inserted to
void computeCentroid(coord_t *centroid) const
compute the centroid of the box
void updateMinMax(scalar_t newBoundary, int isMax, int dimInd)
function to update the boundary of the box.
part_t getpId() const
function to get the part id
coordinateModelPartBox(part_t pid, int dim_)
Constructor.
void setpId(part_t pid)
function to set the part id
std::set< part_t > * getNeighbors()
function to get the indices of the neighboring parts.
void addNeighbor(part_t nIndex)
function to add a new neighbor to the neighbor list.
coordinateModelPartBox(const coordinateModelPartBox &other)
Copy Constructor deep copy of the maximum and minimum boundaries.
int getDim() const
function to set the dimension
void print()
function to print the boundaries.
void writeGnuPlot(std::ofstream &file, std::ofstream &mm)
function for visualization.
coord_t * getlmins() const
function to get minimum values along all dimensions
Created by mbenlioglu on Aug 31, 2020.
SparseMatrixAdapter_t::part_t part_t