Zoltan2
Zoltan2_MatcherHelper.hpp
Go to the documentation of this file.
1 #ifndef _ZOLTAN2_MatcherHelper_hpp_
2 #define _ZOLTAN2_MatcherHelper_hpp_
3 
4 
5 #ifdef ZOLTAN2ND_HAVE_OMP
6 #include <omp.h>
7 #endif
8 
9 #include <iostream>
10 #include <fstream>
11 #include <string>
12 #include <vector>
13 #include <algorithm>
14 #include <cstdlib>
15 #include <set>
16 #include <queue>
17 #include <cassert>
18 
19 // PPF Matching code copied from Isorropia. Siva has reservations about algorithm
20 // robustness. It uses recursion, which can cause problems on the stack.
21 // Probably should rewrite at least some of these algorithms, eventually.
22 
23 
24 namespace Zoltan2 {
25 
34 class Matcher
36 {
37 private:
38 
39  int *mCRS_rowPtrs;
40  int *mCRS_cols;
41 
42  // Number of vertices in u set (num row vertices)
43  int numU;
44 
45  // Number of vertices in v set (num col vertices)
46  int numV;
47 
48  // For each u vertex, the matching v vertex
49  std::vector<int> matchedVertexForU;
50 
51  // For each v vertex, the matching u vertex
52  std::vector<int> matchedVertexForV;
53 
54  std::vector<int> queue;
55  int* LV_;
56  int* LU_;
57  int* unmatchedU_;
58  int *parent_;
59  int *lookahead_;
60  bool finish_;
61  int numE,avgDegU_,k_star_,icm_,BFSInd_,numThread_,Qst_,Qend_,numMatched;
62 
63  int *visitedV_;
64 
65  void delete_matched_v();
66  int SGM();
67  int match_dfs();
68  int match_hk();
69  int construct_layered_graph();
70  int find_set_del_M();
71  int recursive_path_finder(int, int);
72  int dfs_path_finder(int);
73  int dfs_augment();
74  int augment_matching(int);
75  int DW_phase();
76 
77 public:
85  Matcher(int *_rowPtr, int *_cols, int _numU, int _numV, int _numE);
86 
87 
88 
89 
93  virtual ~Matcher();
94 
95  /* @ingroup matching_grp
96  Returns the number of matched vertices from Maximum Cardinality
97  Matching set
98  */
100  {
101  int count=0;
102  for(unsigned int i=0;i<matchedVertexForU.size(); i++)
103  {
104  if(matchedVertexForU[i]!=-1)
105  {
106  count++;
107  }
108  }
109  return count;
110  }
111 
112 
113 
114  const std::vector<int> &getVertexUMatches() {return matchedVertexForU;};
115  const std::vector<int> &getVertexVMatches() {return matchedVertexForV;};
116 
120  int match();
121 };
123 
126 Matcher::Matcher(int *_rowPtr, int *_cols, int _numU, int _numV, int _numE)
127  :mCRS_rowPtrs(_rowPtr),mCRS_cols(_cols),numU(_numU),numV(_numV),numE(_numE)
128 {
129  finish_=false;
130  avgDegU_=numE/numU+1;
131  matchedVertexForU.resize(numU);
132  matchedVertexForV.resize(numV);
133 
134  lookahead_=new int[numU];
135  unmatchedU_=new int[numU];
136 
137  for(int i=0;i<numU;i++)
138  {
139  matchedVertexForU[i]=-1;
140 
141  lookahead_[i]=mCRS_rowPtrs[i];
142  unmatchedU_[i]=i;
143  }
144 
145  visitedV_=new int[numV];
146 
147  parent_=new int[numV];
148 
149  for(int i=0;i<numV;i++)
150  {
151  visitedV_[i]=0;
152 
153  matchedVertexForV[i]=-1;
154  parent_[i]=-1;
155  }
156 
157  numThread_=1;
158 }
160 
161 
165 {
166  delete [] lookahead_;
167  delete [] unmatchedU_;
168 
169  if (parent_) delete [] parent_; parent_ = NULL;
170 
171  if(visitedV_)
172  {
173  delete [] visitedV_;
174  }
175 
176 }
178 
180 // This function increases the matching size by one with the help of a
181 // augmenting path. The input is an integer which is the id of the last
182 // columns vertex of the augmenting path. We trace the whole path by
183 // backtraking using parent array. while tracing back we flip the unmatched
184 // edges to matched edges.
186 int Matcher::augment_matching(int tv)
187 {
188  int u,v,t,lnt=1;
189  v=tv;
190  while(true)
191  {
192  u=parent_[v];
193  t=matchedVertexForU[u];
194  matchedVertexForV[v]=u;
195  matchedVertexForU[u]=v;
196  if(t==-1)
197  break;
198  lnt+=2;
199  v=t;
200  }
201 
202  return lnt;
203 }
205 
207 // This function is almost similar to the previous function which is to find
208 // a vertex disjoint path. It is used for the algorithm DFS, PPF and in HKDW. The
209 // difference is that this function operates on the original graph not on the
210 // layerd subgraph. It also does the incorporates the lookahead mechanism and
211 // scanning adjacency list alternately from backward and from forward in
212 // alternate iterations for PPF and HKDW.
214 int Matcher::dfs_path_finder(int u)
215 {
216  int i,ind=-1,res=0;
217 
218  for(i=lookahead_[u];i<mCRS_rowPtrs[u+1];i++) // the lookahead scheme
219  {
220  ind=mCRS_cols[i];
221  assert(ind>=0 && ind<numV);
222  if(matchedVertexForV[ind]==-1)
223  {
224  int lock2=0;
225  if (visitedV_[ind]==0)
226  {
227  visitedV_[ind]=1;
228  lock2=1;
229  }
230 
231  if(lock2>0)
232  {
233  parent_[ind]=u;
234  lookahead_[u]=i+1;
235  return ind;
236  }
237  }
238  }
239 
240 
241  if(icm_%2==1) // odd number iteration so scan the adj list forward dir
242  {
243  for(i=mCRS_rowPtrs[u];i<mCRS_rowPtrs[u+1];i++)
244  {
245  ind=mCRS_cols[i];
246  assert(ind>=0 && ind<numV);
247 
248  int lock2=0;
249  if (visitedV_[ind]==0)
250  {
251  visitedV_[ind]=1;
252  lock2=1;
253  }
254 
255  if(lock2>0)
256  {
257  parent_[ind]=u;
258  res=dfs_path_finder(matchedVertexForV[ind]);
259  if(res!=-1)
260  return res;
261  }
262  }
263  }
264  else // even number iteration so scan from backward
265  {
266  for(i=mCRS_rowPtrs[u+1]-1;i>=mCRS_rowPtrs[u];i--)
267  {
268  ind=mCRS_cols[i];
269  assert(ind>=0 && ind<numV);
270 
271 
272  int lock2=0;
273  if (visitedV_[ind]==0)
274  {
275  visitedV_[ind]=1;
276  lock2=1;
277  }
278 
279  if(lock2>0)
280  {
281  parent_[ind]=u;
282  res=dfs_path_finder(matchedVertexForV[ind]);
283  if(res!=-1)
284  return res;
285  }
286  }
287  }
288 
289 
290  return -1;
291 }
293 
294 // ////////////////////////////////////////////////////////////////////////////////
295 // // This function starts the BFS phase for the HK and HKDW. It starts from the
296 // // layer 0 vertices and tries to find a vertex disjoint path from each of
297 // // these vertices.
298 // ////////////////////////////////////////////////////////////////////////////////
299 // int Matcher::find_set_del_M()
300 // {
301 
302 // int i,j,count=0;
303 // delete_matched_v();
304 
305 // #ifdef ZOLTAN2ND_HAVE_OMP
306 // #pragma omp parallel for private(j)
307 // #endif
308 // for(i=0;i<BFSInd_;i++)
309 // {
310 
311 // j=recursive_path_finder(0,queue[i]);
312 // if(j!=-1)
313 // {
314 // augment_matching(j);
315 
316 // //MMW: Not 100% this is necessary, let this in just in case bug in original code
317 // #ifdef ZOLTAN2ND_HAVE_OMP
318 // #pragma omp atomic
319 // #endif
320 // count++;
321 
322 
323 
324 // }
325 
326 // }
327 // return count;
328 // }
329 // ////////////////////////////////////////////////////////////////////////////////
330 
331 // ////////////////////////////////////////////////////////////////////////////////
332 // // This is the additional Duff and Wiberg phase for the HKDW. This function
333 // // does nothing but first unset the locks and then runs PPF from the
334 // // remaining unmatched row vertices after the BFS phase of HK.
335 // ////////////////////////////////////////////////////////////////////////////////
336 // int Matcher::DW_phase()
337 // {
338 // int i,count=0;
339 
340 // #ifdef ZOLTAN2ND_HAVE_OMP
341 // #pragma omp parallel for
342 // #endif
343 // for(i=0;i<numV;i++)
344 // {
345 // #ifdef ZOLTAN2ND_HAVE_OMP
346 // omp_unset_lock(&scannedV_[i]); // unsetting the locks
347 // #endif
348 // }
349 
350 // #ifdef ZOLTAN2ND_HAVE_OMP
351 // #pragma omp parallel for
352 // #endif
353 // for(i=0;i<BFSInd_;i++) // calling the PPF
354 // {
355 // int u=queue[i];
356 // if(matchedVertexForU[u]==-1)
357 // {
358 // int ind=dfs_path_finder(u);
359 // if(ind!=-1)
360 // {
361 // augment_matching(ind);
362 
363 // //MMW: Not 100% this is necessary, let this in just in case bug in original code
364 // #ifdef ISORROPIA_HAVE_OMP
365 // #pragma omp atomic
366 // #endif
367 // count++;
368 
369 // }
370 // }
371 // }
372 // return count;
373 // }
374 // ////////////////////////////////////////////////////////////////////////////////
375 
377  //This function is the starter function for PPF and DFS. It unsets the locks
378  //the call dfs_path_finder and then call the augment_matching.
380 int Matcher::dfs_augment()
381 {
382  int i,flag=0,flag1=0,count=0,totc=0,index=numU,cur=0;
383  icm_=0;
384 
385  while(true)
386  {
387  icm_++;
388 
389  for(i=0;i<numV;i++)
390  {
391  visitedV_[i]=0;
392  }
393 
394  cur=0;
395  for(i=0;i<numU;i++)
396  {
397  if(matchedVertexForU[i]==-1)
398  unmatchedU_[cur++]=i;
399  }
400  index=cur;
401  flag=flag1=count=0;
402 
403  for(i=0;i<index;i++)
404  {
405  flag=1;
406  int u=unmatchedU_[i];
407  int ind=dfs_path_finder(u);
408  if(ind!=-1)
409  {
410  flag1=1;
411  augment_matching(ind);
412  }
413  }
414 
415  if(flag==0 || flag1==0)
416  break;
417  else
418  {
419  cur=0;
420  for(i=0;i<index;i++)
421  {
422  if(matchedVertexForU[unmatchedU_[i]]==-1)
423  unmatchedU_[cur++]=unmatchedU_[i];
424  }
425  index=cur;
426 
427  }
428  }
429 
430  return totc;
431 }
433 
434 // ////////////////////////////////////////////////////////////////////////////////
435 // ////////////////////////////////////////////////////////////////////////////////
436 // int Matcher::SGM()
437 // {
438 // int i,j,lock,ind,count=0;
439 // #ifdef ZOLTAN2ND_HAVE_OMP
440 // #pragma omp parallel for private(j,ind,lock)
441 // #endif
442 // for(i=0;i<numU;i++)
443 // {
444 // for(j=mCRS_rowPtrs[i];j<mCRS_rowPtrs[i+1];j++)
445 // {
446 // ind=mCRS_cols[j];
447 // #ifdef ZOLTAN2ND_HAVE_OMP
448 // lock=omp_test_lock(&scannedV_[ind]);
449 // #else
450 // // mfh 07 Feb 2013: lock wasn't getting initialized if
451 // // ZOLTAN2ND_HAVE_OMP was not defined. omp_test_lock()
452 // // returns nonzero if the thread successfully acquired the
453 // // lock. If there's only one thread, that thread will
454 // // always be successful, so the default value of lock
455 // // should be nonzero.
456 // lock = 1;
457 // #endif
458 // if(lock>0)
459 // {
460 // matchedVertexForU[i]=ind;
461 // matchedVertexForV[ind]=i;
462 // break;
463 // }
464 // }
465 // }
466 // return count;
467 // }
468 // ////////////////////////////////////////////////////////////////////////////////
469 
470 
473 int Matcher::SGM()
474 {
475  int i,j,ind,count=0;
476 
477  for(i=0;i<numU;i++)
478  {
479  for(j=mCRS_rowPtrs[i];j<mCRS_rowPtrs[i+1];j++)
480  {
481  ind=mCRS_cols[j];
482 
483  int lock2=0;
484  if (visitedV_[ind]==0)
485  {
486  visitedV_[ind]=1;
487  lock2=1;
488  }
489 
490  if(lock2>0)
491  {
492  matchedVertexForU[i]=ind;
493  matchedVertexForV[ind]=i;
494  break;
495  }
496  }
497  }
498  return count;
499 }
501 
502 
503 
505 // Forking function for DFS based algorithm
507 int Matcher::match_dfs()
508 {
509  int totc=0;
510  icm_=0;
511 
512  totc=dfs_augment();
513 
514  return totc;
515 }
517 
521 {
522  // User interface function for the matching..
523 
524  int totc=0;
525  totc=SGM();
526  totc+=match_dfs();
527 
528  numMatched = totc;
529 
530  return 0;
531 }
533 
534 
536 // Calculate vertex cover (which is vertex separator) from matching
537 // VC = (U-L) union (B intersection L)
538 //
539 // Let X be the exposed nodes in U (those without a match)
540 //
541 // E':
542 // Matched edges should be directed VtoU -- just use vertVMatches array
543 // All other edges should be directed UtoV -- use modified biGraphCRScols
544 //
545 // L is the set of vertices in (U union V, E') that can be reached from X
546 //
547 //
548 // Perhaps I should use internal class members in this function
549 // However, I decided not to do this for now since I'm changing some of these
550 // member variables in this function
552 void getVCfromMatching(const std::vector<int> &bigraphCRSRowPtr,
553  std::vector<int> &bigraphCRSCols,
554  const std::vector<int> &vertUMatches,
555  const std::vector<int> &vertVMatches,
556  const std::vector<int> &bigraphVMapU,
557  const std::vector<int> &bigraphVMapV,
558  std::vector<int> &VC)
559 {
560  int numU = vertUMatches.size();
561  int numV = vertVMatches.size();
562 
564  // Build X
566  std::set<int> X;
567  for(int i=0;i<numU; i++)
568  {
569  if(vertUMatches[i]==-1)
570  {
571  X.insert(i);
572  }
573  }
575 
577  // Non matched edges should be directed UtoV
578  // Removed matches edges from bipartite graph (set to -1)
579  // -- perhaps replace this with something more efficient/compact
581  for (int uID=0;uID<numU;uID++)
582  {
583  for (int eIdx=bigraphCRSRowPtr[uID];eIdx<bigraphCRSRowPtr[uID+1];eIdx++)
584  {
585  int vID = bigraphCRSCols[eIdx];
586 
587  if (vertUMatches[uID]==vID)
588  {
589  bigraphCRSCols[eIdx]=-1;
590  }
591  }
592 
593  }
595 
597  // Calculate L - set of vertices in (U union V, E') that can be reached from X
599  std::set<int> L;
600 
601  std::queue<int> queue;
602 
603  std::copy(X.begin(), X.end(), std::inserter( L, L.begin() ) );
604 
605  for(std::set<int>::const_iterator iter = X.begin(); iter != X.end(); ++iter)
606  {
607  L.insert(bigraphVMapU[*iter]); // L contains all vertices in X
608 
609  queue.push(*iter); // copy on to queue
610  }
611 
612  int iteration=0;
613  while(queue.size()!=0)
614  {
615  //Number of active vertices on this side of bipartite graph
616  int nstarts=queue.size();
617 
618  for (int i=0; i<nstarts; i++)
619  {
620  int start = queue.front();
621  queue.pop();
622 
623  //Traverse from U to V
624  if(iteration%2==0)
625  {
626  //Traverse edges starting from vertex "start"
627  for (int eIdx=bigraphCRSRowPtr[start];eIdx<bigraphCRSRowPtr[start+1];eIdx++)
628  {
629  int newV = bigraphCRSCols[eIdx];
630 
631  //Edge is in correct direction (U to V) => newV is valid
632  if (newV!=-1)
633  {
634  // If this vertex has not been reached
635  if(L.find(bigraphVMapV[newV])==L.end())
636  {
637  L.insert(bigraphVMapV[newV]);
638  queue.push(newV);
639  }
640  }
641  }
642 
643  }
644 
645  //Traverse from V to U
646  else
647  {
648  int newU = vertVMatches[start];
649 
650  // If this vertex has not been reached
651  if(L.find(bigraphVMapU[newU])==L.end())
652  {
653  L.insert(bigraphVMapU[newU]);
654  queue.push(newU);
655  }
656  }
657  } // for
658 
659  iteration++;
660  } //while
662 
664  // Calc VC = (U-L) union (B intersection L)
666  for(int uID=0;uID<numU;uID++)
667  {
668  // if vertex not in L, it is in VC
669  if(L.find(bigraphVMapU[uID])==L.end())
670  {
671  VC.push_back(bigraphVMapU[uID]);
672  }
673  }
674 
675  for(int vID=0;vID<numV;vID++)
676  {
677  // if vertex is in L, it is in VC
678  if(L.find(bigraphVMapV[vID])!=L.end())
679  {
680  VC.push_back(bigraphVMapV[vID]);
681  }
682  }
684 
685 
686 
687 }
689 
690 
691 
692 
693 
694 
695 
696 
697 
698 
699 
700 
701 
702 } //Zoltan2 namespace
703 #endif //ifdef
void getVCfromMatching(const std::vector< int > &bigraphCRSRowPtr, std::vector< int > &bigraphCRSCols, const std::vector< int > &vertUMatches, const std::vector< int > &vertVMatches, const std::vector< int > &bigraphVMapU, const std::vector< int > &bigraphVMapV, std::vector< int > &VC)
int match()
Computes the maximum cardinality matching.
const std::vector< int > & getVertexVMatches()
virtual ~Matcher()
Destructor.
Matcher(int *_rowPtr, int *_cols, int _numU, int _numV, int _numE)
Constructor.
const std::vector< int > & getVertexUMatches()