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