COMBINATORIAL_BLAS 1.6
 
Loading...
Searching...
No Matches
BPMaximalMatching.cpp
Go to the documentation of this file.
1#include "../CombBLAS.h"
2#include <mpi.h>
3#include <sys/time.h>
4#include <iostream>
5#include <functional>
6#include <algorithm>
7#include <vector>
8#include <string>
9#include <sstream>
10#include "BPMaximalMatching.h"
11
12
13#ifdef THREADED
14 #ifndef _OPENMP
15 #define _OPENMP
16 #endif
17
18 #include <omp.h>
19 int cblas_splits = 1;
20#endif
21
22using namespace std;
23using namespace combblas;
24
25
26bool prune, mvInvertMate, randMM, moreSplit;
27int init;
28bool randMaximal;
29bool fewexp;
30
31template <typename PARMAT>
32void Symmetricize(PARMAT & A)
33{
34 // boolean addition is practically a "logical or"
35 // therefore this doesn't destruct any links
36 PARMAT AT = A;
37 AT.Transpose();
38 A += AT;
39}
40
41
42struct VertexType
43{
44public:
45 VertexType(int64_t p=-1, int64_t r=-1, int16_t pr=0){parent=p; root = r; prob = pr;};
46
47 friend bool operator<(const VertexType & vtx1, const VertexType & vtx2 )
48 {
49 if(vtx1.prob==vtx2.prob) return vtx1.parent<vtx2.parent;
50 else return vtx1.prob<vtx2.prob;
51 };
52 friend bool operator==(const VertexType & vtx1, const VertexType & vtx2 ){return vtx1.parent==vtx2.parent;};
53 friend ostream& operator<<(ostream& os, const VertexType & vertex ){os << "(" << vertex.parent << "," << vertex.root << ")"; return os;};
54 //private:
55 int64_t parent;
56 int64_t root;
57 int16_t prob; // probability of selecting an edge
58
59};
60
61
62
63
64typedef SpParMat < int64_t, bool, SpDCCols<int64_t,bool> > PSpMat_Bool;
65
66/*
67 Remove isolated vertices and purmute
68 */
69void removeIsolated(PSpMat_Bool & A)
70{
71
72 int nprocs, myrank;
73 MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
74 MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
75
76
77 FullyDistVec<int64_t, int64_t> * ColSums = new FullyDistVec<int64_t, int64_t>(A.getcommgrid());
78 FullyDistVec<int64_t, int64_t> * RowSums = new FullyDistVec<int64_t, int64_t>(A.getcommgrid());
79 FullyDistVec<int64_t, int64_t> nonisoRowV; // id's of non-isolated (connected) Row vertices
80 FullyDistVec<int64_t, int64_t> nonisoColV; // id's of non-isolated (connected) Col vertices
81 FullyDistVec<int64_t, int64_t> nonisov; // id's of non-isolated (connected) vertices
82
83 A.Reduce(*ColSums, Column, plus<int64_t>(), static_cast<int64_t>(0));
84 A.Reduce(*RowSums, Row, plus<int64_t>(), static_cast<int64_t>(0));
85
86 // this steps for general graph
87 /*
88 ColSums->EWiseApply(*RowSums, plus<int64_t>()); not needed for bipartite graph
89 nonisov = ColSums->FindInds(bind2nd(greater<int64_t>(), 0));
90 nonisov.RandPerm(); // so that A(v,v) is load-balanced (both memory and time wise)
91 A.operator()(nonisov, nonisov, true); // in-place permute to save memory
92 */
93
94 // this steps for bipartite graph
95 nonisoColV = ColSums->FindInds(bind2nd(greater<int64_t>(), 0));
96 nonisoRowV = RowSums->FindInds(bind2nd(greater<int64_t>(), 0));
97 delete ColSums;
98 delete RowSums;
99
100
101 {
102 nonisoColV.RandPerm();
103 nonisoRowV.RandPerm();
104 }
105
106
107 int64_t nrows1=A.getnrow(), ncols1=A.getncol(), nnz1 = A.getnnz();
108 double avgDeg1 = (double) nnz1/(nrows1+ncols1);
109
110
111 A.operator()(nonisoRowV, nonisoColV, true);
112
113 int64_t nrows2=A.getnrow(), ncols2=A.getncol(), nnz2 = A.getnnz();
114 double avgDeg2 = (double) nnz2/(nrows2+ncols2);
115
116
117 if(myrank == 0)
118 {
119 cout << "ncol nrows nedges deg \n";
120 cout << nrows1 << " " << ncols1 << " " << nnz1 << " " << avgDeg1 << " \n";
121 cout << nrows2 << " " << ncols2 << " " << nnz2 << " " << avgDeg2 << " \n";
122 }
123
124 MPI_Barrier(MPI_COMM_WORLD);
125
126
127}
128
129
130
131
132
133
134
135void ShowUsage()
136{
137 int myrank;
138 MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
139 if(myrank == 0)
140 {
141 cout << "\n-------------- usage --------------\n";
142 cout << "Usage (random matrix): ./maximal <er|g500|ssca> <Scale> <EDGEFACTOR> <algo><rand><moreSplit>\n";
143 cout << "Usage (input matrix): ./maximal <input> <matrix> <algo><rand><moreSplit>\n\n";
144
145 cout << " \n-------------- meaning of arguments ----------\n";
146 cout << "** er: Erdos-Renyi, g500: Graph500 benchmark, ssca: SSCA benchmark\n";
147 cout << "** scale: matrix dimention is 2^scale\n";
148 cout << "** edgefactor: average degree of vertices\n";
149 cout << "** algo : maximal matching algorithm used to initialize\n ";
150 cout << " greedy: greedy init , ks: Karp-Sipser, dmd: dynamic mindegree\n";
151 cout << " default: dynamic mindegree\n";
152 cout << "** (optional) rand: random parent selection in greedy/Karp-Sipser\n" ;
153 cout << "** (optional) moreSplit: more splitting of Matrix.\n" ;
154 cout << "(order of optional arguments does not matter)\n";
155
156 cout << " \n-------------- examples ----------\n";
157 cout << "Example: mpirun -np 4 ./maximal g500 18 16 ks rand" << endl;
158 cout << "Example: mpirun -np 4 ./maximal input cage12.mtx dmd\n" << endl;
159 }
160}
161
162void GetOptions(char* argv[], int argc)
163{
164 string allArg="";
165 for(int i=0; i<argc; i++)
166 {
167 allArg += string(argv[i]);
168 }
169
170 if(allArg.find("moreSplit")!=string::npos)
171 moreSplit = true;
172 if(allArg.find("randMaximal")!=string::npos)
173 randMaximal = true;
174 if(allArg.find("greedy")!=string::npos)
175 init = GREEDY;
176 else if(allArg.find("ks")!=string::npos)
177 init = KARP_SIPSER;
178 else if(allArg.find("dmd")!=string::npos)
179 init = DMD;
180 else
181 init = DMD;
182
183}
184
185void showCurOptions()
186{
187 ostringstream tinfo;
188 tinfo.str("");
189 tinfo << "\n---------------------------------\n";
190 tinfo << " Maximal matching algorithm options: ";
191 if(init == KARP_SIPSER) tinfo << " Karp-Sipser, ";
192 if(init == DMD) tinfo << " dynamic mindegree, ";
193 if(init == GREEDY) tinfo << " greedy, ";
194 if(randMaximal) tinfo << " random parent selection in greedy/Karp-Sipser, ";
195 if(moreSplit) tinfo << " moreSplit ";
196 tinfo << "\n---------------------------------\n\n";
197 SpParHelper::Print(tinfo.str());
198
199}
200
201void experiment( PSpMat_Bool & A, PSpMat_Bool & AT, FullyDistVec<int64_t, int64_t> degCol)
202{
203 FullyDistVec<int64_t, int64_t> mateRow2Col ( A.getcommgrid(), A.getnrow(), (int64_t) -1);
204 FullyDistVec<int64_t, int64_t> mateCol2Row ( A.getcommgrid(), A.getncol(), (int64_t) -1);
205
206 // best option
207 init = DMD; randMaximal = false;
208 //showCurOptions();
209 MaximalMatching(A, AT, mateRow2Col, mateCol2Row, degCol, init, randMaximal);
210 mateRow2Col.Apply([](int64_t val){return (int64_t) -1;});
211 mateCol2Row.Apply([](int64_t val){return (int64_t) -1;});
212
213 // best option + KS
214 init = KARP_SIPSER; randMaximal = true;
215 //showCurOptions();
216 MaximalMatching(A, AT, mateRow2Col, mateCol2Row, degCol, init, randMaximal);
217 mateRow2Col.Apply([](int64_t val){return (int64_t) -1;});
218 mateCol2Row.Apply([](int64_t val){return (int64_t) -1;});
219
220
221 // best option + Greedy
222 init = GREEDY; randMaximal = true;
223 //showCurOptions();
224 MaximalMatching(A, AT, mateRow2Col, mateCol2Row, degCol, init, randMaximal);
225 mateRow2Col.Apply([](int64_t val){return (int64_t) -1;});
226 mateCol2Row.Apply([](int64_t val){return (int64_t) -1;});
227
228 // best option + KS
229 init = KARP_SIPSER; randMaximal = false;
230 //showCurOptions();
231 MaximalMatching(A, AT, mateRow2Col, mateCol2Row, degCol, init, randMaximal);
232 mateRow2Col.Apply([](int64_t val){return (int64_t) -1;});
233 mateCol2Row.Apply([](int64_t val){return (int64_t) -1;});
234
235
236 // best option + Greedy
237 init = GREEDY; randMaximal = false;
238 //showCurOptions();
239 MaximalMatching(A, AT, mateRow2Col, mateCol2Row, degCol, init, randMaximal);
240 mateRow2Col.Apply([](int64_t val){return (int64_t) -1;});
241 mateCol2Row.Apply([](int64_t val){return (int64_t) -1;});
242
243}
244
245
246
247
248
249int main(int argc, char* argv[])
250{
251
252 // ------------ initialize MPI ---------------
253 int provided;
254 MPI_Init_thread(&argc, &argv, MPI_THREAD_SERIALIZED, &provided);
255 if (provided < MPI_THREAD_SERIALIZED)
256 {
257 printf("ERROR: The MPI library does not have MPI_THREAD_SERIALIZED support\n");
258 MPI_Abort(MPI_COMM_WORLD, 1);
259 }
260 int nprocs, myrank;
261 MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
262 MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
263 if(argc < 3)
264 {
265 ShowUsage();
266 MPI_Finalize();
267 return -1;
268 }
269
270 init = DMD;
271 randMaximal = false;
272 moreSplit = false;
273
274 // ------------ Process input arguments and build matrix ---------------
275 {
276
277 PSpMat_Bool * ABool;
278 ostringstream tinfo;
279 double t01, t02;
280 if(string(argv[1]) == string("input")) // input option
281 {
282 ABool = new PSpMat_Bool();
283
284 string filename(argv[2]);
285 tinfo.str("");
286 tinfo << "**** Reading input matrix: " << filename << " ******* " << endl;
287 SpParHelper::Print(tinfo.str());
288 t01 = MPI_Wtime();
289 ABool->ParallelReadMM(filename, true, maximum<double>());
290 t02 = MPI_Wtime();
291 ABool->PrintInfo();
292 tinfo.str("");
293 tinfo << "Reader took " << t02-t01 << " seconds" << endl;
294 SpParHelper::Print(tinfo.str());
295 GetOptions(argv+3, argc-3);
296
297 }
298 else if(argc < 4)
299 {
300 ShowUsage();
301 MPI_Finalize();
302 return -1;
303 }
304 else
305 {
306 unsigned scale = (unsigned) atoi(argv[2]);
307 unsigned EDGEFACTOR = (unsigned) atoi(argv[3]);
308 double initiator[4];
309 if(string(argv[1]) == string("er"))
310 {
311 initiator[0] = .25;
312 initiator[1] = .25;
313 initiator[2] = .25;
314 initiator[3] = .25;
315 }
316 else if(string(argv[1]) == string("g500"))
317 {
318 initiator[0] = .57;
319 initiator[1] = .19;
320 initiator[2] = .19;
321 initiator[3] = .05;
322 }
323 else if(string(argv[1]) == string("ssca"))
324 {
325 initiator[0] = .6;
326 initiator[1] = .4/3;
327 initiator[2] = .4/3;
328 initiator[3] = .4/3;
329 }
330 else
331 {
332 if(myrank == 0)
333 printf("The input type - %s - is not recognized.\n", argv[2]);
334 MPI_Abort(MPI_COMM_WORLD, 1);
335 }
336
337 SpParHelper::Print("Generating input matrix....\n");
338 t01 = MPI_Wtime();
339 DistEdgeList<int64_t> * DEL = new DistEdgeList<int64_t>();
340 DEL->GenGraph500Data(initiator, scale, EDGEFACTOR, true, true);
341 ABool = new PSpMat_Bool(*DEL, false);
342 delete DEL;
343 t02 = MPI_Wtime();
344 ABool->PrintInfo();
345 tinfo.str("");
346 tinfo << "Generator took " << t02-t01 << " seconds" << endl;
347 SpParHelper::Print(tinfo.str());
348
349 //Symmetricize(*ABool);
350 //removeIsolated(*ABool);
351 //SpParHelper::Print("Generated matrix symmetricized....\n");
352 ABool->PrintInfo();
353
354 GetOptions(argv+4, argc-4);
355
356 }
357
358
359 // randomly permute for load balance
360 SpParHelper::Print("Performing random permuation of matrix.\n");
361 FullyDistVec<int64_t, int64_t> prow(ABool->getcommgrid());
362 FullyDistVec<int64_t, int64_t> pcol(ABool->getcommgrid());
363 prow.iota(ABool->getnrow(), 0);
364 pcol.iota(ABool->getncol(), 0);
365 prow.RandPerm();
366 pcol.RandPerm();
367 (*ABool)(prow, pcol, true);
368 SpParHelper::Print("Performed random permuation of matrix.\n");
369
370
371 PSpMat_Bool A = *ABool;
372 PSpMat_Bool AT = A;
373 if(ABool->getnrow() > ABool->getncol())
374 AT.Transpose();
375 else
376 A.Transpose();
377
378
379 // Reduce is not multithreaded, so I am doing it here
380 FullyDistVec<int64_t, int64_t> degCol(A.getcommgrid());
381 A.Reduce(degCol, Column, plus<int64_t>(), static_cast<int64_t>(0));
382
383 int nthreads;
384#ifdef _OPENMP
385#pragma omp parallel
386 {
387 int splitPerThread = 1;
388 if(moreSplit) splitPerThread = 4;
389 nthreads = omp_get_num_threads();
390 cblas_splits = nthreads*splitPerThread;
391 }
392 tinfo.str("");
393 tinfo << "Threading activated with " << nthreads << " threads, and matrix split into "<< cblas_splits << " parts" << endl;
394 SpParHelper::Print(tinfo.str());
395 A.ActivateThreading(cblas_splits); // note: crash on empty matrix
396 AT.ActivateThreading(cblas_splits);
397#endif
398
399 SpParHelper::Print(" #####################################################\n");
400 SpParHelper::Print(" ################## Run 1 ############################\n");
401 SpParHelper::Print(" #####################################################\n");
402 experiment(A, AT, degCol);
403
404 SpParHelper::Print(" #####################################################\n");
405 SpParHelper::Print(" ################## Run 2 ############################\n");
406 SpParHelper::Print(" #####################################################\n");
407 experiment(A, AT, degCol);
408
409
410 SpParHelper::Print(" #####################################################\n");
411 SpParHelper::Print(" ################## Run 3 ############################\n");
412 SpParHelper::Print(" #####################################################\n");
413 experiment(A, AT, degCol);
414 }
415 MPI_Finalize();
416 return 0;
417}
418
419
moreSplit
bool moreSplit
Definition BPMaximalMatching.cpp:26
main
int main(int argc, char *argv[])
Definition BPMaximalMatching.cpp:249
mvInvertMate
bool mvInvertMate
Definition BPMaximalMatching.cpp:26
PSpMat_Bool
SpParMat< int64_t, bool, SpDCCols< int64_t, bool > > PSpMat_Bool
Definition BPMaximalMatching.cpp:64
randMaximal
bool randMaximal
Definition BPMaximalMatching.cpp:28
fewexp
bool fewexp
Definition BPMaximalMatching.cpp:29
prune
bool prune
Definition BPMaximalMatching.cpp:26
init
int init
Definition BPMaximalMatching.cpp:27
GetOptions
void GetOptions(char *argv[], int argc)
Definition BPMaximalMatching.cpp:162
ShowUsage
void ShowUsage()
Definition BPMaximalMatching.cpp:135
randMM
bool randMM
Definition BPMaximalMatching.cpp:26
showCurOptions
void showCurOptions()
Definition BPMaximalMatching.cpp:185
experiment
void experiment(PSpMat_Bool &A, PSpMat_Bool &AT, FullyDistVec< int64_t, int64_t > degCol)
Definition BPMaximalMatching.cpp:201
removeIsolated
void removeIsolated(PSpMat_Bool &A)
Definition BPMaximalMatching.cpp:69
GREEDY
#define GREEDY
Definition BPMaximalMatching.h:14
KARP_SIPSER
#define KARP_SIPSER
Definition BPMaximalMatching.h:15
DMD
#define DMD
Definition BPMaximalMatching.h:16
cblas_splits
int cblas_splits
Definition DirOptBFS.cpp:72
EDGEFACTOR
#define EDGEFACTOR
Definition DirOptBFS.cpp:81
combblas::DistEdgeList::GenGraph500Data
void GenGraph500Data(double initiator[4], int log_numverts, int edgefactor, bool scramble=false, bool packed=false)
Definition DistEdgeList.cpp:223
combblas::SpParHelper::Print
static void Print(const std::string &s)
Definition SpParHelper.cpp:811
int64_t
long int64_t
Definition compat.h:21
combblas::Column
@ Column
Definition SpDefs.h:114
combblas::Symmetricize
void Symmetricize(PARMAT &A)
Definition ReadMatDist.h:29
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
A
double A
VertexType::operator<
friend bool operator<(const VertexType &vtx1, const VertexType &vtx2)
Definition BPMaximalMatching.cpp:47
VertexType::operator==
friend bool operator==(const VertexType &vtx1, const VertexType &vtx2)
Definition BPMaximalMatching.cpp:52
VertexType::VertexType
VertexType(int64_t p=-1, int64_t r=-1, int16_t pr=0)
Definition BPMaximalMatching.cpp:45
VertexType::operator<<
friend ostream & operator<<(ostream &os, const VertexType &vertex)
Definition BPMaximalMatching.cpp:53