Zoltan2
Zoltan2_AlgND.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Zoltan2: A package of combinatorial algorithms for scientific computing
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Karen Devine (kddevin@sandia.gov)
39 // Erik Boman (egboman@sandia.gov)
40 // Siva Rajamanickam (srajama@sandia.gov)
41 // Michael Wolf (mmwolf@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 
47 //MMW need to specify that this requires Zoltan
48 
49 #ifndef _ZOLTAN2_ALGND_HPP_
50 #define _ZOLTAN2_ALGND_HPP_
51 
54 #include <Zoltan2_Algorithm.hpp>
55 #include <Zoltan2_AlgZoltan.hpp>
56 
58 
59 
60 #include <sstream>
61 #include <string>
62 #include <bitset>
63 
69 void buildPartTree(int level, int leftPart, int splitPart, int rightPart, std::vector<int> &partTree);
70 
71 
72 namespace Zoltan2
73 {
74 
76 
87 template <typename Adapter>
89 class AlgND : public Algorithm<typename Adapter::base_adapter_t>
90 //class AlgND : public Algorithm<Adapter>
91 {
92 
93 private:
94 
95  typedef typename Adapter::part_t part_t;
96 
97  typedef typename Adapter::lno_t lno_t;
98  typedef typename Adapter::gno_t gno_t;
99 
100 
101  const RCP<const Environment> mEnv;
102  const RCP<const Comm<int> > mProblemComm;
103 
104  // const RCP<const GraphModel<Adapter> > mGraphModel;
105  const RCP<GraphModel<typename Adapter::base_adapter_t> > mGraphModel;
106  // const RCP<const CoordinateModel<Adapter> > mIds;
107  const RCP<CoordinateModel<typename Adapter::base_adapter_t> > mIds;
108 
109  //const RCP<const Adapter> mBaseInputAdapter;
110  //const RCP<const Adapter> mInputAdapter;
111  const RCP<const typename Adapter::base_adapter_t> mBaseInputAdapter;
112 
113  void getBoundLayer(int levelIndx, const std::vector<part_t> &partMap,
114  const part_t * parts,
115  const std::set<int> &excVerts,
116  int &bigraphNumS, int &bigraphNumT, int &bigraphNumE,
117  std::vector<int> &bigraphCRSRowPtr, std::vector<int> &bigraphCRSCols,
118  std::vector<int> &bigraphVMapU, std::vector<int> &bigraphVMapV);
119 
120 
121 public:
122  // Constructor
123  AlgND(const RCP<const Environment> &env_,
124  const RCP<const Comm<int> > &problemComm_,
127  const RCP<const typename Adapter::base_adapter_t> baseInputAdapter_
128  )
129  :mEnv(env_), mProblemComm(problemComm_), mGraphModel(gModel_),
130  mIds(cModel_), mBaseInputAdapter(baseInputAdapter_)
131  {
132 #ifndef INCLUDE_ZOLTAN2_EXPERIMENTAL
133  Z2_THROW_EXPERIMENTAL("Zoltan2 AlgND is strictly experimental software ")
134 #endif
135 
136 #ifndef INCLUDE_ZOLTAN2_EXPERIMENTAL_WOLF
137  Z2_THROW_EXPERIMENTAL_WOLF("Zoltan2 algND is strictly experimental software ")
138 #endif
139 
140  if(mProblemComm->getSize()!=1)
141  {
142  Z2_THROW_SERIAL("Zoltan2 AlgND is strictly serial!");
143  }
144 
145  }
146 
147  // Ordering method
148  int order(const RCP<OrderingSolution<lno_t, gno_t> > &solution_);
149 
150 };
152 
155 template <typename Adapter>
157 {
158  // typedef typename Adapter::lno_t lno_t; // local ids
159  // typedef typename Adapter::gno_t gno_t; // global ids
160  // typedef typename Adapter::scalar_t scalar_t; // scalars
161 
162  mEnv->debug(DETAILED_STATUS, std::string("Entering AlgND"));
163 
165  // First, let's partition with RCB using Zoltan. Eventually, we will change this
166  // to use PHG
168 
169  RCP<PartitioningSolution<Adapter> > partSoln;
170  int nUserWts=0;
171 
172  std::cout << "HERE1" << std::endl;
173 
174  partSoln =
175  RCP<PartitioningSolution<Adapter> > (new PartitioningSolution<Adapter>(this->mEnv, mProblemComm, nUserWts));
176 
177  AlgZoltan<Adapter> algZoltan(this->mEnv, mProblemComm, this->mBaseInputAdapter);
178 
179  std::cout << "HERE2" << std::endl;
180 
181  algZoltan.partition(partSoln);
182 
183  std::cout << "HERE3" << std::endl;
184 
185 
186  size_t numGlobalParts = partSoln->getTargetGlobalNumberOfParts();
187 
188  const part_t *parts = partSoln->getPartListView();
190 
192  // Build up tree that represents partitioning subproblems, which will
193  // be used for determining separators at each level
194  // -- for now, this is built up artificially
195  // -- eventually this will be obtained from PHG output
196  //
197  // Each separator i is represented by 4 integers/part_t? in the partTree
198  // structure: partTree[4*i], partTree[4*i+1], partTree[4*i+2], partTree[4*i+3]
199  // These 4 integers are level of separator, smallest part in 1st half of separator,
200  // smallest part in 2nd half of separator, largest part in 2nd half of separator + 1
202  // change int to something, part_t?
203 
204  std::cout << "HERE4" << std::endl;
205 
206 
207  std::vector<int> partTree;
208 
209  buildPartTree( 0, 0, (numGlobalParts-1)/2 + 1, numGlobalParts, partTree);
210  unsigned int numSeparators = partTree.size() / 4;
211 
212  for(unsigned int i=0;i<partTree.size(); i++)
213  {
214  std::cout << "partTree: " << partTree[i] << std::endl;
215  }
216  std::cout << "NumSeparators: " << numSeparators << std::endl;
217 
218  std::cout << "HERE5" << std::endl;
219 
221 
222 
224  // Create a map that maps each part number to a new number based on
225  // the level of the hiearchy of the separator tree. This allows us
226  // to easily identify the boundary value vertices
228  std::cout << "HERE6" << std::endl;
229 
230 
231  int numLevels = partTree[4*(numSeparators-1)]+1;
232 
233  std::vector<std::vector<int> > partLevelMap(numLevels,std::vector<int>(numGlobalParts));
234 
235  std::vector<int> sepsInLev(numLevels,0);
236 
237  for(unsigned int i=0;i<numSeparators;i++)
238  {
239  int level = partTree[4*i];
240  int leftPart = partTree[4*i+1];
241  int splitPart = partTree[4*i+2];
242  int rightPart = partTree[4*i+3];
243 
244  for(int part=leftPart; part<splitPart; part++)
245  {
246  partLevelMap[level][part] = 2*sepsInLev[level];
247  }
248 
249  for(int part=splitPart; part<rightPart; part++)
250  {
251  partLevelMap[level][part] = 2*sepsInLev[level]+1;
252  }
253 
254  sepsInLev[level]++;
255  }
256 
257  std::cout << "partLevelMap[0][0] = " << partLevelMap[0][0] << std::endl;
258  std::cout << "partLevelMap[0][1] = " << partLevelMap[0][1] << std::endl;
259 
260  std::cout << "HERE7" << std::endl;
261 
263 
264  // Set of separator vertices. Used to keep track of what vertices are
265  // already in previous calculated separators. These vertices should be
266  // excluded from future separator calculations
267  const std::set<int> sepVerts;
268 
270  // Loop over each cut
271  // 1. Build boundary layer between parts
272  // 2. Build vertex separator from boundary layer
274  std::cout << "HERE8" << std::endl;
275 
276  for(unsigned int level=0;level<numLevels;level++)
277  {
278  for(unsigned int levIndx=0;levIndx<sepsInLev[level];levIndx++)
279  {
281  // Build boundary layer between parts (edge separator)
283  std::cout << "HERE9" << std::endl;
284 
285  int bigraphNumU=0, bigraphNumV=0, bigraphNumE=0;
286  std::vector<int> bigraphVMapU;
287  std::vector<int> bigraphVMapV;
288 
289  std::vector<int> bigraphCRSRowPtr;
290  std::vector<int> bigraphCRSCols;
291 
292 
293  getBoundLayer(levIndx, partLevelMap[level], parts, sepVerts,
294  bigraphNumU,bigraphNumV,bigraphNumE,
295  bigraphCRSRowPtr, bigraphCRSCols,
296  bigraphVMapU,bigraphVMapV);
297 
298  std::cout << "Bipartite graph: " << bigraphNumU << " " << bigraphNumV << " "
299  << bigraphNumE << std::endl;
300 
301  for (unsigned int i=0;i<bigraphVMapU.size();i++)
302  {
303  std::cout << "boundVertU: " << bigraphVMapU[i] << std::endl;
304  }
305 
306  for (unsigned int i=0;i<bigraphVMapV.size();i++)
307  {
308  std::cout << "boundVertV: " << bigraphVMapV[i] << std::endl;
309  }
310 
311 
312 
313  for (int rownum=0;rownum<bigraphNumU;rownum++)
314  {
315 
316  for (int eIdx=bigraphCRSRowPtr[rownum];eIdx<bigraphCRSRowPtr[rownum+1];eIdx++)
317  {
318  std::cout << "bipartite E: " << bigraphVMapU[rownum] << ", " << bigraphVMapV[ bigraphCRSCols[eIdx]]
319  << " ( " << rownum << "," << bigraphCRSCols[eIdx] << " )" << std::endl;
320  }
321 
322  }
323  std::cout << "HERE10" << std::endl;
325 
327  // Calculate bipartite matching from boundary layer
329  Matcher bpMatch(bigraphCRSRowPtr.data(), bigraphCRSCols.data(), bigraphNumU, bigraphNumV, bigraphNumE);
330  bpMatch.match();
331 
332  const std::vector<int> &vertUMatches = bpMatch.getVertexUMatches();
333  const std::vector<int> &vertVMatches = bpMatch.getVertexVMatches();
335 
337  // Calculate vertex cover (which is vertex separator) from matching
339  std::vector<int> VC;
340 
341  getVCfromMatching(bigraphCRSRowPtr,bigraphCRSCols,vertUMatches,vertVMatches,
342  bigraphVMapU,bigraphVMapV,VC);
343 
344  for(unsigned int i=0;i<VC.size();i++)
345  {
346  std::cout << "VC: " << VC[i] << std::endl;
347  }
349 
350 
351 
352  }
353  }
354 
355  std::cout << "HERE20" << std::endl;
356 
358 
359  // //TODO: calculate vertex separator for each layer,
360  // //TODO: using vertex separators, compute new ordering and store in solution
361  // //TODO: move to ordering directory
362 
363  mEnv->debug(DETAILED_STATUS, std::string("Exiting AlgND"));
364  return 0;
365 }
367 
369 // Create boundary layer of vertices between 2 partitions
370 //
371 // Could improve the efficiency here by first creating a boundary layer graph
372 // between all parts
374 template <typename Adapter>
375 void AlgND<Adapter>::getBoundLayer(int levelIndx, const std::vector<part_t> &partMap,
376  const part_t * parts,
377  const std::set<int> &excVerts,
378  int &bigraphNumS, int &bigraphNumT, int &bigraphNumE,
379  std::vector<int> &bigraphCRSRowPtr, std::vector<int> &bigraphCRSCols,
380  std::vector<int> &bigraphVMapS, std::vector<int> &bigraphVMapT)
381 {
382  std::cout << "HI1" << std::endl;
383 
384  typedef typename Adapter::lno_t lno_t; // local ids
385  typedef typename Adapter::scalar_t scalar_t; // scalars
386  typedef StridedData<lno_t, scalar_t> input_t;
387 
388  int numVerts = mGraphModel->getLocalNumVertices();
389 
390  //Teuchos ArrayView
391  ArrayView< const lno_t > eIDs;
392  ArrayView< const lno_t > vOffsets;
393  ArrayView< input_t > wgts;
394 
395  // For some reason getLocalEdgeList seems to be returning empty eIDs
396  //size_t numEdges = ( (GraphModel<typename Adapter::base_adapter_t>) *mGraphModel).getLocalEdgeList(eIDs, vOffsets, wgts);
397 
398  //size_t numEdges = ( (GraphModel<typename Adapter::base_adapter_t>) *mGraphModel).getEdgeList(eIDs, vOffsets, wgts);
399  ( (GraphModel<typename Adapter::base_adapter_t>) *mGraphModel).getEdgeList(eIDs, vOffsets, wgts);
400 
401 
402  std::map<int,std::set<int> > bigraphEs;
403  std::set<int> vSetS;
404  std::set<int> vSetT;
405 
406  bigraphNumE=0;
407 
408  for(int v1=0;v1<numVerts;v1++)
409  {
410 
411  part_t vpart1 = partMap[parts[v1]];
412 
413  bool correctBL = (vpart1 >= 2*levelIndx && vpart1 < 2*(levelIndx+1) );
414 
415  // if vertex is not in the correct range of parts, it cannot be a member of
416  // this boundary layer
417  if(!correctBL)
418  {
419  continue;
420  }
421 
422  // Ignore vertices that belong to set of vertices to exclude
423  if(excVerts.find(v1)!=excVerts.end())
424  {
425  continue;
426  }
427 
428  //Loop over edges connected to v1
429  //MMW figure out how to get this from Zoltan2
430  for(int j=vOffsets[v1];j<vOffsets[v1+1];j++)
431  {
432 
433  int v2 = eIDs[j];
434 
435  part_t vpart2 = partMap[parts[v2]];
436 
437  correctBL = (vpart2 >= 2*levelIndx && vpart2 < 2*(levelIndx+1) );
438 
439  // if vertex is not in the correct range of parts, it cannot be a member of
440  // this boundary layer
441  if(!correctBL)
442  {
443  continue;
444  }
445 
446  // Ignore vertices that belong to set of vertices to exclude
447  if(excVerts.find(v2)!=excVerts.end())
448  {
449  continue;
450  }
451 
452  if ( vpart1 != vpart2 )
453  {
454  // Vertex added to 1st set of boundary vertices
455  if(vpart1<vpart2)
456  {
457  vSetS.insert(v1);
458 
459  // v1, v2
460  if(bigraphEs.find(v1)==bigraphEs.end())
461  {
462  bigraphEs[v1] = std::set<int>();
463  }
464  bigraphEs[v1].insert(v2);
465  bigraphNumE++;
466 
467  }
468  // Vertex added to 2nd set of boundary vertices
469  else
470  {
471  vSetT.insert(v1);
472  }
473  }
474 
475  }
476  }
477 
479  // Set size of two vertex sets for bipartite graph
481  bigraphNumS = vSetS.size();
482  bigraphNumT = vSetT.size();
484 
487 
488  bigraphVMapS.resize(bigraphNumS);
489 
490  std::map<int,int> glob2LocTMap;
491 
492  unsigned int indx=0;
493  for(std::set<int>::const_iterator iter=vSetS.begin(); iter!=vSetS.end(); ++iter)
494  {
495  bigraphVMapS[indx] = *iter;
496  indx++;
497  }
498 
499 
500  bigraphVMapT.resize(bigraphNumT);
501  indx=0;
502  for(std::set<int>::const_iterator iter=vSetT.begin();iter!=vSetT.end();++iter)
503  {
504  bigraphVMapT[indx] = *iter;
505  glob2LocTMap[*iter]=indx;
506  indx++;
507  }
509 
510 
512  // Set sizes for bipartite graph data structures
514  bigraphCRSRowPtr.resize(bigraphNumS+1);
515  bigraphCRSCols.resize(bigraphNumE);
517 
519  // Copy bipartite graph edges into CRS format
521  bigraphCRSRowPtr[0]=0;
522 
523  unsigned int rownum=0;
524  unsigned int nzIndx=0;
525  std::map<int,std::set<int> >::const_iterator iterM;
526  for (iterM=bigraphEs.begin();iterM!=bigraphEs.end();++iterM)
527  {
528  bigraphCRSRowPtr[rownum+1] = bigraphCRSRowPtr[rownum] + (*iterM).second.size();
529 
530  for(std::set<int>::const_iterator iter=(*iterM).second.begin(); iter!=(*iterM).second.end(); ++iter)
531  {
532  bigraphCRSCols[nzIndx] = glob2LocTMap[(*iter)];
533 
534  nzIndx++;
535  }
536 
537  rownum++;
538  }
540 
541 }
543 
544 
546 // Create boundary layer of vertices between 2 partitions
548 // template <typename Adapter>
549 // void AlgND<Adapter>::getBoundLayer(int levelIndx, const std::vector<part_t> &partMap,
550 // const part_t * parts,
551 // const std::set<int> &excVerts,
552 // std::vector<int> &boundVerts,
553 // std::vector<std::vector<int> > &boundVertsST)
554 
555 // {
556 // std::cout << "HI1" << std::endl;
557 
558 // typedef typename Adapter::lno_t lno_t; // local ids
559 // typedef typename Adapter::scalar_t scalar_t; // scalars
560 // typedef StridedData<lno_t, scalar_t> input_t;
561 
562 // int numVerts = mGraphModel->getLocalNumVertices();
563 
564 // std::cout << "NumVerts: " << numVerts << std::endl;
565 
566 // //Teuchos ArrayView
567 // ArrayView< const lno_t > eIDs;
568 // ArrayView< const lno_t > vOffsets;
569 // ArrayView< input_t > wgts;
570 
571 // // For some reason getLocalEdgeList seems to be returning empty eIDs
572 // //size_t numEdges = ( (GraphModel<typename Adapter::base_adapter_t>) *mGraphModel).getLocalEdgeList(eIDs, vOffsets, wgts);
573 
574 // size_t numEdges = ( (GraphModel<typename Adapter::base_adapter_t>) *mGraphModel).getEdgeList(eIDs, vOffsets, wgts);
575 
576 // for(int v1=0;v1<numVerts;v1++)
577 // {
578 
579 // part_t vpart1 = partMap[parts[v1]];
580 
581 // bool correctBL = (vpart1 >= 2*levelIndx && vpart1 < 2*(levelIndx+1) );
582 
583 // // if vertex is not in the correct range of parts, it cannot be a member of
584 // // this boundary layer
585 // if(!correctBL)
586 // {
587 // continue;
588 // }
589 
590 // // Ignore vertices that belong to set of vertices to exclude
591 // if(excVerts.find(v1)!=excVerts.end())
592 // {
593 // continue;
594 // }
595 
596 // //Loop over edges connected to v1
597 // //MMW figure out how to get this from Zoltan2
598 // for(int j=vOffsets[v1];j<vOffsets[v1+1];j++)
599 // {
600 
601 // int v2 = eIDs[j];
602 
603 // part_t vpart2 = partMap[parts[v2]];
604 
605 // correctBL = (vpart2 >= 2*levelIndx && vpart2 < 2*(levelIndx+1) );
606 
607 // // if vertex is not in the correct range of parts, it cannot be a member of
608 // // this boundary layer
609 // if(!correctBL)
610 // {
611 // continue;
612 // }
613 
614 // // Ignore vertices that belong to set of vertices to exclude
615 // if(excVerts.find(v2)!=excVerts.end())
616 // {
617 // continue;
618 // }
619 
620 // if ( vpart1 != vpart2 )
621 // {
622 // // Vertex added to set of all boundary vertices
623 // boundVerts.push_back(v1);
624 
625 // // Vertex added to 1st set of boundary vertices
626 // if(vpart1<vpart2)
627 // {
628 // boundVertsST[0].push_back(v1);
629 // }
630 // // Vertex added to 2nd set of boundary vertices
631 // else
632 // {
633 // boundVertsST[1].push_back(v1);
634 // }
635 // break;
636 // }
637 
638 // }
639 // }
640 
641 // }
643 
644 } // namespace Zoltan2
645 
646 
649 void buildPartTree(int level, int leftPart, int splitPart, int rightPart, std::vector<int> &partTree)
650 {
651  // Insert information for this separator
652  partTree.push_back(level);
653  partTree.push_back(leftPart);
654  partTree.push_back(splitPart);
655  partTree.push_back(rightPart);
656 
657  // Recurse down left side of tree
658  if(splitPart-leftPart > 1)
659  {
660  int newSplit = leftPart+(splitPart-leftPart-1)/2 + 1;
661  buildPartTree(level+1,leftPart,newSplit,splitPart,partTree);
662  }
663 
664  // Recurse down right side of tree
665  if(rightPart-splitPart>1)
666  {
667  int newSplit = splitPart+(rightPart-splitPart-1)/2 + 1;
668  buildPartTree(level+1,splitPart,newSplit,rightPart,partTree);
669  }
670 }
672 
673 
674 
675 
676 #endif
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.
#define Z2_THROW_SERIAL(mystr)
Throw an error when code is run on more than one processor.
#define Z2_THROW_EXPERIMENTAL_WOLF(mystr)
Throw an error when wolf experimental code is requested but not compiled.
interface to the Zoltan package
Defines the PartitioningSolution class.
An implementation of the Matcher interface that operates on Epetra matrices and Graphs.
sub-steps, each method&#39;s entry and exit
This class provides geometric coordinates with optional weights to the Zoltan2 algorithm.
void partition(const RCP< PartitioningSolution< Adapter > > &solution)
Partitioning method.
Defines the IdentifierModel interface.
void buildPartTree(int level, int leftPart, int splitPart, int rightPart, std::vector< int > &partTree)
A PartitioningSolution is a solution to a partitioning problem.
The StridedData class manages lists of weights or coordinates.
lno_t size() const
Return the length of the strided array.
Algorithm defines the base class for all algorithms.
#define Z2_THROW_EXPERIMENTAL(mystr)
Throw an error when experimental code is requested but not compiled.
GraphModel defines the interface required for graph models.
AlgND(const RCP< const Environment > &env_, const RCP< const Comm< int > > &problemComm_, const RCP< GraphModel< typename Adapter::base_adapter_t > > &gModel_, const RCP< CoordinateModel< typename Adapter::base_adapter_t > > &cModel_, const RCP< const typename Adapter::base_adapter_t > baseInputAdapter_)
The class containing ordering solutions.
int order(const RCP< OrderingSolution< lno_t, gno_t > > &solution_)
Ordering method.