1 #ifndef RESTRICTION_OP_H
2 #define RESTRICTION_OP_H
3
4 //#define DETERMINISTIC
5 #include "CombBLAS/CombBLAS.h"
6 #ifdef THREADED
7 #ifndef _OPENMP
8 #define _OPENMP
9 #endif
10 #include <omp.h>
11 #endif
12
13
14
15 #ifdef DETERMINISTIC
16 MTRand GlobalMT(1);
17 #else
18 MTRand GlobalMT; // generate random numbers with Mersenne Twister
19 #endif
20
21
22 namespace combblas {
23
24 template <typename T1, typename T2>
25 struct Select2ndMinSR
26 {
idSelect2ndMinSR27 static T2 id(){ return T2(); };
returnedSAIDSelect2ndMinSR28 static bool returnedSAID() { return false; }
mpi_opSelect2ndMinSR29 static MPI_Op mpi_op() { return MPI_MIN; };
30
addSelect2ndMinSR31 static T2 add(const T2 & arg1, const T2 & arg2)
32 {
33 return std::min(arg1, arg2);
34 }
35
multiplySelect2ndMinSR36 static T2 multiply(const T1 & arg1, const T2 & arg2)
37 {
38 return arg2;
39 }
40
axpySelect2ndMinSR41 static void axpy(const T1 a, const T2 & x, T2 & y)
42 {
43 y = add(y, multiply(a, x));
44 }
45 };
46
47 template <typename T>
48 struct VertexType
49 {
50 public:
VertexTypeVertexType51 VertexType(){parent=-1; prob=0.0;};
VertexTypeVertexType52 VertexType(T pa, double pr){parent=pa; prob=pr;};
53 friend bool operator==(const VertexType & vtx1, const VertexType & vtx2 ){return vtx1.parent==vtx2.parent;};
54 friend bool operator<(const VertexType & vtx1, const VertexType & vtx2 ){return vtx1.parent<vtx2.parent;};
55 friend std::ostream& operator<<(std::ostream& os, const VertexType & vertex ){os << "(" << vertex.parent << "," << vertex.prob << ")"; return os;};
VertexTypeVertexType56 VertexType(T pa){parent=pa; prob=0.0;};
57 T parent;
58 double prob;
59 };
60
61
62
63 template <typename T1, typename T2>
64 struct Select2ndRandSR
65 {
idSelect2ndRandSR66 static VertexType<T2> id(){ return VertexType<T2>(); };
returnedSAIDSelect2ndRandSR67 static bool returnedSAID() { return false; }
68 //static MPI_Op mpi_op() { return MPI_MIN; };
69
addSelect2ndRandSR70 static VertexType<T2> add(const VertexType<T2> & arg1, const VertexType<T2> & arg2)
71 {
72 if((arg1.prob) < (arg2.prob)) return arg1;
73 else return arg2;
74 }
multiplySelect2ndRandSR75 static VertexType<T2> multiply(const T1 & arg1, const VertexType<T2> & arg2)
76 {
77 return arg2;
78 }
79
axpySelect2ndRandSR80 static void axpy(T1 a, const VertexType<T2> & x, VertexType<T2> & y)
81 {
82 y = add(y, multiply(a, x));
83 }
84 };
85
86
87
88
89 template <typename T1, typename T2>
90 struct MIS2verifySR // identical to Select2ndMinSR except for the printout in add()
91 {
idMIS2verifySR92 static T2 id(){ return T2(); };
returnedSAIDMIS2verifySR93 static bool returnedSAID() { return false; }
mpi_opMIS2verifySR94 static MPI_Op mpi_op() { return MPI_MIN; };
95
addMIS2verifySR96 static T2 add(const T2 & arg1, const T2 & arg2)
97 {
98 std::cout << "This should have never been executed for MIS-2 to be correct" << std::endl;
99 return std::min(arg1, arg2);
100 }
101
multiplyMIS2verifySR102 static T2 multiply(const T1 & arg1, const T2 & arg2)
103 {
104 return arg2;
105 }
106
axpyMIS2verifySR107 static void axpy(const T1 a, const T2 & x, T2 & y)
108 {
109 y = add(y, multiply(a, x));
110 }
111 };
112
113
114
115
116 // second hop MIS (i.e., MIS on A \Union A^2)
117 template <typename ONT, typename IT, typename INT, typename DER>
MIS2(SpParMat<IT,INT,DER> A)118 FullyDistSpVec<IT, ONT> MIS2(SpParMat < IT, INT, DER> A)
119 {
120 IT nvert = A.getncol();
121
122 //# the final result set. S[i] exists and is 1 if vertex i is in the MIS
123 FullyDistSpVec<IT, ONT> mis ( A.getcommgrid(), nvert);
124
125
126 //# the candidate set. initially all vertices are candidates.
127 //# If cand[i] exists, then i is a candidate. The value cand[i] is i's random number for this iteration.
128 FullyDistSpVec<IT, double> cand(A.getcommgrid());
129 cand.iota(nvert, 1.0); // any value is fine since we randomize it later
130 FullyDistSpVec<IT, double> min_neighbor_r ( A.getcommgrid(), nvert);
131 FullyDistSpVec<IT, double> min_neighbor2_r ( A.getcommgrid(), nvert);
132
133
134 FullyDistSpVec<IT, ONT> new_S_members ( A.getcommgrid(), nvert);
135 FullyDistSpVec<IT, ONT> new_S_neighbors ( A.getcommgrid(), nvert);
136 FullyDistSpVec<IT, ONT> new_S_neighbors2 ( A.getcommgrid(), nvert);
137
138
139 while (cand.getnnz() > 0)
140 {
141
142 //# label each vertex in cand with a random value (in what range, [0,1])
143 cand.Apply([](const double & ignore){return (double) GlobalMT.rand();});
144
145 //# find the smallest random value among a vertex's 1 and 2-hop neighbors
146 SpMV<Select2ndMinSR<INT, double>>(A, cand, min_neighbor_r, false);
147 SpMV<Select2ndMinSR<INT, double>>(A, min_neighbor_r, min_neighbor2_r, false);
148
149 FullyDistSpVec<IT, double> min_neighbor_r_union = EWiseApply<double>(min_neighbor2_r, min_neighbor_r,
150 [](double x, double y){return std::min(x,y);},
151 [](double x, double y){return true;}, // do_op is totalogy
152 true, true, 2.0, 2.0, true); // we allow nulls for both V and W
153
154
155 //# The vertices to be added to S this iteration are those whose random value is
156 //# smaller than those of all its 1-hop and 2-hop neighbors:
157 // **** if cand has isolated vertices, they will be included in new_S_members ******
158 new_S_members = EWiseApply<ONT>(min_neighbor_r_union, cand,
159 [](double x, double y){return (ONT)1;},
160 [](double x, double y){return y<=x;}, // equality is for back edges since we are operating on A^2
161 true, false, 2.0, 2.0, true);
162
163 //# new_S_members are no longer candidates, so remove them from cand
164 cand = EWiseApply<double>(cand, new_S_members,
165 [](double x, ONT y){return x;},
166 [](double x, ONT y){return true;},
167 false, true, 0.0, (ONT) 0, false);
168
169 //# find 1-hop and 2-hop neighbors of new_S_members
170 SpMV<Select2ndMinSR<INT, ONT>>(A, new_S_members, new_S_neighbors, false);
171 SpMV<Select2ndMinSR<INT, ONT>>(A, new_S_neighbors, new_S_neighbors2, false);
172
173 FullyDistSpVec<IT, ONT> new_S_neighbors_union = EWiseApply<ONT>(new_S_neighbors, new_S_neighbors2,
174 [](ONT x, ONT y){return x;}, // in case of intersection, doesn't matter which one to propagate
175 [](ONT x, ONT y){return true;},
176 true, true, (ONT) 1, (ONT) 1, true);
177
178
179 //# remove 1-hop and 2-hop neighbors of new_S_members from cand, because they cannot be part of the MIS anymore
180 cand = EWiseApply<double>(cand, new_S_neighbors_union,
181 [](double x, ONT y){return x;},
182 [](double x, ONT y){return true;},
183 false, true, 0.0, (ONT) 0, false);
184
185 //# add new_S_members to mis
186 mis = EWiseApply<ONT>(mis, new_S_members,
187 [](ONT x, ONT y){return x;},
188 [](ONT x, ONT y){return true;},
189 true, true, (ONT) 1, (ONT) 1, true);
190 }
191
192 return mis;
193 }
194
195
196 template <typename IT, typename NT>
RestrictionOp(CCGrid & CMG,SpDCCols<IT,NT> * localmat,SpDCCols<IT,NT> * & R,SpDCCols<IT,NT> * & RT)197 void RestrictionOp( CCGrid & CMG, SpDCCols<IT, NT> * localmat, SpDCCols<IT, NT> *& R, SpDCCols<IT, NT> *& RT)
198 {
199 if(CMG.layer_grid == 0)
200 {
201 SpDCCols<IT, bool> *A = new SpDCCols<IT, bool>(*localmat);
202
203 SpParMat < IT, bool, SpDCCols < IT, bool >> B (A, CMG.layerWorld);
204
205 B.RemoveLoops();
206
207 SpParMat < IT, bool, SpDCCols < IT, bool >> BT = B;
208 BT.Transpose();
209 B += BT;
210 B.PrintInfo();
211
212
213 FullyDistSpVec<IT,IT>mis2 = MIS2<IT>(B);
214 mis2.setNumToInd();
215 mis2.PrintInfo("mis2");
216 FullyDistSpVec<IT,IT> mis2neigh = SpMV<MIS2verifySR<bool, IT>>(B, mis2, false);
217
218 // union of mis2 and mis2neigh
219 mis2neigh = EWiseApply<IT>(mis2neigh, mis2,
220 [](IT x, IT y){return x==-1?y:x;},
221 [](IT x, IT y){return true;},
222 true, true, (IT) -1, (IT) -1, true);
223
224 // mis2neigh with a probability
225 FullyDistSpVec<IT, VertexType<IT>> mis2neigh_p(mis2neigh.getcommgrid(), mis2neigh.TotalLength());
226 mis2neigh_p = EWiseApply<VertexType<IT>>(mis2neigh, mis2neigh_p,
227 [](IT x, VertexType<IT> y){return VertexType<IT>(x,GlobalMT.rand());},
228 [](IT x, VertexType<IT> y){return true;},
229 false, true, (IT) -1, VertexType<IT>(), false);
230
231 // mis2neigh2 with a probability
232 FullyDistSpVec<IT, VertexType<IT>> mis2neigh2_p(mis2neigh.getcommgrid(), mis2neigh.TotalLength());
233 SpMV<Select2ndRandSR<bool, IT>>(B, mis2neigh_p, mis2neigh2_p, false);
234
235 // mis2neigh2 without probability
236 FullyDistSpVec<IT,IT> mis2neigh2(mis2neigh.getcommgrid(), mis2neigh.TotalLength());
237 mis2neigh2 = EWiseApply<IT>(mis2neigh2, mis2neigh2_p,
238 [](IT x, VertexType<IT> y){return y.parent;},
239 [](IT x, VertexType<IT> y){return true;},
240 true, false, (IT) -1, VertexType<IT>(), false);
241
242
243 // union of mis2 and mis2neigh and mis2neigh2
244 FullyDistSpVec<IT,IT> mis2neighUnion = EWiseApply<IT>(mis2neigh, mis2neigh2,
245 [](IT x, IT y){return x==-1?y:x;},
246 [](IT x, IT y){return true;},
247 true, true, (IT) -1, (IT) -1, true);
248
249 mis2neighUnion.PrintInfo("mis2neighUnion");
250 if(mis2neighUnion.getnnz() != mis2neighUnion.TotalLength())
251 {
252 SpParHelper::Print(" !!!! Error: mis2neighUnion does not include all rows/columns. !!!! ");
253 }
254
255 // At first, create nxn matrix
256 FullyDistVec<IT, IT> ci = mis2neighUnion;
257 FullyDistVec<IT,IT> ri = mis2neighUnion.FindInds([](IT x){return true;}); // this should be equivalent to iota
258 SpParMat<IT,NT,SpDCCols<IT,NT>> Rop(B.getnrow(), ci.TotalLength(), ri, ci, (NT)1, false);
259
260 // next, select nonempty columns
261 FullyDistVec<IT, IT> cimis2 = mis2.FindInds([](IT x){return true;}); // nonzero columns
262 Rop(ri,cimis2,true);
263 SpParHelper::Print("Rop final (before normalization)... ");
264 Rop.PrintInfo();
265
266 // permute for load balance
267 float balance_before = Rop.LoadImbalance();
268 FullyDistVec<IT, IT> perm_row(Rop.getcommgrid()); // permutation vector defined on layers
269 FullyDistVec<IT, IT> perm_col(Rop.getcommgrid()); // permutation vector defined on layers
270
271 perm_row.iota(Rop.getnrow(), 0); // don't permute rows because they represent the IDs of "fine" vertices
272 perm_col.iota(Rop.getncol(), 0); // CAN permute columns because they define the IDs of new aggregates
273 perm_col.RandPerm(); // permuting columns for load balance
274
275 Rop(perm_row, perm_col, true); // in place permute
276 float balance_after = Rop.LoadImbalance();
277
278 std::ostringstream outs;
279 outs << "Load balance (before): " << balance_before << std::endl;
280 outs << "Load balance (after): " << balance_after << std::endl;
281 SpParHelper::Print(outs.str());
282
283
284 SpParMat<IT,NT,SpDCCols<IT,NT>> RopT = Rop;
285 RopT.Transpose();
286
287 R = new SpDCCols<IT,NT>(Rop.seq()); // deep copy
288 RT = new SpDCCols<IT,NT>(RopT.seq()); // deep copy
289
290 }
291 }
292
293 // with added column
294 /*
295 template <typename IT, typename NT>
296 void RestrictionOp( CCGrid & CMG, SpDCCols<IT, NT> * localmat, SpDCCols<IT, NT> *& R, SpDCCols<IT, NT> *& RT)
297 {
298 if(CMG.layer_grid == 0)
299 {
300 SpDCCols<IT, bool> *A = new SpDCCols<IT, bool>(*localmat);
301
302 SpParMat < IT, bool, SpDCCols < IT, bool >> B (A, CMG.layerWorld);
303
304 B.RemoveLoops();
305
306 SpParMat < IT, bool, SpDCCols < IT, bool >> BT = B;
307 BT.Transpose();
308 B += BT;
309
310 // ------------ compute MIS-2 ----------------------------
311 FullyDistSpVec<IT, IT> mis2 (B.getcommgrid(), B.getncol()); // values of the mis2 vector are just "ones"
312 mis2 = MIS2<IT>(B);
313 mis2.PrintInfo("MIS original");
314 // ------------ Obtain restriction matrix from mis2 ----
315 FullyDistVec<IT,IT> ri = mis2.FindInds([](IT x){return true;});
316
317 // find the vertices that are not covered by mis2 AND its one hop neighborhood
318 FullyDistSpVec<IT,IT> mis2neigh = SpMV<MIS2verifySR<bool, IT>>(B, mis2, false);
319 mis2neigh.PrintInfo("MIS neighbors");
320
321 // ABAB: mis2 and mis2neigh should be independent, because B doesn't have any loops.
322 FullyDistSpVec<IT,IT> isection = EWiseApply<IT>(mis2neigh, mis2,
323 [](IT x, IT y){return x;},
324 [](IT x, IT y){return true;},
325 false, false, (IT) 1, (IT) 1, true); /// allowVNulls and allowWNulls are both false
326 isection.PrintInfo("intersection of mis2neigh and mis2");
327
328
329 // find the union of mis2neigh and mis2 (ABAB: this function to be wrapped & called "SetUnion")
330 mis2neigh = EWiseApply<IT>(mis2neigh, mis2,
331 [](IT x, IT y){return x;},
332 [](IT x, IT y){return true;},
333 true, true, (IT) 1, (IT) 1, true);
334 mis2neigh.PrintInfo("MIS original+neighbors");
335
336
337 // FullyDistVec<IT, NT>::FullyDistVec ( shared_ptr<CommGrid> grid, IT globallen, NT initval)
338 // : FullyDist<IT,NT,typename CombBLAS::disable_if< CombBLAS::is_boolean<NT>::value, NT >::type>(grid,globallen)
339
340 FullyDistVec<IT, IT> denseones(mis2neigh.getcommgrid(), B.getncol(), 1); // calls the default constructor... why?
341 FullyDistSpVec<IT,IT> spones (denseones);
342
343 // subtract the entries of mis2neigh from all vertices (ABAB: this function to be wrapped & called "SetDiff")
344 spones = EWiseApply<IT>(spones, mis2neigh,
345 [](IT x, IT y){return x;}, // binop
346 [](IT x, IT y){return true;}, // doop
347 false, true, (IT) 1, (IT) 1, false); // allowintersect=false (all joint entries are removed)
348
349 spones.PrintInfo("Leftovers (singletons)");
350
351
352 FullyDistVec<IT, IT> ci(B.getcommgrid());
353 ci.iota(mis2.getnnz(), (IT)0);
354 SpParMat<IT,NT,SpDCCols<IT,NT>> M(B.getnrow(), ci.TotalLength(), ri, ci, (NT)1, false);
355
356 SpParHelper::Print("M matrix... ");
357 M.PrintInfo();
358
359 SpParMat<IT,NT,SpDCCols<IT,NT>> Rop = PSpGEMM<PlusTimesSRing<bool, NT>>(B,M);
360
361 SpParHelper::Print("R (minus M) matrix... ");
362 Rop.PrintInfo();
363
364 Rop += M;
365
366 SpParHelper::Print("R without singletons... ");
367 Rop.PrintInfo();
368
369 FullyDistVec<IT,IT> rrow(Rop.getcommgrid());
370 FullyDistVec<IT,IT> rcol(Rop.getcommgrid());
371 FullyDistVec<IT,NT> rval(Rop.getcommgrid());
372 Rop.Find(rrow, rcol, rval);
373
374 FullyDistVec<IT, IT> extracols(Rop.getcommgrid());
375 extracols.iota(spones.getnnz(), ci.TotalLength()); // one column per singleton
376
377 // Returns a dense vector of nonzero global indices for which the predicate is satisfied on values
378 FullyDistVec<IT,IT> extrarows = spones.FindInds([](IT x){return true;}); // dense leftovers array is the extra rows
379
380 // Resize Rop
381 SpParMat<IT,NT,SpDCCols<IT,NT>> RopFull1(Rop.getnrow(), Rop.getncol() +extracols.TotalLength(), rrow, rcol, rval, false);
382
383 SpParHelper::Print("RopFull1... ");
384 RopFull1.PrintInfo();
385
386 SpParMat<IT,NT,SpDCCols<IT,NT>> RopFull2(Rop.getnrow(), Rop.getncol() +extracols.TotalLength(), extrarows, extracols, (NT)1, false);
387
388 SpParHelper::Print("RopFull2... ");
389 RopFull2.PrintInfo();
390
391 RopFull1 += RopFull2;
392
393 SpParHelper::Print("RopFull final (before normalization)... ");
394 RopFull1.PrintInfo();
395
396 float balance_before = RopFull1.LoadImbalance();
397
398 FullyDistVec<IT, IT> perm_row(RopFull1.getcommgrid()); // permutation vector defined on layers
399 FullyDistVec<IT, IT> perm_col(RopFull1.getcommgrid()); // permutation vector defined on layers
400
401 perm_row.iota(RopFull1.getnrow(), 0); // don't permute rows because they represent the IDs of "fine" vertices
402 perm_col.iota(RopFull1.getncol(), 0); // CAN permute columns because they define the IDs of new aggregates
403 perm_col.RandPerm(); // permuting columns for load balance
404
405 RopFull1(perm_row, perm_col, true); // in place permute
406 float balance_after = RopFull1.LoadImbalance();
407
408 ostringstream outs;
409 outs << "Load balance (before): " << balance_before << endl;
410 outs << "Load balance (after): " << balance_after << endl;
411 SpParHelper::Print(outs.str());
412
413
414 SpParMat<IT,NT,SpDCCols<IT,NT>> RopT = RopFull1;
415 RopT.Transpose();
416
417 R = new SpDCCols<IT,NT>(RopFull1.seq()); // deep copy
418 RT = new SpDCCols<IT,NT>(RopT.seq()); // deep copy
419
420 }
421 }
422 */
423
424 }
425
426 #endif
427
428