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