Zoltan2
Zoltan2_AlgParMETIS.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_ALGPARMETIS_HPP_
46 #define _ZOLTAN2_ALGPARMETIS_HPP_
47 
48 #include <Zoltan2_GraphModel.hpp>
49 #include <Zoltan2_Algorithm.hpp>
51 #include <Zoltan2_Util.hpp>
52 
57 
58 #ifndef HAVE_ZOLTAN2_PARMETIS
59 
60 // Error handling for when ParMETIS is requested
61 // but Zoltan2 not built with ParMETIS.
62 
63 namespace Zoltan2 {
64 template <typename Adapter>
65 class AlgParMETIS : public Algorithm<Adapter>
66 {
67 public:
68  AlgParMETIS(const RCP<const Environment> &env,
69  const RCP<const Comm<int> > &problemComm,
71  )
72  {
73  throw std::runtime_error(
74  "BUILD ERROR: ParMETIS requested but not compiled into Zoltan2.\n"
75  "Please set CMake flag Zoltan2_ENABLE_ParMETIS:BOOL=ON.");
76  }
77 };
78 }
79 
80 #endif
81 
84 
85 #ifdef HAVE_ZOLTAN2_PARMETIS
86 
87 #ifndef HAVE_ZOLTAN2_MPI
88 
89 // ParMETIS requires compilation with MPI.
90 // If MPI is not available, make compilation fail.
91 #error "TPL ParMETIS requires compilation with MPI. Configure with -DTPL_ENABLE_MPI:BOOL=ON or -DZoltan2_ENABLE_ParMETIS:BOOL=OFF"
92 
93 #else
94 
95 extern "C" {
96 #include "parmetis.h"
97 }
98 
99 #if (PARMETIS_MAJOR_VERSION < 4)
100 
101 // Zoltan2 requires ParMETIS v4.x.
102 // Make compilation fail for earlier versions of ParMETIS.
103 #error "Specified version of ParMETIS is not compatible with Zoltan2; upgrade to ParMETIS v4 or later, or build Zoltan2 without ParMETIS."
104 
105 #else
106 
107 // MPI and ParMETIS version requirements are met. Proceed.
108 
109 namespace Zoltan2 {
110 
111 template <typename Adapter>
112 class AlgParMETIS : public Algorithm<Adapter>
113 {
114 public:
115 
116  typedef GraphModel<typename Adapter::base_adapter_t> graphModel_t;
117  typedef typename Adapter::lno_t lno_t;
118  typedef typename Adapter::gno_t gno_t;
119  typedef typename Adapter::scalar_t scalar_t;
120  typedef typename Adapter::part_t part_t;
121 
122  typedef idx_t pm_idx_t;
123  typedef real_t pm_real_t;
124 
135  AlgParMETIS(const RCP<const Environment> &env__,
136  const RCP<const Comm<int> > &problemComm__,
137  const RCP<graphModel_t> &model__) :
138  env(env__), problemComm(problemComm__),
139  model(model__)
140  { }
141 
142  void partition(const RCP<PartitioningSolution<Adapter> > &solution);
143 
144 private:
145 
146  const RCP<const Environment> env;
147  const RCP<const Comm<int> > problemComm;
148  const RCP<GraphModel<typename Adapter::base_adapter_t> > model;
149 
150  void scale_weights(size_t n, ArrayView<StridedData<lno_t, scalar_t> > &fwgts,
151  pm_idx_t *iwgts);
152 };
153 
154 
156 template <typename Adapter>
158  const RCP<PartitioningSolution<Adapter> > &solution
159 )
160 {
161  HELLO;
162 
163  size_t numGlobalParts = solution->getTargetGlobalNumberOfParts();
164 
165  int me = problemComm->getRank();
166  int np = problemComm->getSize();
167 
168  // Get vertex info
169  ArrayView<const gno_t> vtxgnos;
170  ArrayView<StridedData<lno_t, scalar_t> > vwgts;
171  int nVwgt = model->getNumWeightsPerVertex();
172  size_t nVtx = model->getVertexList(vtxgnos, vwgts);
173  pm_idx_t pm_nVtx;
175 
176  pm_idx_t *pm_vwgts = NULL;
177  if (nVwgt) {
178  pm_vwgts = new pm_idx_t[nVtx*nVwgt];
179  scale_weights(nVtx, vwgts, pm_vwgts);
180  }
181 
182  // Get edge info
183  ArrayView<const gno_t> adjgnos;
184  ArrayView<const lno_t> offsets;
185  ArrayView<StridedData<lno_t, scalar_t> > ewgts;
186  int nEwgt = model->getNumWeightsPerEdge();
187  size_t nEdge = model->getEdgeList(adjgnos, offsets, ewgts);
188 
189  pm_idx_t *pm_ewgts = NULL;
190  if (nEwgt) {
191  pm_ewgts = new pm_idx_t[nEdge*nEwgt];
192  scale_weights(nEdge, ewgts, pm_ewgts);
193  }
194 
195  // Convert index types for edges, if needed
196  pm_idx_t *pm_offsets;
198  pm_idx_t *pm_adjs;
199  pm_idx_t pm_dummy_adj;
200  if (nEdge)
202  else
203  pm_adjs = &pm_dummy_adj; // ParMETIS does not like NULL pm_adjs;
204 
205 
206  // Build vtxdist
207  pm_idx_t *pm_vtxdist;
208  ArrayView<size_t> vtxdist;
209  model->getVertexDist(vtxdist);
210  TPL_Traits<pm_idx_t,size_t>::ASSIGN_ARRAY(&pm_vtxdist, vtxdist);
211 
212  // ParMETIS does not like processors having no vertices.
213  // Inspect vtxdist and remove from communicator procs that have no vertices
214  RCP<Comm<int> > subcomm;
215  MPI_Comm mpicomm; // Note: mpicomm is valid only while subcomm is in scope
216 
217  int nKeep = 0;
218  if (np > 1) {
219  Array<int> keepRanks(np);
220  for (int i = 0; i < np; i++) {
221  if ((pm_vtxdist[i+1] - pm_vtxdist[i]) > 0) {
222  keepRanks[nKeep] = i;
223  pm_vtxdist[nKeep] = pm_vtxdist[i];
224  nKeep++;
225  }
226  }
227  pm_vtxdist[nKeep] = pm_vtxdist[np];
228  if (nKeep < np) {
229  subcomm = problemComm->createSubcommunicator(keepRanks.view(0,nKeep));
230  if (subcomm != Teuchos::null)
231  mpicomm = Teuchos::getRawMpiComm(*subcomm);
232  else
233  mpicomm = MPI_COMM_NULL;
234  }
235  else {
236  mpicomm = Teuchos::getRawMpiComm(*problemComm);
237  }
238  }
239  else {
240  mpicomm = Teuchos::getRawMpiComm(*problemComm);
241  }
242 
243  // Create array for ParMETIS to return results in.
244  pm_idx_t *pm_partList = NULL;
245  if (nVtx) pm_partList = new pm_idx_t[nVtx];
246  for (size_t i = 0; i < nVtx; i++) pm_partList[i] = 0;
247  int pm_return = METIS_OK;
248 
249  if (mpicomm != MPI_COMM_NULL) {
250  // If in ParMETIS' communicator (i.e., have vertices), call ParMETIS
251 
252  // Get target part sizes
253  pm_idx_t pm_nCon = (nVwgt == 0 ? 1 : pm_idx_t(nVwgt));
254  pm_real_t *pm_partsizes = new pm_real_t[numGlobalParts*pm_nCon];
255  for (pm_idx_t dim = 0; dim < pm_nCon; dim++) {
256  if (!solution->criteriaHasUniformPartSizes(dim))
257  for (size_t i=0; i<numGlobalParts; i++)
258  pm_partsizes[i*pm_nCon+dim] =
259  pm_real_t(solution->getCriteriaPartSize(dim,i));
260  else
261  for (size_t i=0; i<numGlobalParts; i++)
262  pm_partsizes[i*pm_nCon+dim] = pm_real_t(1.)/pm_real_t(numGlobalParts);
263  }
264 
265  // Get imbalance tolerances
266  double tolerance = 1.1;
267  const Teuchos::ParameterList &pl = env->getParameters();
268  const Teuchos::ParameterEntry *pe = pl.getEntryPtr("imbalance_tolerance");
269  if (pe) tolerance = pe->getValue<double>(&tolerance);
270 
271  // ParMETIS requires tolerance to be greater than 1.0;
272  // fudge it if condition is not met
273  if (tolerance <= 1.0) {
274  if (me == 0)
275  std::cerr << "Warning: ParMETIS requires imbalance_tolerance > 1.0; "
276  << "to comply, Zoltan2 reset imbalance_tolerance to 1.01."
277  << std::endl;
278  tolerance = 1.01;
279  }
280 
281  pm_real_t *pm_imbTols = new pm_real_t[pm_nCon];
282  for (pm_idx_t dim = 0; dim < pm_nCon; dim++)
283  pm_imbTols[dim] = pm_real_t(tolerance);
284 
285  std::string parmetis_method("PARTKWAY");
286  pe = pl.getEntryPtr("partitioning_approach");
287  if (pe){
288  std::string approach;
289  approach = pe->getValue<std::string>(&approach);
290  if ((approach == "repartition") || (approach == "maximize_overlap")) {
291  if (nKeep > 1)
292  // ParMETIS_V3_AdaptiveRepart requires two or more processors
293  parmetis_method = "ADAPTIVE_REPART";
294  else
295  // Probably best to do PartKway if nKeep == 1;
296  // I think REFINE_KWAY won't give a good answer in most use cases
297  // parmetis_method = "REFINE_KWAY";
298  parmetis_method = "PARTKWAY";
299  }
300  }
301 
302  // Other ParMETIS parameters?
303  pm_idx_t pm_wgtflag = 2*(nVwgt > 0) + (nEwgt > 0);
304  pm_idx_t pm_numflag = 0;
305  pm_idx_t pm_edgecut = -1;
306  pm_idx_t pm_options[METIS_NOPTIONS];
307  pm_options[0] = 1; // Use non-default options for some ParMETIS options
308  for (int i = 0; i < METIS_NOPTIONS; i++)
309  pm_options[i] = 0; // Default options
310  pm_options[2] = 15; // Matches default value used in Zoltan
311 
312  pm_idx_t pm_nPart;
313  TPL_Traits<pm_idx_t,size_t>::ASSIGN(pm_nPart, numGlobalParts);
314 
315  if (parmetis_method == "PARTKWAY") {
316  pm_return = ParMETIS_V3_PartKway(pm_vtxdist, pm_offsets, pm_adjs,
317  pm_vwgts, pm_ewgts, &pm_wgtflag,
318  &pm_numflag, &pm_nCon, &pm_nPart,
319  pm_partsizes, pm_imbTols, pm_options,
320  &pm_edgecut, pm_partList, &mpicomm);
321  }
322  else if (parmetis_method == "ADAPTIVE_REPART") {
323  // Get object sizes: pm_vsize
324  // TODO: get pm_vsize info from input adapter or graph model
325  // TODO: This is just a placeholder
326  pm_idx_t *pm_vsize = new pm_idx_t[nVtx];
327  for (size_t i = 0; i < nVtx; i++) pm_vsize[i] = 1;
328 
329  pm_real_t itr = 100.; // Same default as in Zoltan
330  pm_return = ParMETIS_V3_AdaptiveRepart(pm_vtxdist, pm_offsets, pm_adjs,
331  pm_vwgts,
332  pm_vsize, pm_ewgts, &pm_wgtflag,
333  &pm_numflag, &pm_nCon, &pm_nPart,
334  pm_partsizes, pm_imbTols,
335  &itr, pm_options, &pm_edgecut,
336  pm_partList, &mpicomm);
337  delete [] pm_vsize;
338  }
339  // else if (parmetis_method == "REFINE_KWAY") {
340  // We do not currently have an execution path that calls REFINE_KWAY.
341  // pm_return = ParMETIS_V3_RefineKway(pm_vtxdist, pm_offsets, pm_adjs,
342  // pm_vwgts, pm_ewgts, &pm_wgtflag,
343  // &pm_numflag, &pm_nCon, &pm_nPart,
344  // pm_partsizes, pm_imbTols, pm_options,
345  // &pm_edgecut, pm_partList, &mpicomm);
346  // }
347  else {
348  // We should not reach this condition.
349  throw std::logic_error("\nInvalid ParMETIS method requested.\n");
350  }
351 
352  // Clean up
353  delete [] pm_partsizes;
354  delete [] pm_imbTols;
355  }
356 
357  // Load answer into the solution.
358 
359  ArrayRCP<part_t> partList;
360  if (nVtx)
361  TPL_Traits<part_t, pm_idx_t>::SAVE_ARRAYRCP(&partList, pm_partList, nVtx);
363 
364  solution->setParts(partList);
365 
366  env->memory("Zoltan2-ParMETIS: After creating solution");
367 
368  // Clean up copies made due to differing data sizes.
371  if (nEdge)
373 
374  if (nVwgt) delete [] pm_vwgts;
375  if (nEwgt) delete [] pm_ewgts;
376 
377  if (pm_return != METIS_OK) {
378  throw std::runtime_error(
379  "\nParMETIS returned an error; no valid partition generated.\n"
380  "Look for 'PARMETIS ERROR' in your output for more details.\n");
381  }
382 }
383 
385 // Scale and round scalar_t weights (typically float or double) to
386 // ParMETIS' idx_t (typically int or long).
387 // subject to sum(weights) <= max_wgt_sum.
388 // Scale only if deemed necessary.
389 //
390 // Note that we use ceil() instead of round() to avoid
391 // rounding to zero weights.
392 // Based on Zoltan's scale_round_weights, mode 1
393 
394 template <typename Adapter>
395 void AlgParMETIS<Adapter>::scale_weights(
396  size_t n,
397  ArrayView<StridedData<typename Adapter::lno_t,
398  typename Adapter::scalar_t> > &fwgts,
399  pm_idx_t *iwgts
400 )
401 {
402  const double INT_EPSILON = 1e-5;
403  const int nWgt = fwgts.size();
404 
405  int *nonint_local = new int[nWgt+nWgt];
406  int *nonint = nonint_local + nWgt;
407 
408  double *sum_wgt_local = new double[nWgt*4];
409  double *max_wgt_local = sum_wgt_local + nWgt;
410  double *sum_wgt = max_wgt_local + nWgt;
411  double *max_wgt = sum_wgt + nWgt;
412 
413  for (int i = 0; i < nWgt; i++) {
414  nonint_local[i] = 0;
415  sum_wgt_local[i] = 0.;
416  max_wgt_local[i] = 0;
417  }
418 
419  // Compute local sums of the weights
420  // Check whether all weights are integers
421  for (int j = 0; j < nWgt; j++) {
422  for (size_t i = 0; i < n; i++) {
423  double fw = double(fwgts[j][i]);
424  if (!nonint_local[j]) {
425  pm_idx_t tmp = (pm_idx_t) floor(fw + .5); /* Nearest int */
426  if (fabs((double)tmp-fw) > INT_EPSILON) {
427  nonint_local[j] = 1;
428  }
429  }
430  sum_wgt_local[j] += fw;
431  if (fw > max_wgt_local[j]) max_wgt_local[j] = fw;
432  }
433  }
434 
435  Teuchos::reduceAll<int,int>(*problemComm, Teuchos::REDUCE_MAX, nWgt,
436  nonint_local, nonint);
437  Teuchos::reduceAll<int,double>(*problemComm, Teuchos::REDUCE_SUM, nWgt,
438  sum_wgt_local, sum_wgt);
439  Teuchos::reduceAll<int,double>(*problemComm, Teuchos::REDUCE_MAX, nWgt,
440  max_wgt_local, max_wgt);
441 
442  const double max_wgt_sum = double(std::numeric_limits<pm_idx_t>::max()/8);
443  for (int j = 0; j < nWgt; j++) {
444  double scale = 1.;
445 
446  // Scaling needed if weights are not integers or weights'
447  // range is not sufficient
448  if (nonint[j] || (max_wgt[j]<=INT_EPSILON) || (sum_wgt[j]>max_wgt_sum)) {
449  /* Calculate scale factor */
450  if (sum_wgt[j] != 0.) scale = max_wgt_sum/sum_wgt[j];
451  }
452 
453  /* Convert weights to positive integers using the computed scale factor */
454  for (size_t i = 0; i < n; i++)
455  iwgts[i*nWgt+j] = (pm_idx_t) ceil(double(fwgts[j][i])*scale);
456  }
457  delete [] nonint_local;
458  delete [] sum_wgt_local;
459 }
460 
461 } // namespace Zoltan2
462 
463 #endif // PARMETIS VERSION 4 OR HIGHER CHECK
464 
465 #endif // HAVE_ZOLTAN2_MPI
466 
467 #endif // HAVE_ZOLTAN2_PARMETIS
468 
469 #endif
#define HELLO
Defines the PartitioningSolution class.
static void SAVE_ARRAYRCP(ArrayRCP< first_t > *a, second_t *b, size_t size)
AlgParMETIS(const RCP< const Environment > &env, const RCP< const Comm< int > > &problemComm, const RCP< GraphModel< typename Adapter::base_adapter_t > > &model)
Adapter::scalar_t scalar_t
Adapter::part_t part_t
Algorithm defines the base class for all algorithms.
static void ASSIGN(first_t &a, second_t b)
virtual void partition(const RCP< PartitioningSolution< Adapter > > &solution)
Partitioning method.
GraphModel defines the interface required for graph models.
static void ASSIGN_ARRAY(first_t **a, ArrayView< second_t > &b)
Defines the GraphModel interface.
A gathering of useful namespace methods.
static void DELETE_ARRAY(first_t **a)