Zoltan2
Zoltan2_AlgPuLP.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_ALGPULP_HPP_
46 #define _ZOLTAN2_ALGPULP_HPP_
47 
48 #include <Zoltan2_GraphModel.hpp>
49 #include <Zoltan2_Algorithm.hpp>
51 #include <Zoltan2_Util.hpp>
52 #include <Zoltan2_TPLTraits.hpp>
53 
57 
59 #ifndef HAVE_ZOLTAN2_PULP
60 
61 namespace Zoltan2 {
62 // Error handling for when PuLP is requested
63 // but Zoltan2 not built with PuLP.
64 
65 template <typename Adapter>
66 class AlgPuLP : public Algorithm<Adapter>
67 {
68 public:
70  AlgPuLP(const RCP<const Environment> &env,
71  const RCP<const Comm<int> > &problemComm,
72  const RCP<const base_adapter_t> &adapter
73  )
74  {
75  throw std::runtime_error(
76  "BUILD ERROR: PuLP requested but not compiled into Zoltan2.\n"
77  "Please set CMake flag Zoltan2_ENABLE_PuLP:BOOL=ON.");
78  }
79 
82  static void getValidParameters(ParameterList & pl)
83  {
84  pl.set("pulp_vert_imbalance", 1.1, "vertex imbalance tolerance, ratio of "
85  "maximum load over average load",
87 
88  pl.set("pulp_edge_imbalance", 1.1, "edge imbalance tolerance, ratio of "
89  "maximum load over average load",
91 
92  // bool parameter
93  pl.set("pulp_lp_init", false, "perform label propagation-based "
94  "initialization", Environment::getBoolValidator() );
95 
96  // bool parameter
97  pl.set("pulp_minimize_maxcut", false, "perform per-part max cut "
98  "minimization", Environment::getBoolValidator() );
99 
100  // bool parameter
101  pl.set("pulp_verbose", false, "verbose output",
103 
104  // bool parameter
105  pl.set("pulp_do_repart", false, "perform repartitioning",
107 
108  pl.set("pulp_seed", 0, "set pulp seed", Environment::getAnyIntValidator());
109  }
110 };
111 
112 }
113 #endif
114 
117 
118 #ifdef HAVE_ZOLTAN2_PULP
119 
120 namespace Zoltan2 {
121 
122 extern "C" {
123 // TODO: XtraPuLP
124 #ifndef HAVE_ZOLTAN2_MPI
125 #include "pulp.h"
126 #else
127 #include "xtrapulp.h"
128 #endif
129 }
130 
131 
132 
133 
134 template <typename Adapter>
135 class AlgPuLP : public Algorithm<Adapter>
136 {
137 public:
138  typedef typename Adapter::base_adapter_t base_adapter_t;
139  typedef typename Adapter::lno_t lno_t;
140  typedef typename Adapter::gno_t gno_t;
141  typedef typename Adapter::scalar_t scalar_t;
142  typedef typename Adapter::part_t part_t;
143  typedef typename Adapter::user_t user_t;
144  typedef typename Adapter::userCoord_t userCoord_t;
145 
156  AlgPuLP(const RCP<const Environment> &env__,
157  const RCP<const Comm<int> > &problemComm__,
158  const RCP<const IdentifierAdapter<user_t> > &adapter__) :
159  env(env__), problemComm(problemComm__), adapter(adapter__)
160  {
161  std::string errStr = "cannot build GraphModel from IdentifierAdapter, ";
162  errStr += "PuLP requires Graph, Matrix, or Mesh Adapter";
163  throw std::runtime_error(errStr);
164  }
165 
166  AlgPuLP(const RCP<const Environment> &env__,
167  const RCP<const Comm<int> > &problemComm__,
168  const RCP<const VectorAdapter<user_t> > &adapter__) :
169  env(env__), problemComm(problemComm__), adapter(adapter__)
170  {
171  std::string errStr = "cannot build GraphModel from VectorAdapter, ";
172  errStr += "PuLP requires Graph, Matrix, or Mesh Adapter";
173  throw std::runtime_error(errStr);
174  }
175 
176  AlgPuLP(const RCP<const Environment> &env__,
177  const RCP<const Comm<int> > &problemComm__,
178  const RCP<const GraphAdapter<user_t,userCoord_t> > &adapter__) :
179  env(env__), problemComm(problemComm__), adapter(adapter__)
180  {
181  modelFlag_t flags;
182  flags.reset();
183 
184  buildModel(flags);
185  }
186 
187  AlgPuLP(const RCP<const Environment> &env__,
188  const RCP<const Comm<int> > &problemComm__,
189  const RCP<const MatrixAdapter<user_t,userCoord_t> > &adapter__) :
190  env(env__), problemComm(problemComm__), adapter(adapter__)
191  {
192  modelFlag_t flags;
193  flags.reset();
194 
195  const ParameterList &pl = env->getParameters();
196  const Teuchos::ParameterEntry *pe;
197 
198  std::string defString("default");
199  std::string objectOfInterest(defString);
200  pe = pl.getEntryPtr("objects_to_partition");
201  if (pe)
202  objectOfInterest = pe->getValue<std::string>(&objectOfInterest);
203 
204  if (objectOfInterest == defString ||
205  objectOfInterest == std::string("matrix_rows") )
206  flags.set(VERTICES_ARE_MATRIX_ROWS);
207  else if (objectOfInterest == std::string("matrix_columns"))
208  flags.set(VERTICES_ARE_MATRIX_COLUMNS);
209  else if (objectOfInterest == std::string("matrix_nonzeros"))
210  flags.set(VERTICES_ARE_MATRIX_NONZEROS);
211 
212  buildModel(flags);
213  }
214 
215  AlgPuLP(const RCP<const Environment> &env__,
216  const RCP<const Comm<int> > &problemComm__,
217  const RCP<const MeshAdapter<user_t> > &adapter__) :
218  env(env__), problemComm(problemComm__), adapter(adapter__)
219  {
220  modelFlag_t flags;
221  flags.reset();
222 
223  const ParameterList &pl = env->getParameters();
224  const Teuchos::ParameterEntry *pe;
225 
226  std::string defString("default");
227  std::string objectOfInterest(defString);
228  pe = pl.getEntryPtr("objects_to_partition");
229  if (pe)
230  objectOfInterest = pe->getValue<std::string>(&objectOfInterest);
231 
232  if (objectOfInterest == defString ||
233  objectOfInterest == std::string("mesh_nodes") )
234  flags.set(VERTICES_ARE_MESH_NODES);
235  else if (objectOfInterest == std::string("mesh_elements"))
236  flags.set(VERTICES_ARE_MESH_ELEMENTS);
237 
238  buildModel(flags);
239  }
240 
241  void partition(const RCP<PartitioningSolution<Adapter> > &solution);
242 
243 private:
244 
245  void buildModel(modelFlag_t &flags);
246 
247  void scale_weights(size_t n, StridedData<lno_t, scalar_t> &fwgts,
248  int *iwgts);
249 
250  const RCP<const Environment> env;
251  const RCP<const Comm<int> > problemComm;
252  const RCP<const base_adapter_t> adapter;
253  RCP<const GraphModel<base_adapter_t> > model;
254 };
255 
256 
258 template <typename Adapter>
259 void AlgPuLP<Adapter>::buildModel(modelFlag_t &flags)
260 {
261  const ParameterList &pl = env->getParameters();
262  const Teuchos::ParameterEntry *pe;
263 
264  std::string defString("default");
265  std::string symParameter(defString);
266  pe = pl.getEntryPtr("symmetrize_graph");
267  if (pe){
268  symParameter = pe->getValue<std::string>(&symParameter);
269  if (symParameter == std::string("transpose"))
270  flags.set(SYMMETRIZE_INPUT_TRANSPOSE);
271  else if (symParameter == std::string("bipartite"))
272  flags.set(SYMMETRIZE_INPUT_BIPARTITE); }
273 
274  bool sgParameter = false;
275  pe = pl.getEntryPtr("subset_graph");
276  if (pe)
277  sgParameter = pe->getValue(&sgParameter);
278  if (sgParameter)
279  flags.set(BUILD_SUBSET_GRAPH);
280 
281  flags.set(REMOVE_SELF_EDGES);
282  flags.set(GENERATE_CONSECUTIVE_IDS);
283 #ifndef HAVE_ZOLTAN2_MPI
284  flags.set(BUILD_LOCAL_GRAPH);
285 #endif
286  this->env->debug(DETAILED_STATUS, " building graph model");
287  this->model = rcp(new GraphModel<base_adapter_t>(this->adapter, this->env,
288  this->problemComm, flags));
289  this->env->debug(DETAILED_STATUS, " graph model built");
290 }
291 
292 /*
293 NOTE:
294  Assumes installed PuLP library is version pulp-0.2
295 */
296 template <typename Adapter>
298  const RCP<PartitioningSolution<Adapter> > &solution
299 )
300 {
301  HELLO;
302 
303  size_t numGlobalParts = solution->getTargetGlobalNumberOfParts();
304 
305  int num_parts = (int)numGlobalParts;
306  //TPL_Traits<int, size_t>::ASSIGN(num_parts, numGlobalParts, env);
307 
308  //#ifdef HAVE_ZOLTAN2_MPI
309  // TODO: XtraPuLP
310 
311  int ierr = 0;
312  int np = problemComm->getSize();
313 
314  // Get number of vertices and edges
315  const size_t modelVerts = model->getLocalNumVertices();
316  const size_t modelEdges = model->getLocalNumEdges();
317  int num_verts = (int)modelVerts;
318  long num_edges = (long)modelEdges;
319  //TPL_Traits<int, size_t>::ASSIGN(num_verts, modelVerts, env);
320  //TPL_Traits<long, size_t>::ASSIGN(num_edges, modelEdges, env);
321 
322  // Get vertex info
323  ArrayView<const gno_t> vtxIDs;
324  ArrayView<StridedData<lno_t, scalar_t> > vwgts;
325  size_t nVtx = model->getVertexList(vtxIDs, vwgts);
326  int nVwgts = model->getNumWeightsPerVertex();
327  if (nVwgts > 1) {
328  std::cerr << "Warning: NumWeightsPerVertex is " << nVwgts
329  << " but PuLP allows only one weight. "
330  << " Zoltan2 will use only the first weight per vertex."
331  << std::endl;
332  }
333 
334  int* vertex_weights = NULL;
335  long vertex_weights_sum = 0;
336  if (nVwgts) {
337  vertex_weights = new int[nVtx];
338  scale_weights(nVtx, vwgts[0], vertex_weights);
339  for (int i = 0; i < num_verts; ++i)
340  vertex_weights_sum += vertex_weights[i];
341  }
342 
343  // Get edge info
344  ArrayView<const gno_t> adjs;
345  ArrayView<const lno_t> offsets;
346  ArrayView<StridedData<lno_t, scalar_t> > ewgts;
347  size_t nEdge = model->getEdgeList(adjs, offsets, ewgts);
348  int nEwgts = model->getNumWeightsPerEdge();
349  if (nEwgts > 1) {
350  std::cerr << "Warning: NumWeightsPerEdge is " << nEwgts
351  << " but PuLP allows only one weight. "
352  << " Zoltan2 will use only the first weight per edge."
353  << std::endl;
354  }
355 
356  int* edge_weights = NULL;
357  if (nEwgts) {
358  edge_weights = new int[nEdge];
359  scale_weights(nEdge, ewgts[0], edge_weights);
360  }
361 
362 #ifndef HAVE_ZOLTAN2_MPI
363  // Create PuLP's graph structure
364  int* out_edges = NULL;
365  long* out_offsets = NULL;
367  TPL_Traits<long, const lno_t>::ASSIGN_ARRAY(&out_offsets, offsets);
368 
369  pulp_graph_t g = {num_verts, num_edges,
370  out_edges, out_offsets,
371  vertex_weights, edge_weights, vertex_weights_sum};
372 
373 #else
374  // Create XtraPuLP's graph structure
375  unsigned long* out_edges = NULL;
376  unsigned long* out_offsets = NULL;
379 
380  const size_t modelVertsGlobal = model->getGlobalNumVertices();
381  const size_t modelEdgesGlobal = model->getGlobalNumEdges();
382  unsigned long num_verts_global = (unsigned long)modelVertsGlobal;
383  unsigned long num_edges_global = (unsigned long)modelEdgesGlobal;
384 
385  unsigned long* global_ids = NULL;
387 
388  ArrayView<size_t> vtxDist;
389  model->getVertexDist(vtxDist);
390  unsigned long* verts_per_rank = new unsigned long[np+1];
391  for (int i = 0; i < np+1; ++i)
392  verts_per_rank[i] = vtxDist[i];
393 
394  dist_graph_t g;
395  create_xtrapulp_dist_graph(&g, num_verts_global, num_edges_global,
396  (unsigned long)num_verts, (unsigned long)num_edges,
397  out_edges, out_offsets, global_ids, verts_per_rank,
398  vertex_weights, edge_weights);
399 #endif
400 
401 
402  // Create array for PuLP to return results in.
403  // Or write directly into solution parts array
404  ArrayRCP<part_t> partList(new part_t[num_verts], 0, num_verts, true);
405  int* parts = NULL;
406  if (num_verts && (sizeof(int) == sizeof(part_t))) {
407  // Can write directly into the solution's memory
408  parts = (int *) partList.getRawPtr();
409  }
410  else {
411  // Can't use solution memory directly; will have to copy later.
412  parts = new int[num_verts];
413  }
414 
415  // TODO
416  // Implement target part sizes
417 
418  // Grab options from parameter list
419  const Teuchos::ParameterList &pl = env->getParameters();
420  const Teuchos::ParameterEntry *pe;
421 
422  // figure out which parts of the algorithm we're going to run
423  // Default to PuLP with BFS init
424  // PuLP - do_edge_min = false, do_maxcut_min = false
425  // PuLP-M - do_edge_bal = true, do_maxcut_min = false
426  // PuLP-MM - do_edge_bal = true/false, do_maxcut_min = true
427  bool do_lp_init = false;
428  bool do_bfs_init = true;
429  bool do_edge_bal = false;
430  bool do_repart = false;
431  bool do_maxcut_min = false;
432  bool verbose_output = false;
433 
434  // Do label propagation initialization instead of bfs?
435  pe = pl.getEntryPtr("pulp_lp_init");
436  if (pe) do_lp_init = pe->getValue(&do_lp_init);
437  if (do_lp_init) do_bfs_init = false;
438 
439  // Now look at additional objective
440  pe = pl.getEntryPtr("pulp_minimize_maxcut");
441  if (pe) {
442  do_maxcut_min = pe->getValue(&do_maxcut_min);
443  // If we're doing the secondary objective,
444  // set the additional constraint as well
445  if (do_maxcut_min) do_edge_bal = true;
446  }
447 
448  pe = pl.getEntryPtr("pulp_do_repart");
449  if (pe) {
450  do_repart = pe->getValue(&do_repart);
451  // Do repartitioning with input parts
452  do_bfs_init = false;
453  do_lp_init = false;
454  // TODO: read in current parts
455  // for (int i = 0; i < num_verts; ++i)
456  // parts[i] = something;
457  }
458 
459  // Now grab vertex and edge imbalances, defaults at 10%
460  double vert_imbalance = 1.1;
461  double edge_imbalance = 1.1;
462 
463  pe = pl.getEntryPtr("pulp_vert_imbalance");
464  if (pe) vert_imbalance = pe->getValue<double>(&vert_imbalance);
465  pe = pl.getEntryPtr("pulp_edge_imbalance");
466  if (pe) {
467  edge_imbalance = pe->getValue<double>(&edge_imbalance);
468  // if manually set edge imbalance, add do_edge_bal flag to true
469  do_edge_bal = 1;
470  }
471 
472  if (vert_imbalance < 1.0)
473  throw std::runtime_error("pulp_vert_imbalance must be '1.0' or greater.");
474  if (edge_imbalance < 1.0)
475  throw std::runtime_error("pulp_edge_imbalance must be '1.0' or greater.");
476 
477  // verbose output?
478  // TODO: fully implement verbose flag throughout PuLP
479  pe = pl.getEntryPtr("pulp_verbose");
480  if (pe) verbose_output = pe->getValue(&verbose_output);
481 
482  // using pulp seed?
483  int pulp_seed = rand();
484  pe = pl.getEntryPtr("pulp_seed");
485  if (pe) pulp_seed = pe->getValue(&pulp_seed);
486 
487  // Create PuLP's partitioning data structure
488  pulp_part_control_t ppc = {vert_imbalance, edge_imbalance,
489  do_lp_init, do_bfs_init, do_repart,
490  do_edge_bal, do_maxcut_min,
491  verbose_output, pulp_seed};
492 
493 
494  if (verbose_output) {
495  printf("procid: %d, n: %i, m: %li, vb: %lf, eb: %lf, p: %i\n",
496  problemComm->getRank(),
497  num_verts, num_edges, vert_imbalance, edge_imbalance, num_parts);
498  }
499 
500  // Call partitioning; result returned in parts array
501 #ifndef HAVE_ZOLTAN2_MPI
502  ierr = pulp_run(&g, &ppc, parts, num_parts);
503 
504  env->globalInputAssertion(__FILE__, __LINE__, "pulp_run",
505  !ierr, BASIC_ASSERTION, problemComm);
506 #else
507  ierr = xtrapulp_run(&g, &ppc, parts, num_parts);
508  env->globalInputAssertion(__FILE__, __LINE__, "xtrapulp_run",
509  !ierr, BASIC_ASSERTION, problemComm);
510 #endif
511 
512 
513 
514  // Load answer into the solution if necessary
515  if ((sizeof(int) != sizeof(part_t)) || (num_verts == 0)) {
516  for (int i = 0; i < num_verts; i++) partList[i] = parts[i];
517  delete [] parts;
518  }
519 
520  solution->setParts(partList);
521 
522  env->memory("Zoltan2-(Xtra)PuLP: After creating solution");
523 
524  // Clean up copies made due to differing data sizes.
525 #ifndef HAVE_ZOLTAN2_MPI
528 #else
532 #endif
533 
534 
535 //#endif // DO NOT HAVE_MPI
536 }
537 
539 // Scale and round scalar_t weights (typically float or double) to
540 // PuLP int
541 // subject to sum(weights) <= max_wgt_sum.
542 // Scale only if deemed necessary.
543 //
544 // Note that we use ceil() instead of round() to avoid
545 // rounding to zero weights.
546 // Based on Zoltan's scale_round_weights, mode 1
547 template <typename Adapter>
548 void AlgPuLP<Adapter>::scale_weights(
549  size_t n,
550  StridedData<typename Adapter::lno_t, typename Adapter::scalar_t> &fwgts,
551  int *iwgts
552 )
553 {
554  const double INT_EPSILON = 1e-5;
555  const double MAX_NUM = 1e9;
556 
557  int nonint = 0;
558  double sum_wgt = 0.0;
559  double max_wgt = 0.0;
560 
561  // Compute local sums of the weights
562  // Check whether all weights are integers
563  for (size_t i = 0; i < n; i++) {
564  double fw = double(fwgts[i]);
565  if (!nonint){
566  int tmp = (int) floor(fw + .5); /* Nearest int */
567  if (fabs((double)tmp-fw) > INT_EPSILON) {
568  nonint = 1;
569  }
570  }
571  sum_wgt += fw;
572  if (fw > max_wgt) max_wgt = fw;
573  }
574 
575  // Scaling needed if weights are not integers or weights'
576  // range is not sufficient
577  double scale = 1.0;
578  if (nonint || (max_wgt <= INT_EPSILON) || (sum_wgt > MAX_NUM)) {
579  /* Calculate scale factor */
580  if (sum_wgt != 0.0) scale = MAX_NUM/sum_wgt;
581  }
582 
583  /* Convert weights to positive integers using the computed scale factor */
584  for (size_t i = 0; i < n; i++)
585  iwgts[i] = (int) ceil(double(fwgts[i])*scale);
586 
587 }
588 
589 
590 } // namespace Zoltan2
591 
592 #endif // HAVE_ZOLTAN2_PULP
593 
595 
596 
597 #endif
598 
599 
600 
601 
602 
603 
use columns as graph vertices
#define HELLO
Zoltan2::BaseAdapter< userTypes_t > base_adapter_t
ignore invalid neighbors
use mesh nodes as vertices
fast typical checks for valid arguments
Adapter::base_adapter_t base_adapter_t
static RCP< Teuchos::BoolParameterEntryValidator > getBoolValidator()
Exists to make setting up validators less cluttered.
algorithm requires consecutive ids
std::bitset< NUM_MODEL_FLAGS > modelFlag_t
model must symmetrize input
model must symmetrize input
Defines the PartitioningSolution class.
static RCP< Teuchos::AnyNumberParameterEntryValidator > getAnyDoubleValidator()
Exists to make setting up validators less cluttered.
sub-steps, each method&#39;s entry and exit
use matrix rows as graph vertices
algorithm requires no self edges
Adapter::scalar_t scalar_t
use nonzeros as graph vertices
Adapter::part_t part_t
Algorithm defines the base class for all algorithms.
static RCP< Teuchos::AnyNumberParameterEntryValidator > getAnyIntValidator()
Exists to make setting up validators less cluttered.
virtual void partition(const RCP< PartitioningSolution< Adapter > > &solution)
Partitioning method.
Traits class to handle conversions between gno_t/lno_t and TPL data types (e.g., ParMETIS&#39;s idx_t...
static void getValidParameters(ParameterList &pl)
Set up validators specific to this algorithm.
use mesh elements as vertices
static void ASSIGN_ARRAY(first_t **a, ArrayView< second_t > &b)
Defines the GraphModel interface.
AlgPuLP(const RCP< const Environment > &env, const RCP< const Comm< int > > &problemComm, const RCP< const base_adapter_t > &adapter)
A gathering of useful namespace methods.
static void DELETE_ARRAY(first_t **a)
model represents graph within only one rank