Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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::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 
316 private:
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 
334 template <typename Adapter>
335 void 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 
372 template <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;
408  TPL_Traits<SCOTCH_Num, size_t>::ASSIGN(vertlocnbr, nVtx);
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;
419  TPL_Traits<SCOTCH_Num, size_t>::ASSIGN(edgelocnbr, nEdge);
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 
581 template <typename Adapter>
582 void 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 
631 template <typename Adapter>
632 int 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 
712 template <typename Adapter>
714  const RCP<GlobalOrderingSolution<gno_t> > &solution) {
715  throw std::logic_error("AlgPTScotch does not yet support global ordering.");
716 }
717 
718 template <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 
846  if (TPL_Traits<lno_t, SCOTCH_Num>::OK_TO_CAST()) {
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;
868  TPL_Traits<lno_t, SCOTCH_Num>::ASSIGN(nSepBlocks, cblk);
869  solution->setNumSeparatorBlocks(nSepBlocks);
870 
871  if (!TPL_Traits<lno_t, SCOTCH_Num>::OK_TO_CAST()) {
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
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
virtual void partition(const RCP< PartitioningSolution< Adapter > > &)
Partitioning method.
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)
virtual int globalOrder(const RCP< GlobalOrderingSolution< gno_t > > &)
Ordering method.
sub-steps, each method&#39;s entry and exit
SparseMatrixAdapter_t::part_t part_t
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
AlgPTScotch(const RCP< const Environment > &, const RCP< const Comm< int > > &, const RCP< const base_adapter_t > &)
Algorithm defines the base class for all algorithms.
static void ASSIGN(first_t &a, second_t b)
static void getValidParameters(ParameterList &pl)
Set up validators specific to this algorithm.
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)
Defines the GraphModel interface.
virtual int localOrder(const RCP< LocalOrderingSolution< lno_t > > &)
Ordering method.
A gathering of useful namespace methods.
static void DELETE_ARRAY(first_t **a)
model represents graph within only one rank
Zoltan2::BasicUserTypes< zscalar_t, zlno_t, zgno_t > user_t
Definition: Metric.cpp:74