Zoltan2
Loading...
Searching...
No Matches
Zoltan2_AlgScotch.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Zoltan2: A package of combinatorial algorithms for scientific computing
4//
5// Copyright 2012 NTESS and the Zoltan2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef _ZOLTAN2_ALGSCOTCH_HPP_
11#define _ZOLTAN2_ALGSCOTCH_HPP_
12
14#include <Zoltan2_Algorithm.hpp>
16#include <Zoltan2_OrderingSolution.hpp> // BDD: needed by ordering method
17#include <Zoltan2_Util.hpp>
18#include <Zoltan2_TPLTraits.hpp>
19
23
25
26namespace Zoltan2 {
27
28// this function called by the two scotch types below
29static inline void getScotchParameters(Teuchos::ParameterList & pl)
30{
31 // bool parameter
32 pl.set("scotch_verbose", false, "verbose output",
34
35 RCP<Teuchos::EnhancedNumberValidator<int>> scotch_level_Validator =
36 Teuchos::rcp( new Teuchos::EnhancedNumberValidator<int>(0, 1000, 1, 0) );
37 pl.set("scotch_level", 0, "scotch ordering - Level of the subgraph in the "
38 "separators tree for the initial graph at the root of the tree",
39 scotch_level_Validator);
40
41 pl.set("scotch_imbalance_ratio", 0.2, "scotch ordering - Dissection "
42 "imbalance ratio", Environment::getAnyDoubleValidator());
43
44 // bool parameter
45 pl.set("scotch_ordering_default", true, "use default scotch ordering "
46 "strategy", Environment::getBoolValidator());
47
48 pl.set("scotch_ordering_strategy", "", "scotch ordering - Dissection "
49 "imbalance ratio");
50}
51
52} // Zoltan2 namespace
53
54#ifndef HAVE_ZOLTAN2_SCOTCH
55
56namespace Zoltan2 {
57
58// Error handling for when Scotch is requested
59// but Zoltan2 not built with Scotch.
60
61template <typename Adapter>
62class AlgPTScotch : public Algorithm<Adapter>
63{
64public:
65 typedef typename Adapter::base_adapter_t base_adapter_t;
66 //AlgPTScotch(const RCP<const Environment> &env,
67 // const RCP<const Comm<int> > &problemComm,
68 // const RCP<GraphModel<typename Adapter::base_adapter_t> > &model
69 //) BDD: old inteface for reference
70 AlgPTScotch(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: Scotch requested but not compiled into Zoltan2.\n"
77 "Please set CMake flag Zoltan2_ENABLE_Scotch:BOOL=ON.");
78 }
79
82 static void getValidParameters(ParameterList & pl)
83 {
85 }
86};
87}
88#endif
89
92
93#ifdef HAVE_ZOLTAN2_SCOTCH
94
95namespace Zoltan2 {
96
97// stdint.h for int64_t in scotch header
98#include <stdint.h>
99extern "C" {
100#ifndef HAVE_ZOLTAN2_MPI
101#include "scotch.h"
102#else
103#include "ptscotch.h"
104#endif
105}
106
107#ifdef SHOW_ZOLTAN2_SCOTCH_MEMORY
108//
109// Scotch keeps track of memory high water mark, but doesn't
110// provide a way to get that number. So add this function:
111// "size_t SCOTCH_getMemoryMax() { return memorymax;}"
112// to src/libscotch/common_memory.c
113//
114// and this macro:
115// "#define HAVE_SCOTCH_GETMEMORYMAX
116// to include/ptscotch.h
117//
118// and compile scotch with -DCOMMON_MEMORY_TRACE
119//
120#ifdef HAVE_SCOTCH_GETMEMORYMAX
121 extern "C"{
122 extern size_t SCOTCH_getMemoryMax();
123 }
124#else
125#error "Either turn off ZOLTAN2_ENABLE_SCOTCH_MEMORY_REPORT in cmake configure, or see SHOW_ZOLTAN2_SCOTCH_MEMORY comment in Zoltan2_AlgScotch.hpp"
126#endif // HAVE_SCOTCH_GETMEMORYMAX
127#endif // SHOW_ZOLTAN2_SCOTCH_MEMORY
128
129template <typename Adapter>
130class AlgPTScotch : public Algorithm<Adapter>
131{
132public:
133
134 typedef typename Adapter::base_adapter_t base_adapter_t;
135 typedef GraphModel<typename Adapter::base_adapter_t> graphModel_t;
136 typedef typename Adapter::lno_t lno_t;
137 typedef typename Adapter::gno_t gno_t;
138 typedef typename Adapter::offset_t offset_t;
139 typedef typename Adapter::scalar_t scalar_t;
140 typedef typename Adapter::part_t part_t;
141 typedef typename Adapter::user_t user_t;
142 typedef typename Adapter::userCoord_t userCoord_t;
143
144// /*! Scotch constructor
145// * \param env parameters for the problem and library configuration
146// * \param problemComm the communicator for the problem
147// * \param model a graph
148// *
149// * Preconditions: The parameters in the environment have been processed.
150// * TODO: THIS IS A MINIMAL CONSTRUCTOR FOR NOW.
151// * TODO: WHEN ADD SCOTCH ORDERING OR MAPPING, MOVE SCOTCH GRAPH CONSTRUCTION
152// * TODO: TO THE CONSTRUCTOR SO THAT CODE MAY BE SHARED.
153// */
154// AlgPTScotch(const RCP<const Environment> &env__,
155// const RCP<const Comm<int> > &problemComm__,
156// const RCP<graphModel_t> &model__) :
157// env(env__), problemComm(problemComm__),
158//#ifdef HAVE_ZOLTAN2_MPI
159// mpicomm(Teuchos::getRawMpiComm(*problemComm__)),
160//#endif
161// model(model__) BDD: olde interface for reference
162
173 AlgPTScotch(const RCP<const Environment> &env__,
174 const RCP<const Comm<int> > &problemComm__,
175 const RCP<const IdentifierAdapter<user_t> > &adapter__) :
176 env(env__), problemComm(problemComm__),
177#ifdef HAVE_ZOLTAN2_MPI
178 mpicomm(Teuchos::getRawMpiComm(*problemComm__)),
179#endif
180 adapter(adapter__)
181 {
182 std::string errStr = "cannot build GraphModel from IdentifierAdapter, ";
183 errStr += "Scotch requires Graph, Matrix, or Mesh Adapter";
184 throw std::runtime_error(errStr);
185 }
186
187 AlgPTScotch(const RCP<const Environment> &env__,
188 const RCP<const Comm<int> > &problemComm__,
189 const RCP<const VectorAdapter<user_t> > &adapter__) :
190 env(env__), problemComm(problemComm__),
191#ifdef HAVE_ZOLTAN2_MPI
192 mpicomm(Teuchos::getRawMpiComm(*problemComm__)),
193#endif
194 adapter(adapter__)
195 {
196 std::string errStr = "cannot build GraphModel from IdentifierAdapter, ";
197 errStr += "Scotch requires Graph, Matrix, or Mesh Adapter";
198 throw std::runtime_error(errStr);
199 }
200
201 AlgPTScotch(const RCP<const Environment> &env__,
202 const RCP<const Comm<int> > &problemComm__,
203 const RCP<const GraphAdapter<user_t, userCoord_t> > &adapter__) :
204 env(env__), problemComm(problemComm__),
205#ifdef HAVE_ZOLTAN2_MPI
206 mpicomm(Teuchos::getRawMpiComm(*problemComm__)),
207#endif
208 adapter(adapter__), model_flags()
209 {
210 this->model_flags.reset();
211 }
212
213 AlgPTScotch(const RCP<const Environment> &env__,
214 const RCP<const Comm<int> > &problemComm__,
215 const RCP<const MatrixAdapter<user_t, userCoord_t> > &adapter__) :
216 env(env__), problemComm(problemComm__),
217#ifdef HAVE_ZOLTAN2_MPI
218 mpicomm(Teuchos::getRawMpiComm(*problemComm__)),
219#endif
220 adapter(adapter__), model_flags()
221 {
222 this->model_flags.reset();
223
224 const ParameterList &pl = env->getParameters();
225 const Teuchos::ParameterEntry *pe;
226
227 std::string defString("default");
228 std::string objectOfInterest(defString);
229 pe = pl.getEntryPtr("objects_to_partition");
230 if (pe)
231 objectOfInterest = pe->getValue<std::string>(&objectOfInterest);
232
233 if (objectOfInterest == defString ||
234 objectOfInterest == std::string("matrix_rows") )
235 this->model_flags.set(VERTICES_ARE_MATRIX_ROWS);
236 else if (objectOfInterest == std::string("matrix_columns"))
237 this->model_flags.set(VERTICES_ARE_MATRIX_COLUMNS);
238 else if (objectOfInterest == std::string("matrix_nonzeros"))
239 this->model_flags.set(VERTICES_ARE_MATRIX_NONZEROS);
240 }
241
242 AlgPTScotch(const RCP<const Environment> &env__,
243 const RCP<const Comm<int> > &problemComm__,
244 const RCP<const MeshAdapter<user_t> > &adapter__) :
245 env(env__), problemComm(problemComm__),
246#ifdef HAVE_ZOLTAN2_MPI
247 mpicomm(Teuchos::getRawMpiComm(*problemComm__)),
248#endif
249 adapter(adapter__), model_flags()
250 {
251 this->model_flags.reset();
252
253 const ParameterList &pl = env->getParameters();
254 const Teuchos::ParameterEntry *pe;
255
256 std::string defString("default");
257 std::string objectOfInterest(defString);
258 pe = pl.getEntryPtr("objects_to_partition");
259 if (pe)
260 objectOfInterest = pe->getValue<std::string>(&objectOfInterest);
261
262 if (objectOfInterest == defString ||
263 objectOfInterest == std::string("mesh_nodes") )
264 this->model_flags.set(VERTICES_ARE_MESH_NODES);
265 else if (objectOfInterest == std::string("mesh_elements"))
266 this->model_flags.set(VERTICES_ARE_MESH_ELEMENTS);
267 }
268
271 static void getValidParameters(ParameterList & pl)
272 {
274 }
275
276 void partition(const RCP<PartitioningSolution<Adapter> > &solution);
277
278 int localOrder(const RCP<LocalOrderingSolution<lno_t> > &solution);
279 int globalOrder(const RCP<GlobalOrderingSolution<gno_t> > &solution);
280
281private:
282
283 const RCP<const Environment> env;
284 const RCP<const Comm<int> > problemComm;
285#ifdef HAVE_ZOLTAN2_MPI
286 const MPI_Comm mpicomm;
287#endif
288 const RCP<const base_adapter_t> adapter;
289 modelFlag_t model_flags;
290 RCP<graphModel_t > model; // BDD:to be constructed
291
292 void buildModel(bool isLocal);
293 void scale_weights(size_t n, StridedData<lno_t, scalar_t> &fwgts,
294 SCOTCH_Num *iwgts);
295 static int setStrategy(SCOTCH_Strat * c_strat_ptr, const ParameterList &pList, int procRank);
296};
297
299template <typename Adapter>
300void AlgPTScotch<Adapter>::buildModel(bool isLocal) {
301 HELLO;
302 const ParameterList &pl = env->getParameters();
303 const Teuchos::ParameterEntry *pe;
304
305 std::string defString("default");
306 std::string symParameter(defString);
307 pe = pl.getEntryPtr("symmetrize_graph");
308 if (pe){
309 symParameter = pe->getValue<std::string>(&symParameter);
310 if (symParameter == std::string("transpose"))
311 this->model_flags.set(SYMMETRIZE_INPUT_TRANSPOSE);
312 else if (symParameter == std::string("bipartite"))
313 this->model_flags.set(SYMMETRIZE_INPUT_BIPARTITE); }
314
315 bool sgParameter = false;
316 pe = pl.getEntryPtr("subset_graph");
317 if (pe)
318 sgParameter = pe->getValue(&sgParameter);
319 if (sgParameter)
320 this->model_flags.set(BUILD_SUBSET_GRAPH);
321
322 this->model_flags.set(REMOVE_SELF_EDGES);
323 this->model_flags.set(GENERATE_CONSECUTIVE_IDS);
324 if (isLocal) {
325 HELLO; // only for ordering!
326 this->model_flags.set(BUILD_LOCAL_GRAPH);
327 }
328
329 this->env->debug(DETAILED_STATUS, " building graph model");
330 this->model = rcp(new graphModel_t(this->adapter,
331 this->env,
332 this->problemComm,
333 this->model_flags));
334 this->env->debug(DETAILED_STATUS, " graph model built");
335}
336
337template <typename Adapter>
339 const RCP<PartitioningSolution<Adapter> > &solution
340)
341{
342 HELLO;
343 this->buildModel(false);
344
345 size_t numGlobalParts = solution->getTargetGlobalNumberOfParts();
346
347 SCOTCH_Num partnbr=0;
348 TPL_Traits<SCOTCH_Num, size_t>::ASSIGN(partnbr, numGlobalParts);
349
350#ifdef HAVE_ZOLTAN2_MPI
351 int ierr = 0;
352 int me = problemComm->getRank();
353
354 const SCOTCH_Num baseval = 0; // Base value for array indexing.
355 // GraphModel returns GNOs from base 0.
356
357 SCOTCH_Strat stratstr; // Strategy string
358 // TODO: Set from parameters
359 SCOTCH_stratInit(&stratstr);
360
361 // Allocate and initialize PTScotch Graph data structure.
362 SCOTCH_Dgraph *gr = SCOTCH_dgraphAlloc(); // Scotch distributed graph
363 ierr = SCOTCH_dgraphInit(gr, mpicomm);
364
365 env->globalInputAssertion(__FILE__, __LINE__, "SCOTCH_dgraphInit",
366 !ierr, BASIC_ASSERTION, problemComm);
367
368 // Get vertex info
369 ArrayView<const gno_t> vtxID;
370 ArrayView<StridedData<lno_t, scalar_t> > vwgts;
371 size_t nVtx = model->getVertexList(vtxID, vwgts);
372 SCOTCH_Num vertlocnbr=0;
374 SCOTCH_Num vertlocmax = vertlocnbr; // Assumes no holes in global nums.
375
376 // Get edge info
377 ArrayView<const gno_t> edgeIds;
378 ArrayView<const offset_t> offsets;
379 ArrayView<StridedData<lno_t, scalar_t> > ewgts;
380
381 size_t nEdge = model->getEdgeList(edgeIds, offsets, ewgts);
382
383 SCOTCH_Num edgelocnbr=0;
385 const SCOTCH_Num edgelocsize = edgelocnbr; // Assumes adj array is compact.
386
387 SCOTCH_Num *vertloctab; // starting adj/vtx
389
390 SCOTCH_Num *edgeloctab; // adjacencies
392
393 // We don't use these arrays, but we need them as arguments to Scotch.
394 SCOTCH_Num *vendloctab = NULL; // Assume consecutive storage for adj
395 SCOTCH_Num *vlblloctab = NULL; // Vertex label array
396 SCOTCH_Num *edgegsttab = NULL; // Array for ghost vertices
397
398 // Get weight info.
399 SCOTCH_Num *velotab = NULL; // Vertex weights
400 SCOTCH_Num *edlotab = NULL; // Edge weights
401
402 int nVwgts = model->getNumWeightsPerVertex();
403 int nEwgts = model->getNumWeightsPerEdge();
404 if (nVwgts > 1 && me == 0) {
405 std::cerr << "Warning: NumWeightsPerVertex is " << nVwgts
406 << " but Scotch allows only one weight. "
407 << " Zoltan2 will use only the first weight per vertex."
408 << std::endl;
409 }
410 if (nEwgts > 1 && me == 0) {
411 std::cerr << "Warning: NumWeightsPerEdge is " << nEwgts
412 << " but Scotch allows only one weight. "
413 << " Zoltan2 will use only the first weight per edge."
414 << std::endl;
415 }
416
417 if (nVwgts) {
418 velotab = new SCOTCH_Num[nVtx+1]; // +1 since Scotch wants all procs
419 // to have non-NULL arrays
420 scale_weights(nVtx, vwgts[0], velotab);
421 }
422
423 if (nEwgts) {
424 edlotab = new SCOTCH_Num[nEdge+1]; // +1 since Scotch wants all procs
425 // to have non-NULL arrays
426 scale_weights(nEdge, ewgts[0], edlotab);
427 }
428
429 // Build PTScotch distributed data structure
430 ierr = SCOTCH_dgraphBuild(gr, baseval, vertlocnbr, vertlocmax,
431 vertloctab, vendloctab, velotab, vlblloctab,
432 edgelocnbr, edgelocsize,
433 edgeloctab, edgegsttab, edlotab);
434
435 env->globalInputAssertion(__FILE__, __LINE__, "SCOTCH_dgraphBuild",
436 !ierr, BASIC_ASSERTION, problemComm);
437
438 // Create array for Scotch to return results in.
439 SCOTCH_Num *partloctab = new SCOTCH_Num[nVtx+1];
440 // Note: Scotch does not like NULL arrays, so add 1 to always have non-null.
441 // ParMETIS has this same "feature." See Zoltan bug 4299.
442
443 // Get target part sizes
444 float *partsizes = new float[numGlobalParts];
445 if (!solution->criteriaHasUniformPartSizes(0))
446 for (size_t i=0; i<numGlobalParts; i++)
447 partsizes[i] = solution->getCriteriaPartSize(0, i);
448 else
449 for (size_t i=0; i<numGlobalParts; i++)
450 partsizes[i] = 1.0 / float(numGlobalParts);
451
452 // Allocate and initialize PTScotch target architecture data structure
453 SCOTCH_Arch archdat;
454 SCOTCH_archInit(&archdat);
455
456 SCOTCH_Num velosum = 0;
457 SCOTCH_dgraphSize (gr, &velosum, NULL, NULL, NULL);
458 SCOTCH_Num *goalsizes = new SCOTCH_Num[partnbr];
459 // TODO: The goalsizes are set as in Zoltan; not sure it is correct there
460 // or here.
461 // It appears velosum is global NUMBER of vertices, not global total
462 // vertex weight. I think we should use the latter.
463 // Fix this when we add vertex weights.
464 for (SCOTCH_Num i = 0; i < partnbr; i++)
465 goalsizes[i] = SCOTCH_Num(ceil(velosum * partsizes[i]));
466 delete [] partsizes;
467
468 SCOTCH_archCmpltw(&archdat, partnbr, goalsizes);
469
470 // Call partitioning; result returned in partloctab.
471 ierr = SCOTCH_dgraphMap(gr, &archdat, &stratstr, partloctab);
472
473 env->globalInputAssertion(__FILE__, __LINE__, "SCOTCH_dgraphMap",
474 !ierr, BASIC_ASSERTION, problemComm);
475
476 SCOTCH_archExit(&archdat);
477 delete [] goalsizes;
478
479 // TODO - metrics
480
481#ifdef SHOW_ZOLTAN2_SCOTCH_MEMORY
482 int me = env->comm_->getRank();
483#endif
484
485#ifdef HAVE_SCOTCH_ZOLTAN2_GETMEMORYMAX
486 if (me == 0){
487 size_t scotchBytes = SCOTCH_getMemoryMax();
488 std::cout << "Rank " << me << ": Maximum bytes used by Scotch: ";
489 std::cout << scotchBytes << std::endl;
490 }
491#endif
492
493 // Clean up PTScotch
494 SCOTCH_dgraphExit(gr);
495 free(gr);
496 SCOTCH_stratExit(&stratstr);
497
498 // Load answer into the solution.
499
500 ArrayRCP<part_t> partList;
501 if (nVtx > 0)
502 TPL_Traits<part_t, SCOTCH_Num>::SAVE_ARRAYRCP(&partList, partloctab, nVtx);
504
505 solution->setParts(partList);
506
507 env->memory("Zoltan2-Scotch: After creating solution");
508
509 // Clean up copies made due to differing data sizes.
512
513 if (nVwgts) delete [] velotab;
514 if (nEwgts) delete [] edlotab;
515
516#else // DO NOT HAVE MPI
517
518 // TODO: Handle serial case with calls to Scotch.
519 // TODO: For now, assign everything to rank 0 and assume only one part.
520 // TODO: Can probably use the code above for loading solution,
521 // TODO: instead of duplicating it here.
522 // TODO
523 // TODO: Actual logic should call Scotch when number of processes == 1.
524 ArrayView<const gno_t> vtxID;
525 ArrayView<StridedData<lno_t, scalar_t> > vwgts;
526 size_t nVtx = model->getVertexList(vtxID, vwgts);
527
528 ArrayRCP<part_t> partList(new part_t[nVtx], 0, nVtx, true);
529 for (size_t i = 0; i < nVtx; i++) partList[i] = 0;
530
531 solution->setParts(partList);
532
533#endif // DO NOT HAVE MPI
534}
535
537// Scale and round scalar_t weights (typically float or double) to
538// SCOTCH_Num (typically int or long).
539// subject to sum(weights) <= max_wgt_sum.
540// Only scale if deemed necessary.
541//
542// Note that we use ceil() instead of round() to avoid
543// rounding to zero weights.
544// Based on Zoltan's scale_round_weights, mode 1.
545
546template <typename Adapter>
547void AlgPTScotch<Adapter>::scale_weights(
548 size_t n,
549 StridedData<typename Adapter::lno_t, typename Adapter::scalar_t> &fwgts,
550 SCOTCH_Num *iwgts
551)
552{
553 const double INT_EPSILON = 1e-5;
554
555 SCOTCH_Num nonint, nonint_local = 0;
556 double sum_wgt, sum_wgt_local = 0.;
557 double max_wgt, max_wgt_local = 0.;
558
559 // Compute local sums of the weights
560 // Check whether all weights are integers
561 for (size_t i = 0; i < n; i++) {
562 double fw = double(fwgts[i]);
563 if (!nonint_local){
564 SCOTCH_Num tmp = (SCOTCH_Num) floor(fw + .5); /* Nearest int */
565 if (fabs((double)tmp-fw) > INT_EPSILON) {
566 nonint_local = 1;
567 }
568 }
569 sum_wgt_local += fw;
570 if (fw > max_wgt_local) max_wgt_local = fw;
571 }
572
573 Teuchos::reduceAll<int,SCOTCH_Num>(*problemComm, Teuchos::REDUCE_MAX, 1,
574 &nonint_local, &nonint);
575 Teuchos::reduceAll<int,double>(*problemComm, Teuchos::REDUCE_SUM, 1,
576 &sum_wgt_local, &sum_wgt);
577 Teuchos::reduceAll<int,double>(*problemComm, Teuchos::REDUCE_MAX, 1,
578 &max_wgt_local, &max_wgt);
579
580 double scale = 1.;
581 const double max_wgt_sum = double(SCOTCH_NUMMAX/8);
582
583 // Scaling needed if weights are not integers or weights'
584 // range is not sufficient
585 if (nonint || (max_wgt <= INT_EPSILON) || (sum_wgt > max_wgt_sum)) {
586 /* Calculate scale factor */
587 if (sum_wgt != 0.) scale = max_wgt_sum/sum_wgt;
588 }
589
590 /* Convert weights to positive integers using the computed scale factor */
591 for (size_t i = 0; i < n; i++)
592 iwgts[i] = (SCOTCH_Num) ceil(double(fwgts[i])*scale);
593
594}
595
596template <typename Adapter>
597int AlgPTScotch<Adapter>::setStrategy(SCOTCH_Strat * c_strat_ptr, const ParameterList &pList, int procRank) {
598 // get ordering parameters from parameter list
599 bool bPrintOutput = false; // will be true if rank 0 and verbose is true
600 const Teuchos::ParameterEntry *pe;
601
602 if (procRank == 0) {
603 pe = pList.getEntryPtr("scotch_verbose");
604 if (pe) {
605 bPrintOutput = pe->getValue(&bPrintOutput);
606 }
607 }
608
609 // get parameter setting for ordering default true/false
610 bool scotch_ordering_default = true;
611 pe = pList.getEntryPtr("scotch_ordering_default");
612 if (pe) {
613 scotch_ordering_default = pe->getValue(&scotch_ordering_default);
614 }
615 if (bPrintOutput) {
616 std::cout <<
617 "Scotch: Ordering default setting (parameter: scotch_ordering_default): "
618 << (scotch_ordering_default ? "True" : "False" ) << std::endl;
619 }
620
621 // set up error code for return
622 int ierr = 1; // will be set 0 if successful
623
624 if (scotch_ordering_default) {
625 // get parameter scotch_level
626 int scotch_level = 0;
627 pe = pList.getEntryPtr("scotch_level");
628 if (pe) {
629 scotch_level = pe->getValue(&scotch_level);
630 }
631 if (bPrintOutput) {
632 std::cout << "Scotch: Ordering level (parameter: scotch_level): " <<
633 scotch_level << std::endl;
634 }
635
636 // get parameter scotch_imbalance_ratio
637 double scotch_imbalance_ratio = 0.2;
638 pe = pList.getEntryPtr("scotch_imbalance_ratio");
639 if (pe) {
640 scotch_imbalance_ratio = pe->getValue(&scotch_imbalance_ratio);
641 }
642 if (bPrintOutput) {
643 std::cout << "Scotch: Ordering dissection imbalance ratio "
644 "(parameter: scotch_imbalance_ratio): "
645 << scotch_imbalance_ratio << std::endl;
646 }
647
648 SCOTCH_stratInit(c_strat_ptr);
649
650 ierr = SCOTCH_stratGraphOrderBuild( c_strat_ptr,
651 SCOTCH_STRATLEVELMAX | SCOTCH_STRATLEVELMIN |
652 SCOTCH_STRATLEAFSIMPLE | SCOTCH_STRATSEPASIMPLE,
653 scotch_level, scotch_imbalance_ratio); // based on Siva's example
654 }
655 else {
656 // get parameter scotch_ordering_strategy
657 std::string scotch_ordering_strategy_string("");
658 pe = pList.getEntryPtr("scotch_ordering_strategy");
659 if (pe) {
660 scotch_ordering_strategy_string =
661 pe->getValue(&scotch_ordering_strategy_string);
662 }
663 if (bPrintOutput) {
664 std::cout << "Scotch ordering strategy"
665 " (parameter: scotch_ordering_strategy): " <<
666 scotch_ordering_strategy_string << std::endl;
667 }
668
669 SCOTCH_stratInit(c_strat_ptr);
670 ierr = SCOTCH_stratGraphOrder( c_strat_ptr,
671 scotch_ordering_strategy_string.c_str());
672 }
673
674 return ierr;
675}
676
677template <typename Adapter>
679 const RCP<GlobalOrderingSolution<gno_t> > &solution) {
680 throw std::logic_error("AlgPTScotch does not yet support global ordering.");
681}
682
683template <typename Adapter>
685 const RCP<LocalOrderingSolution<lno_t> > &solution) {
686 // TODO: translate teuchos sublist parameters to scotch strategy string
687 // TODO: construct local graph model
688 // TODO: solve with scotch
689 // TODO: set solution fields
690
691 HELLO; // say hi so that we know we have called this method
692 int me = problemComm->getRank();
693 const ParameterList &pl = env->getParameters();
694 const Teuchos::ParameterEntry *pe;
695
696 bool isVerbose = false;
697 pe = pl.getEntryPtr("scotch_verbose");
698 if (pe)
699 isVerbose = pe->getValue(&isVerbose);
700
701 // build a local graph model
702 this->buildModel(true);
703 if (isVerbose && me==0) {
704 std::cout << "Built local graph model." << std::endl;
705 }
706
707 // based off of Siva's example
708 SCOTCH_Strat c_strat_ptr; // strategy
709 SCOTCH_Graph c_graph_ptr; // pointer to scotch graph
710 int ierr;
711
712 ierr = SCOTCH_graphInit( &c_graph_ptr);
713 if ( ierr != 0) {
714 throw std::runtime_error("Failed to initialize Scotch graph!");
715 } else if (isVerbose && me == 0) {
716 std::cout << "Initialized Scotch graph." << std::endl;
717 }
718
719 // Get vertex info
720 ArrayView<const gno_t> vtxID;
721 ArrayView<StridedData<lno_t, scalar_t> > vwgts;
722 size_t nVtx = model->getVertexList(vtxID, vwgts);
723 SCOTCH_Num vertnbr=0;
725
726 // Get edge info
727 ArrayView<const gno_t> edgeIds;
728 ArrayView<const offset_t> offsets;
729 ArrayView<StridedData<lno_t, scalar_t> > ewgts;
730
731 size_t nEdge = model->getEdgeList(edgeIds, offsets, ewgts);
732 SCOTCH_Num edgenbr=0;
734
735 SCOTCH_Num *verttab; // starting adj/vtx
737
738 SCOTCH_Num *edgetab; // adjacencies
740
741 // We don't use these arrays, but we need them as arguments to Scotch.
742 SCOTCH_Num *vendtab = NULL; // Assume consecutive storage for adj
743 //SCOTCH_Num *vendtab = verttab+1; // Assume consecutive storage for adj
744 // Get weight info.
745 SCOTCH_Num *velotab = NULL; // Vertex weights
746 SCOTCH_Num *vlbltab = NULL; // vertes labels
747 SCOTCH_Num *edlotab = NULL; // Edge weights
748
749 int nVwgts = model->getNumWeightsPerVertex();
750 int nEwgts = model->getNumWeightsPerEdge();
751 if (nVwgts > 1 && me == 0) {
752 std::cerr << "Warning: NumWeightsPerVertex is " << nVwgts
753 << " but Scotch allows only one weight. "
754 << " Zoltan2 will use only the first weight per vertex."
755 << std::endl;
756 }
757 if (nEwgts > 1 && me == 0) {
758 std::cerr << "Warning: NumWeightsPerEdge is " << nEwgts
759 << " but Scotch allows only one weight. "
760 << " Zoltan2 will use only the first weight per edge."
761 << std::endl;
762 }
763
764 if (nVwgts) {
765 velotab = new SCOTCH_Num[nVtx+1]; // +1 since Scotch wants all procs
766 // to have non-NULL arrays
767 scale_weights(nVtx, vwgts[0], velotab);
768 }
769
770 if (nEwgts) {
771 edlotab = new SCOTCH_Num[nEdge+1]; // +1 since Scotch wants all procs
772 // to have non-NULL arrays
773 scale_weights(nEdge, ewgts[0], edlotab);
774 }
775
776 // build scotch graph
777 int baseval = 0;
778 ierr = 1;
779 ierr = SCOTCH_graphBuild( &c_graph_ptr, baseval,
780 vertnbr, verttab, vendtab, velotab, vlbltab,
781 edgenbr, edgetab, edlotab);
782 if (ierr != 0) {
783 throw std::runtime_error("Failed to build Scotch graph!");
784 } else if (isVerbose && me == 0) {
785 std::cout << "Built Scotch graph." << std::endl;
786 }
787
788 ierr = SCOTCH_graphCheck(&c_graph_ptr);
789 if (ierr != 0) {
790 throw std::runtime_error("Graph is inconsistent!");
791 } else if (isVerbose && me == 0) {
792 std::cout << "Graph is consistent." << std::endl;
793 }
794
795 // set the strategy
796 ierr = AlgPTScotch<Adapter>::setStrategy(&c_strat_ptr, pl, me);
797
798 if (ierr != 0) {
799 throw std::runtime_error("Can't build strategy!");
800 }else if (isVerbose && me == 0) {
801 std::cout << "Graphing strategy built." << std::endl;
802 }
803
804 // Allocate results
805 SCOTCH_Num cblk = 0;
806 SCOTCH_Num *permtab; // permutation array
807 SCOTCH_Num *peritab; // inverse permutation array
808 SCOTCH_Num *rangetab; // separator range array
809 SCOTCH_Num *treetab; // separator tree
810
812 permtab = reinterpret_cast<SCOTCH_Num*>(solution->getPermutationView(false));
813 peritab = reinterpret_cast<SCOTCH_Num*>(solution->getPermutationView(true));
814 rangetab = reinterpret_cast<SCOTCH_Num*>(solution->getSeparatorRangeView());
815 treetab = reinterpret_cast<SCOTCH_Num*>(solution->getSeparatorTreeView());
816 }
817 else {
818 permtab = new SCOTCH_Num[nVtx];
819 peritab = new SCOTCH_Num[nVtx];
820 rangetab = new SCOTCH_Num[nVtx+1];
821 treetab = new SCOTCH_Num[nVtx];
822 }
823
824 ierr = SCOTCH_graphOrder(&c_graph_ptr, &c_strat_ptr, permtab, peritab,
825 &cblk, rangetab, treetab);
826 if (ierr != 0) {
827 throw std::runtime_error("Could not compute ordering!!");
828 } else if(isVerbose && me == 0) {
829 std::cout << "Ordering computed." << std::endl;
830 }
831
832 lno_t nSepBlocks;
834 solution->setNumSeparatorBlocks(nSepBlocks);
835
837
838 const ArrayRCP<lno_t> arv_perm = solution->getPermutationRCP(false);
839 for (size_t i = 0; i < nVtx; i++)
840 TPL_Traits<lno_t, SCOTCH_Num>::ASSIGN(arv_perm[i], permtab[i]);
841 delete [] permtab;
842
843 const ArrayRCP<lno_t> arv_peri = solution->getPermutationRCP(true);
844 for (size_t i = 0; i < nVtx; i++)
845 TPL_Traits<lno_t, SCOTCH_Num>::ASSIGN(arv_peri[i], peritab[i]);
846 delete [] peritab;
847
848 const ArrayRCP<lno_t> arv_range = solution->getSeparatorRangeRCP();
849 for (lno_t i = 0; i <= nSepBlocks; i++)
850 TPL_Traits<lno_t, SCOTCH_Num>::ASSIGN(arv_range[i], rangetab[i]);
851 delete [] rangetab;
852
853 const ArrayRCP<lno_t> arv_tree = solution->getSeparatorTreeRCP();
854 for (lno_t i = 0; i < nSepBlocks; i++)
855 TPL_Traits<lno_t, SCOTCH_Num>::ASSIGN(arv_tree[i], treetab[i]);
856 delete [] treetab;
857 }
858
859 solution->setHaveSeparator(true);
860
861 // reclaim memory
862 // Clean up copies made due to differing data sizes.
865
866 if (nVwgts) delete [] velotab;
867 if (nEwgts) delete [] edlotab;
868
869 SCOTCH_stratExit(&c_strat_ptr);
870 SCOTCH_graphFree(&c_graph_ptr);
871
872 if (isVerbose && problemComm->getRank() == 0) {
873 std::cout << "Freed Scotch graph!" << std::endl;
874 }
875 return 0;
876}
877
878} // namespace Zoltan2
879
880#endif // HAVE_ZOLTAN2_SCOTCH
881
883
884#endif
Zoltan2::BasicUserTypes< zscalar_t, zlno_t, zgno_t > user_t
Definition Metric.cpp:39
Defines the GraphModel interface.
Defines the OrderingSolution class.
Defines the PartitioningSolution class.
#define HELLO
Traits class to handle conversions between gno_t/lno_t and TPL data types (e.g., ParMETIS's idx_t,...
A gathering of useful namespace methods.
AlgPTScotch(const RCP< const Environment > &, const RCP< const Comm< int > > &, const RCP< const base_adapter_t > &)
Adapter::base_adapter_t base_adapter_t
static void getValidParameters(ParameterList &pl)
Set up validators specific to this algorithm.
Algorithm defines the base class for all algorithms.
virtual int globalOrder(const RCP< GlobalOrderingSolution< gno_t > > &)
Ordering method.
virtual int localOrder(const RCP< LocalOrderingSolution< lno_t > > &)
Ordering method.
virtual void partition(const RCP< PartitioningSolution< Adapter > > &)
Partitioning method.
Adapter::scalar_t scalar_t
static RCP< Teuchos::BoolParameterEntryValidator > getBoolValidator()
Exists to make setting up validators less cluttered.
static RCP< Teuchos::AnyNumberParameterEntryValidator > getAnyDoubleValidator()
Exists to make setting up validators less cluttered.
map_t::local_ordinal_type lno_t
Created by mbenlioglu on Aug 31, 2020.
std::bitset< NUM_MODEL_FLAGS > modelFlag_t
@ DETAILED_STATUS
sub-steps, each method's entry and exit
static void getScotchParameters(Teuchos::ParameterList &pl)
@ BASIC_ASSERTION
fast typical checks for valid arguments
@ VERTICES_ARE_MATRIX_ROWS
use matrix rows as graph vertices
@ VERTICES_ARE_MESH_ELEMENTS
use mesh elements as vertices
@ BUILD_SUBSET_GRAPH
ignore invalid neighbors
@ GENERATE_CONSECUTIVE_IDS
algorithm requires consecutive ids
@ SYMMETRIZE_INPUT_TRANSPOSE
model must symmetrize input
@ SYMMETRIZE_INPUT_BIPARTITE
model must symmetrize input
@ REMOVE_SELF_EDGES
algorithm requires no self edges
@ VERTICES_ARE_MESH_NODES
use mesh nodes as vertices
@ VERTICES_ARE_MATRIX_COLUMNS
use columns as graph vertices
@ VERTICES_ARE_MATRIX_NONZEROS
use nonzeros as graph vertices
@ BUILD_LOCAL_GRAPH
model represents graph within only one rank
SparseMatrixAdapter_t::part_t part_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)