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