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