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