COMBINATORIAL_BLAS 1.6
 
Loading...
Searching...
No Matches
BPMaximalMatching.h
Go to the documentation of this file.
1#ifndef BP_MAXIMAL_MATCHING_H
2#define BP_MAXIMAL_MATCHING_H
3
4#include "CombBLAS.h"
5#include <iostream>
6#include <functional>
7#include <algorithm>
8#include <vector>
9#include <limits>
10#include "Utility.h"
11#include "MatchingDefs.h"
12
13#define NO_INIT 0
14#define GREEDY 1
15#define KARP_SIPSER 2
16#define DMD 3
17
18namespace combblas {
19
20// This is not tested with CSC yet
21// TODO: test with CSC and Setting SPA (similar to Weighted Greedy)
22template <typename Par_DCSC_Bool, typename IT>
23void MaximalMatching(Par_DCSC_Bool & A, Par_DCSC_Bool & AT, FullyDistVec<IT, IT>& mateRow2Col,
24 FullyDistVec<IT, IT>& mateCol2Row, FullyDistVec<IT, IT>& degColRecv, int type, bool rand=true)
25{
26 static MTRand GlobalMT(123); // for reproducible result
27 typedef VertexTypeML < IT, IT> VertexType;
28 int nprocs, myrank;
29 MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
30 MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
31 int nthreads = 1;
32#ifdef _OPENMP
33#pragma omp parallel
34 {
35 nthreads = omp_get_num_threads();
36 }
37#endif
38
39 FullyDistVec<IT, IT> degCol = degColRecv;
40
41 //unmatched row and column vertices
42 FullyDistSpVec<IT, IT> unmatchedRow(mateRow2Col, [](IT mate){return mate==-1;});
43 FullyDistSpVec<IT, IT> degColSG(A.getcommgrid(), A.getncol());
44 //FullyDistVec<IT, IT> degCol(A.getcommgrid());
45 //A.Reduce(degCol, Column, plus<IT>(), static_cast<IT>(0)); // Reduce is not multithreaded
46
47
48 FullyDistSpVec<IT, VertexType> unmatchedCol(A.getcommgrid(), A.getncol());
49 // every veretx is unmatched. keep non-isolated vertices
50 unmatchedCol = EWiseApply<VertexType>(unmatchedCol, degCol,
51 [](VertexType vtx, IT deg){return VertexType();},
52 [](VertexType vtx, IT deg){return deg>0;},
53 true, VertexType());
54
55
56 FullyDistSpVec<IT, VertexType> fringeRow(A.getcommgrid(), A.getnrow());
57 FullyDistSpVec<IT, IT> fringeRow2(A.getcommgrid(), A.getnrow());
58 FullyDistSpVec<IT, VertexType> deg1Col(A.getcommgrid(), A.getncol());
59
60
61 IT curUnmatchedCol = unmatchedCol.getnnz();
62 IT curUnmatchedRow = unmatchedRow.getnnz();
63 IT newlyMatched = 1; // ensure the first pass of the while loop
64 int iteration = 0;
65 double tStart = MPI_Wtime();
66 std::vector<std::vector<double> > timing;
67
68#ifdef DETAIL_STATS
69 if(myrank == 0)
70 {
71 cout << "=======================================================\n";
72 cout << "@@@@@@ Number of processes: " << nprocs << endl;
73 cout << "=======================================================\n";
74 cout << "It | UMRow | UMCol | newlyMatched | Time "<< endl;
75 cout << "=======================================================\n";
76 }
77#endif
78 MPI_Barrier(MPI_COMM_WORLD);
79
80
81 while(curUnmatchedCol !=0 && curUnmatchedRow!=0 && newlyMatched != 0 )
82 {
83 unmatchedCol.ApplyInd([](VertexType vtx, IT idx){return VertexType(idx,idx);});
84 if(type==DMD)
85 {
86 unmatchedCol = EWiseApply<VertexType>(unmatchedCol, degCol,
87 [](VertexType vtx, IT deg){return VertexType(vtx.parent,deg);},
88 [](VertexType vtx, IT deg){return true;},
89 false, VertexType());
90 }
91 else if(rand)
92 {
93 unmatchedCol.Apply([](VertexType vtx){return VertexType(vtx.parent, static_cast<IT>((GlobalMT.rand() * 9999999)+1));});
94 }
95
96 // ======================== step1: One step of BFS =========================
97 std::vector<double> times;
98 double t1 = MPI_Wtime();
99 if(type==GREEDY)
100 {
101 SpMV<Select2ndMinSR<bool, VertexType>>(A, unmatchedCol, fringeRow, false);
102 }
103 else if(type==DMD)
104 {
105 SpMV<Select2ndMinSR<bool, VertexType>>(A, unmatchedCol, fringeRow, false);
106 }
107 else //(type==KARP_SIPSER)
108 {
109 deg1Col = EWiseApply<VertexType>(unmatchedCol, degCol,
110 [](VertexType vtx, IT deg){return vtx;},
111 [](VertexType vtx, IT deg){return deg==1;},
112 false, VertexType());
113
114 if(deg1Col.getnnz()>9)
115 SpMV<Select2ndMinSR<bool, VertexType>>(A, deg1Col, fringeRow, false);
116 else
117 SpMV<Select2ndMinSR<bool, VertexType>>(A, unmatchedCol, fringeRow, false);
118 }
119 // Remove matched row vertices
120 fringeRow = EWiseApply<VertexType>(fringeRow, mateRow2Col,
121 [](VertexType vtx, IT mate){return vtx;},
122 [](VertexType vtx, IT mate){return mate==-1;},
123 false, VertexType());
124
125 if(myrank == 0){times.push_back(MPI_Wtime()-t1); t1 = MPI_Wtime();}
126 // ===========================================================================
127
128
129 // ======================== step2: Update matching =========================
130
131 fringeRow2 = EWiseApply<IT>(fringeRow, mateRow2Col,
132 [](VertexType vtx, IT mate){return vtx.parent;},
133 [](VertexType vtx, IT mate){return true;},
134 false, VertexType());
135
136 FullyDistSpVec<IT, IT> newMatchedCols = fringeRow2.Invert(A.getncol());
137 FullyDistSpVec<IT, IT> newMatchedRows = newMatchedCols.Invert(A.getnrow());
138 mateCol2Row.Set(newMatchedCols);
139 mateRow2Col.Set(newMatchedRows);
140 if(myrank == 0){times.push_back(MPI_Wtime()-t1); t1 = MPI_Wtime();}
141 // ===========================================================================
142
143
144 // =============== step3: Update degree of unmatched columns =================
145 unmatchedRow.Select(mateRow2Col, [](IT mate){return mate==-1;});
146 unmatchedCol.Select(mateCol2Row, [](IT mate){return mate==-1;});
147
148 if(type!=GREEDY)
149 {
150 // update degree
151 newMatchedRows.Apply([](IT val){return 1;}); // needed if the matrix is Boolean since the SR::multiply isn't called
152 SpMV< SelectPlusSR<bool, IT>>(AT, newMatchedRows, degColSG, false); // degree of column vertices to matched rows
153 // subtract degree of column vertices
154 degCol.EWiseApply(degColSG,
155 [](IT old_deg, IT new_deg){return old_deg-new_deg;},
156 [](IT old_deg, IT new_deg){return true;},
157 false, static_cast<IT>(0));
158 // remove isolated vertices
159 unmatchedCol = EWiseApply<VertexType>(unmatchedCol, degCol,
160 [](VertexType vtx, IT deg){return vtx;},
161 [](VertexType vtx, IT deg){return deg>0;},
162 false, VertexType());
163 }
164 if(myrank == 0){times.push_back(MPI_Wtime()-t1); t1 = MPI_Wtime();}
165 // ===========================================================================
166
167
168 ++iteration;
169 newlyMatched = newMatchedCols.getnnz();
170 times.push_back(std::accumulate(times.begin(), times.end(), 0.0));
171 timing.push_back(times);
172#ifdef DETAIL_STATS
173 if(myrank == 0)
174 {
175 printf("%3d %10lld %10lld %10lld %18lf\n", iteration , curUnmatchedRow, curUnmatchedCol, newlyMatched, times.back());
176 }
177#endif
178 curUnmatchedCol = unmatchedCol.getnnz();
179 curUnmatchedRow = unmatchedRow.getnnz();
180 MPI_Barrier(MPI_COMM_WORLD);
181
182 }
183
184 IT cardinality = mateRow2Col.Count([](IT mate){return mate!=-1;});
185 std::vector<double> totalTimes(timing[0].size(),0);
186 for(int i=0; i<timing.size(); i++)
187 {
188 for(int j=0; j<timing[i].size(); j++)
189 {
190 totalTimes[j] += timing[i][j];
191 }
192 }
193
194
195 if(myrank == 0)
196 {
197#ifdef DETAIL_STATS
198 cout << "==========================================================\n";
199 cout << "\n================individual timings =======================\n";
200 cout << " SpMV Update-Match Update-UMC Total "<< endl;
201 cout << "==========================================================\n";
202 for(int i=0; i<timing.size(); i++)
203 {
204 for(int j=0; j<timing[i].size(); j++)
205 {
206 printf("%12.5lf ", timing[i][j]);
207 }
208 cout << endl;
209 }
210
211 cout << "-------------------------------------------------------\n";
212 for(int i=0; i<totalTimes.size(); i++)
213 printf("%12.5lf ", totalTimes[i]);
214 cout << endl;
215#endif
216 std::cout << "****** maximal matching runtime ********\n";
217 std::cout << "nprocesses nthreads ncores algorithm Unmatched-Rows Cardinality Total Time***\n";
218 std::cout << nprocs << " " << nthreads << " " << nprocs * nthreads << " ";
219 if(type == DMD) std::cout << "DMD";
220 else if(type == GREEDY) std::cout << "Greedy";
221 else if(type == KARP_SIPSER) std::cout << "Karp-Sipser";
222 if(rand && (type == KARP_SIPSER || type == GREEDY) ) std::cout << "-rand";
223 std::cout << " ";
224 printf("%lld %lld %lf\n", curUnmatchedRow, cardinality, totalTimes.back());
225 std::cout << "-------------------------------------------------------\n\n";
226 }
227 //isMatching(mateCol2Row, mateRow2Col);
228}
229
230
231
232// Special version of the greedy algorithm (works for both CSC (multithreaded) and DCSC)
233// That uses WeightMaxSR semiring
234// TODO: check if this is 1/2 approx of the weighted matching (probably no)
235// TODO: should we remove degCol?
236// TODO: can be merged with the generalized MaximalMatching
237template <typename Par_MAT_Double, typename IT>
238void WeightedGreedy(Par_MAT_Double & A, FullyDistVec<IT, IT>& mateRow2Col,
239 FullyDistVec<IT, IT>& mateCol2Row, FullyDistVec<IT, IT>& degCol)
240{
241
242 typedef VertexTypeML < IT, double> VertexType;
243 int nthreads=1;
244#ifdef THREADED
245#pragma omp parallel
246 {
247 nthreads = omp_get_num_threads();
248 }
249#endif
250 PreAllocatedSPA<VertexType> SPA(A.seq(), nthreads*4);
251 int nprocs, myrank;
252 MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
253 MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
254
255 double tStart = MPI_Wtime();
256
257 //unmatched row and column vertices
258 FullyDistSpVec<IT, IT> unmatchedRow(mateRow2Col, [](IT mate){return mate==-1;});
259
260
261
262 FullyDistSpVec<IT, VertexType> unmatchedCol(A.getcommgrid(), A.getncol());
263 // every veretx is unmatched. keep non-isolated vertices
264 unmatchedCol = EWiseApply<VertexType>(unmatchedCol, degCol,
265 [](VertexType vtx, IT deg){return VertexType();},
266 [](VertexType vtx, IT deg){return deg>0;},
267 true, VertexType());
268
269
270 FullyDistSpVec<IT, VertexType> fringeRow(A.getcommgrid(), A.getnrow());
271 FullyDistSpVec<IT, IT> fringeRow2(A.getcommgrid(), A.getnrow());
272 FullyDistSpVec<IT, VertexType> fringeRow3(A.getcommgrid(), A.getnrow());
273
274 IT curUnmatchedCol = unmatchedCol.getnnz();
275 IT curUnmatchedRow = unmatchedRow.getnnz();
276 IT newlyMatched = 1; // ensure the first pass of the while loop
277 int iteration = 0;
278
279 std::vector<std::vector<double> > timing;
280
281#ifdef DETAIL_STATS
282 if(myrank == 0)
283 {
284 cout << "=======================================================\n";
285 cout << "@@@@@@ Number of processes: " << nprocs << endl;
286 cout << "=======================================================\n";
287 cout << "It | UMRow | UMCol | newlyMatched | Time "<< endl;
288 cout << "=======================================================\n";
289 }
290#endif
291 MPI_Barrier(MPI_COMM_WORLD);
292
293
294 while(curUnmatchedCol !=0 && curUnmatchedRow!=0 && newlyMatched != 0 )
295 {
296 // anything is fine in the second argument
297 unmatchedCol.ApplyInd([](VertexType vtx, IT idx){return VertexType(idx,idx);});
298
299
300 // ======================== step1: One step of BFS =========================
301 std::vector<double> times;
302 double t1 = MPI_Wtime();
303
304 SpMV<WeightMaxMLSR<double, VertexType>>(A, unmatchedCol, fringeRow, false, SPA);
305
306 // Remove matched row vertices
307 fringeRow = EWiseApply<VertexType>(fringeRow, mateRow2Col,
308 [](VertexType vtx, IT mate){return vtx;},
309 [](VertexType vtx, IT mate){return mate==-1;},
310 false, VertexType());
311
312 if(myrank == 0){times.push_back(MPI_Wtime()-t1); t1 = MPI_Wtime();}
313 // ===========================================================================
314
315
316 // ======================== step2: Update matching =========================
317
318 fringeRow2 = EWiseApply<IT>(fringeRow, mateRow2Col,
319 [](VertexType vtx, IT mate){return vtx.parent;},
320 [](VertexType vtx, IT mate){return true;},
321 false, VertexType());
322
323 FullyDistSpVec<IT, IT> newMatchedCols = fringeRow2.Invert(A.getncol());
324 FullyDistSpVec<IT, IT> newMatchedRows = newMatchedCols.Invert(A.getnrow());
325 mateCol2Row.Set(newMatchedCols);
326 mateRow2Col.Set(newMatchedRows);
327 if(myrank == 0){times.push_back(MPI_Wtime()-t1); t1 = MPI_Wtime();}
328 // ===========================================================================
329
330
331 // =============== step3: Update unmatched columns and rows =================
332
333 unmatchedRow.Select(mateRow2Col, [](IT mate){return mate==-1;});
334 unmatchedCol.Select(mateCol2Row, [](IT mate){return mate==-1;});
335 if(myrank == 0){times.push_back(MPI_Wtime()-t1); t1 = MPI_Wtime();}
336 // ===========================================================================
337
338
339 ++iteration;
340 newlyMatched = newMatchedCols.getnnz();
341 times.push_back(std::accumulate(times.begin(), times.end(), 0.0));
342 timing.push_back(times);
343#ifdef DETAIL_STATS
344 if(myrank == 0)
345 {
346 printf("%3d %10lld %10lld %10lld %18lf\n", iteration , curUnmatchedRow, curUnmatchedCol, newlyMatched, times.back());
347 }
348#endif
349 curUnmatchedCol = unmatchedCol.getnnz();
350 curUnmatchedRow = unmatchedRow.getnnz();
351 MPI_Barrier(MPI_COMM_WORLD);
352
353 }
354
355 IT cardinality = mateRow2Col.Count([](IT mate){return mate!=-1;});
356 std::vector<double> totalTimes(timing[0].size(),0);
357 for(int i=0; i<timing.size(); i++)
358 {
359 for(int j=0; j<timing[i].size(); j++)
360 {
361 totalTimes[j] += timing[i][j];
362 }
363 }
364
365
366 if(myrank == 0)
367 {
368#ifdef DETAIL_STATS
369 cout << "==========================================================\n";
370 cout << "\n================individual timings =======================\n";
371 cout << " SpMV Update-Match Update-UMC Total "<< endl;
372 cout << "==========================================================\n";
373 for(int i=0; i<timing.size(); i++)
374 {
375 for(int j=0; j<timing[i].size(); j++)
376 {
377 printf("%12.5lf ", timing[i][j]);
378 }
379 cout << endl;
380 }
381
382 cout << "-------------------------------------------------------\n";
383 for(int i=0; i<totalTimes.size(); i++)
384 printf("%12.5lf ", totalTimes[i]);
385 cout << endl;
386#endif
387 std::cout << "****** maximal matching runtime ********\n";
388 std::cout << "Unmatched-Rows Cardinality Total Time***\n";
389 printf("%lld %lld %lf\n", curUnmatchedRow, cardinality, totalTimes.back());
390 std::cout << "-------------------------------------------------------\n\n";
391 }
392 //isMatching(mateCol2Row, mateRow2Col);
393}
394
395
396
397
398template <class Par_DCSC_Bool, class IT, class NT>
399bool isMaximalmatching(Par_DCSC_Bool & A, FullyDistVec<IT,NT> & mateRow2Col, FullyDistVec<IT,NT> & mateCol2Row)
400{
401 typedef VertexTypeML < IT, IT> VertexType;
402 int myrank;
403 MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
404 FullyDistSpVec<IT, IT> fringeRow(A.getcommgrid(), A.getnrow());
405 FullyDistSpVec<IT, IT> fringeCol(A.getcommgrid(), A.getncol());
406 FullyDistSpVec<IT, IT> unmatchedRow(mateRow2Col, [](IT mate){return mate==-1;});
407 FullyDistSpVec<IT, IT> unmatchedCol(mateCol2Row, [](IT mate){return mate==-1;});
408 unmatchedRow.setNumToInd();
409 unmatchedCol.setNumToInd();
410
411
412 SpMV<Select2ndMinSR<bool, VertexType>>(A, unmatchedCol, fringeRow, false);
413 fringeRow = EWiseMult(fringeRow, mateRow2Col, true, (IT) -1);
414 if(fringeRow.getnnz() != 0)
415 {
416 if(myrank == 0)
417 std::cout << "Not maximal matching!!\n";
418 return false;
419 }
420
421 Par_DCSC_Bool tA = A;
422 tA.Transpose();
423 SpMV<Select2ndMinSR<bool, VertexType>>(tA, unmatchedRow, fringeCol, false);
424 fringeCol = EWiseMult(fringeCol, mateCol2Row, true, (IT) -1);
425 if(fringeCol.getnnz() != 0)
426 {
427 if(myrank == 0)
428 std::cout << "Not maximal matching**!!\n";
429 return false;
430 }
431 return true;
432}
433
434}
435
436#endif
437
GREEDY
#define GREEDY
Definition BPMaximalMatching.h:14
KARP_SIPSER
#define KARP_SIPSER
Definition BPMaximalMatching.h:15
DMD
#define DMD
Definition BPMaximalMatching.h:16
GlobalMT
MTRand GlobalMT
Definition RestrictionOp.h:18
MTRand::rand
double rand()
Definition MersenneTwister.h:146
size
int size
Definition common.h:20
SpMV
@ SpMV
Definition papi_combblas_globals.h:14
combblas::MaximalMatching
void MaximalMatching(Par_DCSC_Bool &A, Par_DCSC_Bool &AT, FullyDistVec< IT, IT > &mateRow2Col, FullyDistVec< IT, IT > &mateCol2Row, FullyDistVec< IT, IT > &degColRecv, int type, bool rand=true)
Definition BPMaximalMatching.h:23
combblas::isMaximalmatching
bool isMaximalmatching(Par_DCSC_Bool &A, FullyDistVec< IT, NT > &mateRow2Col, FullyDistVec< IT, NT > &mateCol2Row)
Definition BPMaximalMatching.h:399
combblas::WeightedGreedy
void WeightedGreedy(Par_MAT_Double &A, FullyDistVec< IT, IT > &mateRow2Col, FullyDistVec< IT, IT > &mateCol2Row, FullyDistVec< IT, IT > &degCol)
Definition BPMaximalMatching.h:238
combblas::EWiseMult
Dcsc< IU, typename promote_trait< NU1, NU2 >::T_promote > EWiseMult(const Dcsc< IU, NU1 > &A, const Dcsc< IU, NU2 > *B, bool exclude)
Definition Friends.h:694
A
double A