1 #ifndef LIBGEODECOMP_GEOMETRY_PARTITIONS_PTSCOTCHPARTITION_H
2 #define LIBGEODECOMP_GEOMETRY_PARTITIONS_PTSCOTCHPARTITION_H
3 
4 #include <libgeodecomp/config.h>
5 #include <libgeodecomp/geometry/partitions/partition.h>
6 
7 #ifdef LIBGEODECOMP_WITH_SCOTCH
8 #ifdef LIBGEODECOMP_WITH_MPI
9 
10 #include <mpi.h>
11 #include <ptscotch.h>
12 
13 namespace LibGeoDecomp {
14 
15 template<int DIM>
16 class PTScotchPartition : public Partition<DIM>
17 {
18 public:
19     using Partition<DIM>::startOffsets;
20     using Partition<DIM>::weights;
21 
22     inline explicit PTScotchPartition(
23         const Coord<DIM>& origin = Coord<DIM>(),
24         const Coord<DIM>& dimensions = Coord<DIM>(),
25         const long& offset = 0,
26         const std::vector<std::size_t>& weights = std::vector<std::size_t>(2)) :
27         Partition<DIM>(offset, weights),
28         origin(origin),
29         dimensions(dimensions),
30         cellNbr(dimensions.prod())
31     {
32         std::vector<SCOTCH_Num> indices(cellNbr);
33         initIndices(&indices[0]);
34         regions.resize(weights.size());
35         createRegions(&indices[0]);
36     }
37 
getRegion(const std::size_t node)38     Region<DIM> getRegion(const std::size_t node) const
39     {
40         return regions[node];
41     }
42 
43 private:
44     Coord<DIM> origin;
45     Coord<DIM> dimensions;
46     SCOTCH_Num cellNbr;
47     std::vector<Region<DIM> > regions;
48 
initIndices(SCOTCH_Num * indices)49     void initIndices(SCOTCH_Num *indices)
50     {
51         SCOTCH_Arch arch;
52         SCOTCH_archInit(&arch);
53         SCOTCH_Num * velotabArch;
54         SCOTCH_Num const vertnbrArch = weights.size();
55         velotabArch = new SCOTCH_Num[weights.size()];
56         for(int i = 0; i < vertnbrArch; ++i){
57             velotabArch[i] = weights[i];
58         }
59         SCOTCH_archCmpltw(&arch, vertnbrArch, velotabArch);
60 
61         SCOTCH_Dgraph grafdat;
62         SCOTCH_dgraphInit(&grafdat, MPI_COMM_WORLD);
63 
64         SCOTCH_Num edgenbrGra = 2 * (dimensions[0] *
65                                      (dimensions[1] - 1) +
66                                      (dimensions[0] - 1) * dimensions[1]);
67         if(DIM == 3){
68             edgenbrGra = edgenbrGra
69                 * dimensions[2]
70                 + 2 * (dimensions[0] *
71                        dimensions[1] * (dimensions[2] - 1));
72         }
73 
74         SCOTCH_Num *verttabGra;
75         SCOTCH_Num *edgetabGra;
76         verttabGra = new SCOTCH_Num[cellNbr + 1];
77         edgetabGra = new SCOTCH_Num[edgenbrGra];
78 
79         int pointer = 0;
80         int xyArea = dimensions[0] * dimensions[1];
81         for(int i = 0; i < cellNbr; ++i){
82             verttabGra[i] = pointer;
83             if(i % dimensions[0] != 0){
84                 edgetabGra[pointer] = i - 1;
85                 pointer++;
86             }
87             if(i % dimensions[0] != (dimensions[0] - 1)){
88                 edgetabGra[pointer] = i + 1;
89                 pointer++;
90             }
91             if(!((i % xyArea) < dimensions[0])){
92                 edgetabGra[pointer] = i - dimensions[0];
93                 pointer++;
94             }
95             if(!((i % xyArea) >= dimensions[0] * (dimensions[1] - 1))){
96                 edgetabGra[pointer] = i + dimensions[0];
97                 pointer++;
98             }
99             if(DIM == 3 && i >= (dimensions[0] * dimensions[1])){
100                 edgetabGra[pointer] = i - (dimensions[0] * dimensions[1]);
101                 pointer++;
102             }
103             if(DIM == 3 && i < (dimensions[0] *
104                                 dimensions[1]) * (dimensions[2] - 1)){
105                 edgetabGra[pointer] = i + (dimensions[0] * dimensions[1]);
106                 pointer++;
107             }
108         }
109         verttabGra[cellNbr] = pointer;
110 
111         SCOTCH_dgraphBuild(&grafdat,
112                            0,
113                            cellNbr,
114                            cellNbr,
115                            verttabGra,
116                            verttabGra +1,
117                            NULL,
118                            NULL,
119                            edgenbrGra,
120                            edgenbrGra,
121                            edgetabGra,
122                            NULL,
123                            NULL);
124 
125 
126         SCOTCH_Strat *straptr = SCOTCH_stratAlloc();;
127         SCOTCH_stratInit(straptr);
128 
129         SCOTCH_dgraphMap(&grafdat, &arch, straptr, indices);
130 
131         SCOTCH_archExit(&arch);
132         SCOTCH_dgraphExit(&grafdat);
133         SCOTCH_stratExit(straptr);
134         delete[] velotabArch;
135         delete[] verttabGra;
136         delete[] edgetabGra;
137 
138     }
139 
createRegions(SCOTCH_Num * indices)140     void createRegions(SCOTCH_Num *indices)
141     {
142         int rank = indices[0];
143         int length = 0;
144         int start = 0;
145         for(int i = 1; i < cellNbr + 1; ++i){
146             if((rank == indices[i]) &&
147                (i < cellNbr) && (i % dimensions[0] != 0)){
148                 length++;
149             } else {
150                 length++;
151                 Coord<DIM> startCoord;
152                 Coord<DIM> lengthCoord;
153                 lengthCoord[0] = length;
154                 lengthCoord[1] = 1;
155                 if(DIM == 3){
156                     startCoord[0] = origin[0] +
157                         (start % (dimensions[0] * dimensions[1])) %
158                         dimensions[0];
159                     startCoord[1] = origin[1] +
160                         (start % (dimensions[0] * dimensions[1])) /
161                         dimensions[0];
162                     startCoord[2] = origin[2] + start /
163                         (dimensions[0] *
164                          dimensions[1]);
165                     lengthCoord[2] = 1;
166                 } else {
167                     startCoord[0] = origin[0] + start % dimensions[0];
168                     startCoord[1] = origin[1] + start / dimensions[0];
169                 }
170 
171 
172                 regions[rank] <<
173                     CoordBox<DIM>(startCoord, lengthCoord);
174                 rank = indices[i];
175                 start = i;
176                 length = 0;
177             }
178         }
179     }
180  };
181 }
182 
183 #endif
184 
185 #endif
186 
187 #endif
188