Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Zoltan2_AlgZoltan.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 //
42 // ***********************************************************************
43 //
44 // @HEADER
45 #ifndef _ZOLTAN2_ALGZOLTAN_HPP_
46 #define _ZOLTAN2_ALGZOLTAN_HPP_
47 
48 #include <Zoltan2_Standards.hpp>
49 #include <Zoltan2_Algorithm.hpp>
51 #include <Zoltan2_Util.hpp>
52 #include <Zoltan2_TPLTraits.hpp>
53 
54 #include <Zoltan2_Model.hpp>
55 
57 #include <zoltan_cpp.h>
58 #include <zoltan_partition_tree.h>
59 
63 //
64 // This first design templates Zoltan's callback functions on the
65 // input adapter. This approach has the advantage of simplicity and
66 // is most similar to current usage of Zoltan (where the callbacks define
67 // the model).
68 // A better approach might template them on a model,
69 // allowing Zoltan2 greater flexibility in creating models from the input.
70 // Alternatively, different callback implementations could be provided to
71 // represent different models to Zoltan.
73 
74 namespace Zoltan2 {
75 
76 template <typename Adapter>
77 class AlgZoltan : public Algorithm<Adapter>
78 {
79 
80 private:
81 
82  typedef typename Adapter::lno_t lno_t;
83  typedef typename Adapter::gno_t gno_t;
84  typedef typename Adapter::scalar_t scalar_t;
85  typedef typename Adapter::part_t part_t;
86  typedef typename Adapter::user_t user_t;
87  typedef typename Adapter::userCoord_t userCoord_t;
88 
89  const RCP<const Environment> env;
90  const RCP<const Comm<int> > problemComm;
91  const RCP<const typename Adapter::base_adapter_t> adapter;
92  RCP<const Model<Adapter> > model;
93  RCP<Zoltan> zz;
94 
95  MPI_Comm mpicomm;
96 
97  void setMPIComm(const RCP<const Comm<int> > &problemComm__) {
98 # ifdef HAVE_ZOLTAN2_MPI
99  mpicomm = Teuchos::getRawMpiComm(*problemComm__);
100 # else
101  mpicomm = MPI_COMM_WORLD; // taken from siMPI
102 # endif
103  }
104 
105  void zoltanInit() {
106  // call Zoltan_Initialize to make sure MPI_Init is called (in MPI or siMPI).
107  int argc = 0;
108  char **argv = NULL;
109  float ver;
110  Zoltan_Initialize(argc, argv, &ver);
111  }
112 
113  void setCallbacksIDs()
114  {
115  zz->Set_Num_Obj_Fn(zoltanNumObj<Adapter>, (void *) &(*adapter));
116  zz->Set_Obj_List_Fn(zoltanObjList<Adapter>, (void *) &(*adapter));
117 
118  const part_t *myparts;
119  adapter->getPartsView(myparts);
120  if (myparts != NULL)
121  zz->Set_Part_Multi_Fn(zoltanParts<Adapter>, (void *) &(*adapter));
122 
123  char tmp[4];
124  sprintf(tmp, "%d", TPL_Traits<ZOLTAN_ID_PTR, gno_t>::NUM_ID);
125  zz->Set_Param("NUM_GID_ENTRIES", tmp);
126  sprintf(tmp, "%d", TPL_Traits<ZOLTAN_ID_PTR, lno_t>::NUM_ID);
127  zz->Set_Param("NUM_LID_ENTRIES", tmp);
128  }
129 
130  template <typename AdapterWithCoords>
131  void setCallbacksGeom(const AdapterWithCoords *ia)
132  {
133  // Coordinates may be provided by the MeshAdapter or VectorAdapter.
134  // VectorAdapter may be provided directly by user or indirectly through
135  // GraphAdapter or MatrixAdapter. So separate template type is needed.
136  zz->Set_Num_Geom_Fn(zoltanNumGeom<AdapterWithCoords>, (void *) ia);
137  zz->Set_Geom_Multi_Fn(zoltanGeom<AdapterWithCoords>, (void *) ia);
138  }
139 
140  void setCallbacksGraph(
141  const RCP<const GraphAdapter<user_t,userCoord_t> > &/* adp */)
142  {
143  // std::cout << "NotReadyForGraphCallbacksYet" << std::endl;
144  // TODO
145  }
146 
147  void setCallbacksGraph(
148  const RCP<const MatrixAdapter<user_t,userCoord_t> > &adp)
149  {
150  // std::cout << "NotReadyForGraphCallbacksYet" << std::endl;
151  // TODO
152  }
153 
154  void setCallbacksGraph(
155  const RCP<const MeshAdapter<user_t> > &adp)
156  {
157  // std::cout << "NotReadyForGraphCallbacksYet" << std::endl;
158  // TODO
159  }
160 
161  void setCallbacksHypergraph(
162  const RCP<const MatrixAdapter<user_t,userCoord_t> > &adp)
163  {
164  // TODO: If add parameter list to this function, can register
165  // TODO: different callbacks depending on the hypergraph model to use
166 
167  zz->Set_HG_Size_CS_Fn(zoltanHGSizeCS_withMatrixAdapter<Adapter>,
168  (void *) &(*adp));
169  zz->Set_HG_CS_Fn(zoltanHGCS_withMatrixAdapter<Adapter>,
170  (void *) &(*adp));
171 
172  // zz->Set_HG_Size_Edge_Wts_Fn(zoltanHGSizeEdgeWts_withMatrixAdapter<Adapter>,
173  // (void *) &(*adapter));
174  // zz->Set_HG_Edge_Wts_Fn(zoltanHGSizeEdgeWts_withMatrixAdapter<Adapter>,
175  // (void *) &(*adapter));
176  }
177 
178  void setCallbacksHypergraph(
179  const RCP<const GraphAdapter<user_t,userCoord_t> > &adp)
180  {
181  zz->Set_HG_Size_CS_Fn(zoltanHGSizeCS_withGraphAdapter<Adapter>,
182  (void *) &(*adp));
183  zz->Set_HG_CS_Fn(zoltanHGCS_withGraphAdapter<Adapter>,
184  (void *) &(*adp));
185 
186  if (adp->getNumWeightsPerEdge() != 0) {
187  if (adp->getNumWeightsPerEdge() > 1) {
188  std::cout << "Zoltan2 warning: getNumWeightsPerEdge() returned "
189  << adp->getNumWeightsPerEdge() << " but PHG supports only "
190  << " one weight per edge; only first weight will be used."
191  << std::endl;
192  }
193  zz->Set_HG_Size_Edge_Wts_Fn(zoltanHGSizeEdgeWts_withGraphAdapter<Adapter>,
194  (void *) &(*adapter));
195  zz->Set_HG_Edge_Wts_Fn(zoltanHGEdgeWts_withGraphAdapter<Adapter>,
196  (void *) &(*adapter));
197  }
198  }
199 
200  void setCallbacksHypergraph(const RCP<const MeshAdapter<user_t> > &adp)
201  {
202 
203  const Teuchos::ParameterList &pl = env->getParameters();
204 
205  const Teuchos::ParameterEntry *pe = pl.getEntryPtr("hypergraph_model_type");
206  std::string model_type("traditional");
207  if (pe){
208  model_type = pe->getValue<std::string>(&model_type);
209  }
210 
211  if (model_type=="ghosting" ||
212  !adp->areEntityIDsUnique(adp->getPrimaryEntityType())) {
213  Zoltan2::modelFlag_t flags;
215  problemComm, flags,
217  model = rcp(static_cast<const Model<Adapter>* >(mdl),true);
218 
219  zz->Set_Num_Obj_Fn(zoltanHGNumObj_withModel<Adapter>, (void *) &(*mdl));
220  zz->Set_Obj_List_Fn(zoltanHGObjList_withModel<Adapter>, (void *) &(*mdl));
221 
222  zz->Set_HG_Size_CS_Fn(zoltanHGSizeCS_withModel<Adapter>, (void *) &(*mdl));
223  zz->Set_HG_CS_Fn(zoltanHGCS_withModel<Adapter>, (void *) &(*mdl));
224  }
225  else {
226  //If entities are unique we dont need the extra cost of the model
227  zz->Set_HG_Size_CS_Fn(zoltanHGSizeCS_withMeshAdapter<Adapter>,
228  (void *) &(*adp));
229  zz->Set_HG_CS_Fn(zoltanHGCS_withMeshAdapter<Adapter>,
230  (void *) &(*adp));
231  }
232  // zz->Set_HG_Size_Edge_Wts_Fn(zoltanHGSizeEdgeWts_withMeshAdapter<Adapter>,
233  // (void *) &(*adp));
234  // zz->Set_HG_Edge_Wts_Fn(zoltanHGSizeEdgeWts_withMeshAdapter<Adapter>,
235  // (void *) &(*adp));
236  }
237 
239  virtual bool isPartitioningTreeBinary() const
240  {
241  return true;
242  }
243 
245  void rcb_recursive_partitionTree_calculations(
246  part_t arrayIndex,
247  part_t numParts,
248  std::vector<part_t> & splitRangeBeg,
249  std::vector<part_t> & splitRangeEnd) const
250  {
251  // Note the purpose of the recursive method is make sure the children of a
252  // node have updated their values for splitRangeBeg and splitRangeEnd
253  // Then we can set our own values simply based on the union
254  // first load the rcb data for the node
255  int parent = -1;
256  int left_leaf = -1;
257  int right_leaf = -1;
258  int err = Zoltan_RCB_Partition_Tree(zz->Get_C_Handle(),
259  arrayIndex - numParts + 1, // rcb starts as 1 but does not include terminals
260  &parent, &left_leaf, &right_leaf);
261  if(err != 0) {
262  throw std::logic_error( "Zoltan_RCB_Partition_Tree returned in error." );
263  }
264  // check that children both have their ranges set and if not, do those
265  // range first so we can build them to make our range
266  if(left_leaf > 0) { // neg is terminal and always already built
267  rcb_recursive_partitionTree_calculations(left_leaf+numParts-1, numParts,
268  splitRangeBeg, splitRangeEnd);
269  }
270  if(right_leaf > 0) { // neg is terminal and always already built
271  rcb_recursive_partitionTree_calculations(right_leaf+numParts-1, numParts,
272  splitRangeBeg, splitRangeEnd);
273  }
274  // now we can build our ranges from the children
275  // note this exploits the rcb conventions for right and left so we know
276  // that left_leaf will be our smaller indices
277  int leftIndex = (left_leaf > 0) ? (left_leaf-1+numParts) : (-left_leaf);
278  int rightIndex = (right_leaf > 0) ? (right_leaf-1+numParts) : (-right_leaf);
279  splitRangeBeg[arrayIndex] = splitRangeBeg[leftIndex];
280  splitRangeEnd[arrayIndex] = splitRangeEnd[rightIndex];
281  // for debugging sanity check verify left_leaf is a set of indices which
282  // goes continuously into the right_leaf
283  if(splitRangeBeg[rightIndex] != splitRangeEnd[leftIndex]) { // end is non-inclusive
284  throw std::logic_error( "RCB expected left_leaf indices and right leaf"
285  " indices to be continuous but it was not so." );
286  }
287  }
288 
290  void rcb_getPartitionTree(part_t numParts,
291  part_t & numTreeVerts,
292  std::vector<part_t> & permPartNums,
293  std::vector<part_t> & splitRangeBeg,
294  std::vector<part_t> & splitRangeEnd,
295  std::vector<part_t> & treeVertParents) const
296  {
297  // CALCULATE: numTreeVerts
298  // For rcb a tree node always takes 2 nodes and turns them into 1 node
299  // That means it takes numParts - 1 nodes to reduce a tree of numParts to
300  // a single root node - but we do 2 * numParts - 1 because we are currently
301  // treating all of the 'trivial' terminals as tree nodes - something we
302  // discussed we may change later
303  part_t numTreeNodes = 2 * numParts - 1;
304  numTreeVerts = numTreeNodes - 1; // by design convention root not included
305  // CALCULATE: permPartNums
306  permPartNums.resize(numParts);
307  for(part_t n = 0; n < numParts; ++n) {
308  permPartNums[n] = n; // for rcb we can assume this is trivial and in order
309  }
310  // CALCULATE: treeVertParents
311  treeVertParents.resize(numTreeNodes); // allocate space for numTreeNodes array
312  // scan all the non terminal nodes and set all children to have us as parent
313  // that will set all parents except for the root node which we will set to -1
314  // track the children of the root and final node for swapping later. Couple
315  // ways to do this - all seem a bit awkward but at least this is efficient.
316  part_t rootNode = 0; // track index of the root node for swapping
317  // a bit awkward but efficient - save the children of root and final node
318  // for swap at end to satisfy convention that root is highest index node
319  part_t saveRootNodeChildrenA = -1;
320  part_t saveRootNodeChildrenB = -1;
321  part_t saveFinalNodeChildrenA = -1;
322  part_t saveFinalNodeChildrenB = -1;
323  for(part_t n = numParts; n < numTreeNodes; ++n) { // scan and set all parents
324  int parent = -1;
325  int left_leaf = -1;
326  int right_leaf = -1;
327  int err = Zoltan_RCB_Partition_Tree(zz->Get_C_Handle(),
328  static_cast<int>(n - numParts) + 1, // rcb starts as 1 but does not include terminals
329  &parent, &left_leaf, &right_leaf);
330  if(err != 0) {
331  throw std::logic_error("Zoltan_RCB_Partition_Tree returned in error.");
332  }
333  part_t leftIndex = (left_leaf > 0) ? (left_leaf-1+numParts) : (-left_leaf);
334  part_t rightIndex = (right_leaf > 0) ? (right_leaf-1+numParts) : (-right_leaf);
335  treeVertParents[leftIndex] = n;
336  treeVertParents[rightIndex] = n;
337  // save root node for final swap
338  if(parent == 1 || parent == -1) { // is it the root?
339  rootNode = n; // remember I am the root
340  saveRootNodeChildrenA = leftIndex;
341  saveRootNodeChildrenB = rightIndex;
342  }
343  if(n == numTreeNodes-1) {
344  saveFinalNodeChildrenA = leftIndex;
345  saveFinalNodeChildrenB = rightIndex;
346  }
347  }
348  treeVertParents[rootNode] = -1; // convention parent is root -1
349  // splitRangeBeg and splitRangeEnd
350  splitRangeBeg.resize(numTreeNodes);
351  splitRangeEnd.resize(numTreeNodes);
352  // for terminal nodes this is trivial
353  for(part_t n = 0; n < numParts; ++n) {
354  splitRangeBeg[n] = n;
355  splitRangeEnd[n] = n + 1;
356  }
357  if(numParts > 1) { // not relevant for 1 part
358  // build the splitRangeBeg and splitRangeEnd recursively which forces the
359  // children of each node to be calculated before the parent so parent can
360  // just take the union of the two children
361  rcb_recursive_partitionTree_calculations(rootNode, numParts, splitRangeBeg,
362  splitRangeEnd);
363  // now as a final step handle the swap to root is the highest index node
364  // swap the parent of the two nodes
365  std::swap(treeVertParents[rootNode], treeVertParents[numTreeNodes-1]);
366  // get the children of the swapped nodes to have updated parents
367  treeVertParents[saveFinalNodeChildrenA] = rootNode;
368  treeVertParents[saveFinalNodeChildrenB] = rootNode;
369  // handle case where final node is child of the root
370  if(saveRootNodeChildrenA == numTreeNodes - 1) {
371  saveRootNodeChildrenA = rootNode;
372  }
373  if(saveRootNodeChildrenB == numTreeNodes - 1) {
374  saveRootNodeChildrenB = rootNode;
375  }
376  treeVertParents[saveRootNodeChildrenA] = numTreeNodes - 1;
377  treeVertParents[saveRootNodeChildrenB] = numTreeNodes - 1;
378  // update the beg and end indices - simply swap them
379  std::swap(splitRangeBeg[rootNode], splitRangeBeg[numTreeNodes-1]);
380  std::swap(splitRangeEnd[rootNode], splitRangeEnd[numTreeNodes-1]);
381  }
382  }
383 
385  void phg_getPartitionTree(part_t numParts,
386  part_t & numTreeVerts,
387  std::vector<part_t> & permPartNums,
388  std::vector<part_t> & splitRangeBeg,
389  std::vector<part_t> & splitRangeEnd,
390  std::vector<part_t> & treeVertParents) const
391  {
392  // First thing is to get the length of the tree from zoltan.
393  // The tree is a list of pairs (lo,hi) for each node.
394  // Here tree_array_size is the number of pairs.
395  // In phg indexing the first pairt (i=0) is empty garbage.
396  // The second pair (index 1) will be the root.
397  // Some nodes will be empty nodes, determined by hi = -1.
398  int tree_array_size = -1; // will be set by Zoltan_PHG_Partition_Tree_Size
399  int err = Zoltan_PHG_Partition_Tree_Size(
400  zz->Get_C_Handle(), &tree_array_size);
401  if(err != 0) {
402  throw std::logic_error("Zoltan_PHG_Partition_Tree_Size returned error.");
403  }
404  // Determine the number of valid nodes (PHG will have empty slots)
405  // We scan the list of pairs and count each node which does not have hi = -1
406  // During the loop we will also construct mapIndex which maps initial n
407  // to final n due to some conversions we apply to meet the design specs.
408  // The conversions implemented by mapIndex are:
409  // Move all terminals to the beginning (terminals have hi = lo)
410  // Resort the terminals in order (simply map to index lo works)
411  // Move non-terminals after the terminals (they start at index numParts)
412  // Map the first pair (root) to the be last to meet the design spec
413  part_t numTreeNodes = 0;
414  std::vector<part_t> mapIndex(tree_array_size); // maps n to final index
415  part_t trackNonTerminalIndex = numParts; // starts after terminals
416  for(part_t n = 0; n < static_cast<part_t>(tree_array_size); ++n) {
417  part_t phgIndex = n + 1; // phg indexing starts at 1
418  int lo_index = -1;
419  int hi_index = -1;
420  err = Zoltan_PHG_Partition_Tree(
421  zz->Get_C_Handle(), phgIndex, &lo_index, &hi_index);
422  if(hi_index != -1) { // hi -1 means it's an unused node
423  ++numTreeNodes; // increase the total count because this is a real node
424  if(n != 0) { // the root is mapped last but we don't know the length yet
425  mapIndex[n] = (lo_index == hi_index) ? // is it a terminal?
426  lo_index : // terminals map in sequence - lo_index is correct
427  (trackNonTerminalIndex++); // set then bump trackNonTerminalIndex +1
428  }
429  }
430  }
431  // now complete the map by setting root to the length-1 for the design specs
432  mapIndex[0] = numTreeNodes - 1;
433  // CALCULATE: numTreeVerts
434  numTreeVerts = numTreeNodes - 1; // this is the design - root not included
435  // CALCULATE: permPartNums
436  permPartNums.resize(numParts);
437  for(part_t n = 0; n < numParts; ++n) {
438  permPartNums[n] = n; // for phg we can assume this is trivial and in order
439  }
440  // CALCULATE: treeVertParents, splitRangeBeg, splitRangeEnd
441  // we will determine all of these in this second loop using mapIndex
442  // First set the arrays to have the proper length
443  treeVertParents.resize(numTreeNodes);
444  splitRangeBeg.resize(numTreeNodes);
445  splitRangeEnd.resize(numTreeNodes);
446  // Now loop a second time
447  for(part_t n = 0; n < tree_array_size; ++n) {
448  part_t phgIndex = n + 1; // phg indexing starts at 1
449  int lo_index = -1; // zoltan Zoltan_PHG_Partition_Tree will set this
450  int hi_index = -1; // zoltan Zoltan_PHG_Partition_Tree will set this
451  err = Zoltan_PHG_Partition_Tree( // access zoltan phg tree data
452  zz->Get_C_Handle(), phgIndex, &lo_index, &hi_index);
453  if(err != 0) {
454  throw std::logic_error("Zoltan_PHG_Partition_Tree returned in error.");
455  }
456  if(hi_index != -1) { // hi -1 means it's an unused node (a gap)
457  // get final index using mapIndex - so convert from phg to design plan
458  part_t finalIndex = mapIndex[n]; // get the index for the final output
459  // now determine the parent
460  // In the original phg indexing, the parent can be directly calculated
461  // from the pair index using the following rules:
462  // if phgIndex even, then parent is phgIndex/2
463  // here we determine even by ((phgIndex%2) == 0)
464  // if phgIndex odd, then parent is (phgIndex-1)/2
465  // but after getting parentPHGIndex we convert back to this array
466  // indexing by subtracting 1
467  part_t parentPHGIndex =
468  ((phgIndex%2) == 0) ? (phgIndex/2) : ((phgIndex-1)/2);
469  // map the parent phg index to the final parent index
470  // however, for the special case of the root (n=1), set the parent to -1
471  treeVertParents[finalIndex] = (n==0) ? -1 : mapIndex[parentPHGIndex-1];
472  // set begin (inclusive) and end (non-inclusive), so end is hi+1
473  splitRangeBeg[finalIndex] = static_cast<part_t>(lo_index);
474  splitRangeEnd[finalIndex] = static_cast<part_t>(hi_index+1);
475  }
476  }
477  }
478 
480  void getPartitionTree(part_t numParts,
481  part_t & numTreeVerts,
482  std::vector<part_t> & permPartNums,
483  std::vector<part_t> & splitRangeBeg,
484  std::vector<part_t> & splitRangeEnd,
485  std::vector<part_t> & treeVertParents) const
486  {
487  // first check that our parameters requested we keep the tree
488  const Teuchos::ParameterList &pl = env->getParameters();
489  bool keep_partition_tree = false;
490  const Teuchos::ParameterEntry * pe = pl.getEntryPtr("keep_partition_tree");
491  if(pe) {
492  keep_partition_tree = pe->getValue(&keep_partition_tree);
493  if(!keep_partition_tree) {
494  throw std::logic_error(
495  "Requested tree when param keep_partition_tree is off.");
496  }
497  }
498 
499  // now call the appropriate method based on LB_METHOD in the zoltan
500  // parameters list.
501  const Teuchos::ParameterList & zoltan_pl = pl.sublist("zoltan_parameters");
502  std::string lb_method;
503  pe = zoltan_pl.getEntryPtr("LB_METHOD");
504  if(pe) {
505  lb_method = pe->getValue(&lb_method);
506  }
507  if(lb_method == "phg") {
508  phg_getPartitionTree(numParts, numTreeVerts, permPartNums,
509  splitRangeBeg, splitRangeEnd, treeVertParents);
510  }
511  else if(lb_method == "rcb") {
512  rcb_getPartitionTree(numParts, numTreeVerts, permPartNums,
513  splitRangeBeg, splitRangeEnd, treeVertParents);
514  }
515  else {
516  throw std::logic_error("Did not recognize LB_METHOD: " + lb_method);
517  }
518  }
519 
520 public:
521 
527  AlgZoltan(const RCP<const Environment> &env__,
528  const RCP<const Comm<int> > &problemComm__,
529  const RCP<const IdentifierAdapter<user_t> > &adapter__):
530  env(env__), problemComm(problemComm__), adapter(adapter__)
531  {
532  setMPIComm(problemComm__);
533  zoltanInit();
534  zz = rcp(new Zoltan(mpicomm));
535  setCallbacksIDs();
536  }
537 
538  AlgZoltan(const RCP<const Environment> &env__,
539  const RCP<const Comm<int> > &problemComm__,
540  const RCP<const VectorAdapter<user_t> > &adapter__) :
541  env(env__), problemComm(problemComm__), adapter(adapter__)
542  {
543  setMPIComm(problemComm__);
544  zoltanInit();
545  zz = rcp(new Zoltan(mpicomm));
546  setCallbacksIDs();
547  setCallbacksGeom(&(*adapter));
548  }
549 
550  AlgZoltan(const RCP<const Environment> &env__,
551  const RCP<const Comm<int> > &problemComm__,
552  const RCP<const GraphAdapter<user_t,userCoord_t> > &adapter__) :
553  env(env__), problemComm(problemComm__), adapter(adapter__)
554  {
555  setMPIComm(problemComm__);
556  zoltanInit();
557  zz = rcp(new Zoltan(mpicomm));
558  setCallbacksIDs();
559  setCallbacksGraph(adapter);
560  setCallbacksHypergraph(adapter);
561  if (adapter->coordinatesAvailable()) {
562  setCallbacksGeom(adapter->getCoordinateInput());
563  }
564  }
565 
566  AlgZoltan(const RCP<const Environment> &env__,
567  const RCP<const Comm<int> > &problemComm__,
568  const RCP<const MatrixAdapter<user_t,userCoord_t> > &adapter__) :
569  env(env__), problemComm(problemComm__), adapter(adapter__)
570  {
571  setMPIComm(problemComm__);
572  zoltanInit();
573  zz = rcp(new Zoltan(mpicomm));
574  setCallbacksIDs();
575  setCallbacksGraph(adapter);
576  setCallbacksHypergraph(adapter);
577  if (adapter->coordinatesAvailable()) {
578  setCallbacksGeom(adapter->getCoordinateInput());
579  }
580  }
581 
582  AlgZoltan(const RCP<const Environment> &env__,
583  const RCP<const Comm<int> > &problemComm__,
584  const RCP<const MeshAdapter<user_t> > &adapter__) :
585  env(env__), problemComm(problemComm__), adapter(adapter__)
586  {
587  setMPIComm(problemComm__);
588  zoltanInit();
589  zz = rcp(new Zoltan(mpicomm));
590  setCallbacksIDs();
591  setCallbacksGraph(adapter);
592  //TODO:: check parameter list to see if hypergraph is needed. We dont want to build the model
593  // if we don't have to and we shouldn't as it can take a decent amount of time if the
594  // primary entity is copied
595  setCallbacksHypergraph(adapter);
596  setCallbacksGeom(&(*adapter));
597  }
598 
599  void partition(const RCP<PartitioningSolution<Adapter> > &solution);
600  // void color(const RCP<ColoringSolution<Adapter> > &solution);
601 };
602 
604 template <typename Adapter>
606  const RCP<PartitioningSolution<Adapter> > &solution
607 )
608 {
609  HELLO;
610  char paramstr[128];
611 
612  size_t numGlobalParts = solution->getTargetGlobalNumberOfParts();
613 
614  sprintf(paramstr, "%lu", numGlobalParts);
615  zz->Set_Param("NUM_GLOBAL_PARTS", paramstr);
616 
617  int wdim = adapter->getNumWeightsPerID();
618  sprintf(paramstr, "%d", wdim);
619  zz->Set_Param("OBJ_WEIGHT_DIM", paramstr);
620 
621  const Teuchos::ParameterList &pl = env->getParameters();
622 
623  double tolerance;
624  const Teuchos::ParameterEntry *pe = pl.getEntryPtr("imbalance_tolerance");
625  if (pe){
626  char str[30];
627  tolerance = pe->getValue<double>(&tolerance);
628  sprintf(str, "%f", tolerance);
629  zz->Set_Param("IMBALANCE_TOL", str);
630  }
631 
632  pe = pl.getEntryPtr("partitioning_approach");
633  if (pe){
634  std::string approach;
635  approach = pe->getValue<std::string>(&approach);
636  if (approach == "partition")
637  zz->Set_Param("LB_APPROACH", "PARTITION");
638  else
639  zz->Set_Param("LB_APPROACH", "REPARTITION");
640  }
641 
642  pe = pl.getEntryPtr("partitioning_objective");
643  if (pe){
644  std::string strChoice = pe->getValue<std::string>(&strChoice);
645  if (strChoice == std::string("multicriteria_minimize_total_weight"))
646  zz->Set_Param("RCB_MULTICRITERIA_NORM", "1");
647  else if (strChoice == std::string("multicriteria_balance_total_maximum"))
648  zz->Set_Param("RCB_MULTICRITERIA_NORM", "2");
649  else if (strChoice == std::string("multicriteria_minimize_maximum_weight"))
650  zz->Set_Param("RCB_MULTICRITERIA_NORM", "3");
651  }
652 
653  // perhaps make this a bool stored in the AlgZoltan if we want to follow
654  // the pattern of multijagged mj_keep_part_boxes for example. However we can
655  // collect the error straight from Zoltan if we attempt to access the tree
656  // when we never stored it so that may not be necessary
657  bool keep_partition_tree = false;
658  pe = pl.getEntryPtr("keep_partition_tree");
659  if (pe) {
660  keep_partition_tree = pe->getValue(&keep_partition_tree);
661  if (keep_partition_tree) {
662  // need to resolve the organization of this
663  // when do we want to use the zoltan parameters directly versus
664  // using the zoltan2 parameters like this
665  zz->Set_Param("KEEP_CUTS", "1"); // rcb zoltan setting
666  zz->Set_Param("PHG_KEEP_TREE", "1"); // phg zoltan setting
667  }
668  }
669 
670  pe = pl.getEntryPtr("rectilinear");
671  if (pe) {
672  bool val = pe->getValue(&val);
673  if (val)
674  zz->Set_Param("RCB_RECTILINEAR_BLOCKS", "1");
675  }
676 
677  // Look for zoltan_parameters sublist; pass all zoltan parameters to Zoltan
678  try {
679  const Teuchos::ParameterList &zpl = pl.sublist("zoltan_parameters");
680  for (ParameterList::ConstIterator iter = zpl.begin();
681  iter != zpl.end(); iter++) {
682  const std::string &zname = pl.name(iter);
683  // Convert the value to a string to pass to Zoltan
684  std::string zval = pl.entry(iter).getValue(&zval);
685  zz->Set_Param(zname.c_str(), zval.c_str());
686  }
687  }
688  catch (std::exception &e) {
689  // No zoltan_parameters sublist found; no Zoltan parameters to register
690  }
691 
692  // Get target part sizes
693  int pdim = (wdim > 1 ? wdim : 1);
694  for (int d = 0; d < pdim; d++) {
695  if (!solution->criteriaHasUniformPartSizes(d)) {
696  float *partsizes = new float[numGlobalParts];
697  int *partidx = new int[numGlobalParts];
698  int *wgtidx = new int[numGlobalParts];
699  for (size_t i=0; i<numGlobalParts; i++) partidx[i] = i;
700  for (size_t i=0; i<numGlobalParts; i++) wgtidx[i] = d;
701  for (size_t i=0; i<numGlobalParts; i++)
702  partsizes[i] = solution->getCriteriaPartSize(0, i);
703  zz->LB_Set_Part_Sizes(1, numGlobalParts, partidx, wgtidx, partsizes);
704  delete [] partsizes;
705  delete [] partidx;
706  delete [] wgtidx;
707  }
708  }
709 
710  // Make the call to LB_Partition
711  int changed = 0;
714 
715  int nDummy = -1; // Dummy vars to satisfy arglist
716  ZOLTAN_ID_PTR dGids = NULL, dLids = NULL;
717  int *dProcs = NULL, *dParts = NULL;
718  int nObj = -1; // Output vars with new part info
719  ZOLTAN_ID_PTR oGids = NULL, oLids = NULL;
720  int *oProcs = NULL, *oParts = NULL;
721 
722  zz->Set_Param("RETURN_LISTS", "PARTS"); // required format for Zoltan2;
723  // results in last five arguments
724 
725  int ierr = zz->LB_Partition(changed, nGidEnt, nLidEnt,
726  nDummy, dGids, dLids, dProcs, dParts,
727  nObj, oGids, oLids, oProcs, oParts);
728 
729  env->globalInputAssertion(__FILE__, __LINE__, "Zoltan LB_Partition",
730  (ierr==ZOLTAN_OK || ierr==ZOLTAN_WARN), BASIC_ASSERTION, problemComm);
731 
732  int numObjects=nObj;
733  // The number of objects may be larger than zoltan knows due to copies that
734  // were removed by the hypergraph model
735  if (model!=RCP<const Model<Adapter> >() &&
736  dynamic_cast<const HyperGraphModel<Adapter>* >(&(*model)) &&
737  !(dynamic_cast<const HyperGraphModel<Adapter>* >(&(*model))->areVertexIDsUnique())) {
738  numObjects=model->getLocalNumObjects();
739  }
740 
741  // Load answer into the solution.
742  ArrayRCP<part_t> partList(new part_t[numObjects], 0, numObjects, true);
743  for (int i = 0; i < nObj; i++) {
744  lno_t tmp;
745  TPL_Traits<lno_t, ZOLTAN_ID_PTR>::ASSIGN(tmp, &(oLids[i*nLidEnt]));
746  partList[tmp] = oParts[i];
747  }
748 
749  if (model!=RCP<const Model<Adapter> >() &&
750  dynamic_cast<const HyperGraphModel<Adapter>* >(&(*model)) &&
751  !(dynamic_cast<const HyperGraphModel<Adapter>* >(&(*model))->areVertexIDsUnique())) {
752  // Setup the part ids for copied entities removed by ownership in
753  // hypergraph model.
754  const HyperGraphModel<Adapter>* mdl =
755  static_cast<const HyperGraphModel<Adapter>* >(&(*model));
756 
757  typedef typename HyperGraphModel<Adapter>::map_t map_t;
758  Teuchos::RCP<const map_t> mapWithCopies;
759  Teuchos::RCP<const map_t> oneToOneMap;
760  mdl->getVertexMaps(mapWithCopies,oneToOneMap);
761 
762  typedef Tpetra::Vector<scalar_t, lno_t, gno_t> vector_t;
763  vector_t vecWithCopies(mapWithCopies);
764  vector_t oneToOneVec(oneToOneMap);
765 
766  // Set values in oneToOneVec: each entry == rank
767  assert(nObj == lno_t(oneToOneMap->getNodeNumElements()));
768  for (lno_t i = 0; i < nObj; i++)
769  oneToOneVec.replaceLocalValue(i, oParts[i]);
770 
771  // Now import oneToOneVec's values back to vecWithCopies
772  Teuchos::RCP<const Tpetra::Import<lno_t, gno_t> > importer =
773  Tpetra::createImport<lno_t, gno_t>(oneToOneMap, mapWithCopies);
774  vecWithCopies.doImport(oneToOneVec, *importer, Tpetra::REPLACE);
775 
776  // Should see copied vector values when print VEC WITH COPIES
777  lno_t nlocal = lno_t(mapWithCopies->getNodeNumElements());
778  for (lno_t i = 0; i < nlocal; i++)
779  partList[i] = vecWithCopies.getData()[i];
780  }
781 
782  solution->setParts(partList);
783 
784  // Clean up
785  zz->LB_Free_Part(&oGids, &oLids, &oProcs, &oParts);
786 }
787 
788 } // namespace Zoltan2
789 
790 #endif
AlgZoltan(const RCP< const Environment > &env__, const RCP< const Comm< int > > &problemComm__, const RCP< const VectorAdapter< user_t > > &adapter__)
#define HELLO
IdentifierAdapter defines the interface for identifiers.
fast typical checks for valid arguments
MatrixAdapter defines the adapter interface for matrices.
Defines the Model interface.
AlgZoltan(const RCP< const Environment > &env__, const RCP< const Comm< int > > &problemComm__, const RCP< const MatrixAdapter< user_t, userCoord_t > > &adapter__)
GraphAdapter defines the interface for graph-based user data.
MeshAdapter defines the interface for mesh input.
std::bitset< NUM_MODEL_FLAGS > modelFlag_t
Defines the PartitioningSolution class.
SparseMatrixAdapter_t::part_t part_t
void partition(const RCP< PartitioningSolution< Adapter > > &solution)
Partitioning method.
callback functions for the Zoltan package (templated on Adapter)
A PartitioningSolution is a solution to a partitioning problem.
VectorAdapter defines the interface for vector input.
AlgZoltan(const RCP< const Environment > &env__, const RCP< const Comm< int > > &problemComm__, const RCP< const MeshAdapter< user_t > > &adapter__)
Algorithm defines the base class for all algorithms.
static void ASSIGN(first_t &a, second_t b)
Traits class to handle conversions between gno_t/lno_t and TPL data types (e.g., ParMETIS&#39;s idx_t...
Gathering definitions used in software development.
The base class for all model classes.
AlgZoltan(const RCP< const Environment > &env__, const RCP< const Comm< int > > &problemComm__, const RCP< const GraphAdapter< user_t, userCoord_t > > &adapter__)
A gathering of useful namespace methods.
HyperGraphModel defines the interface required for hyper graph models.
void getVertexMaps(Teuchos::RCP< const map_t > &copiesMap, Teuchos::RCP< const map_t > &onetooneMap) const
Sets pointers to the vertex map with copies and the vertex map without copies Note: the pointers will...
AlgZoltan(const RCP< const Environment > &env__, const RCP< const Comm< int > > &problemComm__, const RCP< const IdentifierAdapter< user_t > > &adapter__)
Zoltan2::BasicUserTypes< zscalar_t, zlno_t, zgno_t > user_t
Definition: Metric.cpp:74