MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_Zoltan2Interface_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
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
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_ZOLTAN2INTERFACE_DEF_HPP
47 #define MUELU_ZOLTAN2INTERFACE_DEF_HPP
48 
49 #include <sstream>
50 #include <set>
51 
53 #if defined(HAVE_MUELU_ZOLTAN2) && defined(HAVE_MPI)
54 
55 #include <Zoltan2_XpetraMultiVectorAdapter.hpp>
56 #include <Zoltan2_XpetraCrsGraphAdapter.hpp>
57 #include <Zoltan2_PartitioningProblem.hpp>
58 
59 #include <Teuchos_Utils.hpp>
61 #include <Teuchos_OpaqueWrapper.hpp>
62 
63 #include "MueLu_Level.hpp"
64 #include "MueLu_Exceptions.hpp"
65 #include "MueLu_Monitor.hpp"
66 #include "MueLu_Utilities.hpp"
67 
68 namespace MueLu {
69 
70  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
72  defaultZoltan2Params = rcp(new ParameterList());
73  defaultZoltan2Params->set("algorithm", "multijagged");
74  defaultZoltan2Params->set("partitioning_approach", "partition");
75 
76  // Improve scaling for communication bound algorithms by premigrating
77  // coordinates to a subset of processors.
78  // For more information, see Github issue #1538
79  defaultZoltan2Params->set("mj_premigration_option", 1);
80  }
81 
82  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84  RCP<ParameterList> validParamList = rcp(new ParameterList());
85 
86  validParamList->set< RCP<const FactoryBase> > ("A", Teuchos::null, "Factory of the matrix A");
87  validParamList->set< RCP<const FactoryBase> > ("number of partitions", Teuchos::null, "Instance of RepartitionHeuristicFactory.");
88  validParamList->set< RCP<const FactoryBase> > ("Coordinates", Teuchos::null, "Factory of the coordinates");
89  validParamList->set< RCP<const ParameterList> > ("ParameterList", Teuchos::null, "Zoltan2 parameters");
90  validParamList->set< RCP<const FactoryBase> > ("repartition: heuristic target rows per process", Teuchos::null, "Factory for number of rows per process to use with MultiJagged");
91 
92  return validParamList;
93  }
94 
95 
96  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
98  Input(currentLevel, "A");
99  Input(currentLevel, "number of partitions");
100  const ParameterList& pL = GetParameterList();
101  // We do this dance, because we don't want "ParameterList" to be marked as used.
102  // Is there a better way?
103  Teuchos::ParameterEntry entry = pL.getEntry("ParameterList");
104  RCP<const Teuchos::ParameterList> providedList = Teuchos::any_cast<RCP<const Teuchos::ParameterList> >(entry.getAny(false));
105  if (providedList != Teuchos::null && providedList->isType<std::string>("algorithm")) {
106  const std::string algo = providedList->get<std::string>("algorithm");
107  if (algo == "multijagged") {
108  Input(currentLevel, "Coordinates");
109  Input(currentLevel, "repartition: heuristic target rows per process");
110  } else if (algo == "rcb") {
111  Input(currentLevel, "Coordinates");
112  }
113  } else {
114  Input(currentLevel, "repartition: heuristic target rows per process");
115  Input(currentLevel, "Coordinates");
116  }
117  }
118 
119  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
121  FactoryMonitor m(*this, "Build", level);
122 
123  typedef typename Teuchos::ScalarTraits<SC>::coordinateType real_type;
124  typedef typename Xpetra::MultiVector<real_type,LO,GO,NO> RealValuedMultiVector;
125 
126  RCP<Matrix> A = Get<RCP<Matrix> >(level, "A");
127  RCP<const Map> rowMap = A->getRowMap();
128  LO blkSize = A->GetFixedBlockSize();
129 
130  int numParts = Get<int>(level, "number of partitions");
131  if (numParts == 1 || numParts == -1) {
132  // Single processor, decomposition is trivial: all zeros
134  Set(level, "Partition", decomposition);
135  return;
136  }/* else if (numParts == -1) {
137  // No repartitioning
138  RCP<Xpetra::Vector<GO,LO,GO,NO> > decomposition = Teuchos::null; //Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
139  //decomposition->putScalar(Teuchos::as<Scalar>(rowMap->getComm()->getRank()));
140  Set(level, "Partition", decomposition);
141  return;
142  }*/
143 
144  const ParameterList& pL = GetParameterList();
145 
146  RCP<const ParameterList> providedList = pL.get<RCP<const ParameterList> >("ParameterList");
147  ParameterList Zoltan2Params;
148  if (providedList != Teuchos::null)
149  Zoltan2Params = *providedList;
150 
151  // Merge defalt Zoltan2 parameters with user provided
152  // If default and user parameters contain the same parameter name, user one is always preferred
153  for (ParameterList::ConstIterator param = defaultZoltan2Params->begin(); param != defaultZoltan2Params->end(); param++) {
154  const std::string& pName = defaultZoltan2Params->name(param);
155  if (!Zoltan2Params.isParameter(pName))
156  Zoltan2Params.setEntry(pName, defaultZoltan2Params->getEntry(pName));
157  }
158  Zoltan2Params.set("num_global_parts", Teuchos::as<int>(numParts));
159 
160  GetOStream(Runtime0) << "Zoltan2 parameters:\n----------\n" << Zoltan2Params << "----------" << std::endl;
161 
162  const std::string& algo = Zoltan2Params.get<std::string>("algorithm");
163 
164  if (algo == "multijagged" || algo == "rcb") {
165 
166  RCP<RealValuedMultiVector> coords = Get<RCP<RealValuedMultiVector> >(level, "Coordinates");
167  RCP<const Map> map = coords->getMap();
168  GO numElements = map->getNodeNumElements();
169 
170  // Check that the number of local coordinates is consistent with the #rows in A
171  TEUCHOS_TEST_FOR_EXCEPTION(rowMap->getNodeNumElements()/blkSize != coords->getLocalLength(), Exceptions::Incompatible,
172  "Coordinate vector length (" + toString(coords->getLocalLength()) << " is incompatible with number of block rows in A ("
173  + toString(rowMap->getNodeNumElements()/blkSize) + "The vector length should be the same as the number of mesh points.");
174 #ifdef HAVE_MUELU_DEBUG
175  GO indexBase = rowMap->getIndexBase();
176  GetOStream(Runtime0) << "Checking consistence of row and coordinates maps" << std::endl;
177  // Make sure that logical blocks in row map coincide with logical nodes in coordinates map
178  ArrayView<const GO> rowElements = rowMap->getNodeElementList();
179  ArrayView<const GO> coordsElements = map ->getNodeElementList();
180  for (LO i = 0; i < Teuchos::as<LO>(numElements); i++)
181  TEUCHOS_TEST_FOR_EXCEPTION((coordsElements[i]-indexBase)*blkSize + indexBase != rowElements[i*blkSize],
182  Exceptions::RuntimeError, "i = " << i << ", coords GID = " << coordsElements[i]
183  << ", row GID = " << rowElements[i*blkSize] << ", blkSize = " << blkSize << std::endl);
184 #endif
185 
186  typedef Zoltan2::XpetraMultiVectorAdapter<RealValuedMultiVector> InputAdapterType;
187  typedef Zoltan2::PartitioningProblem<InputAdapterType> ProblemType;
188 
189  Array<real_type> weightsPerRow(numElements);
190  for (LO i = 0; i < numElements; i++) {
191  weightsPerRow[i] = 0.0;
192 
193  for (LO j = 0; j < blkSize; j++) {
194  weightsPerRow[i] += A->getNumEntriesInLocalRow(i*blkSize+j);
195  }
196  }
197 
198  // MultiJagged: Grab the target rows per process from the Heuristic to use unless the Zoltan2 list says otherwise
199  if(algo == "multijagged" && !Zoltan2Params.isParameter("mj_premigration_coordinate_count")) {
200  LO heuristicTargetRowsPerProcess = Get<LO>(level,"repartition: heuristic target rows per process");
201  Zoltan2Params.set("mj_premigration_coordinate_count", heuristicTargetRowsPerProcess);
202  }
203 
204  std::vector<int> strides;
205  std::vector<const real_type*> weights(1, weightsPerRow.getRawPtr());
206 
207  RCP<const Teuchos::MpiComm<int> > dupMpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(rowMap->getComm()->duplicate());
208  RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > zoltanComm = dupMpiComm->getRawMpiComm();
209 
210  InputAdapterType adapter(coords, weights, strides);
211  RCP<ProblemType> problem(new ProblemType(&adapter, &Zoltan2Params, (*zoltanComm)()));
212 
213  {
214  SubFactoryMonitor m1(*this, "Zoltan2 " + toString(algo), level);
215  problem->solve();
216  }
217 
219  ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
220 
221  const typename InputAdapterType::part_t * parts = problem->getSolution().getPartListView();
222 
223  for (GO i = 0; i < numElements; i++) {
224  int partNum = parts[i];
225 
226  for (LO j = 0; j < blkSize; j++)
227  decompEntries[i*blkSize + j] = partNum;
228  }
229 
230  Set(level, "Partition", decomposition);
231 
232  } else {
233 
234  GO numElements = rowMap->getNodeNumElements();
235 
236  typedef Zoltan2::XpetraCrsGraphAdapter<CrsGraph> InputAdapterType;
237  typedef Zoltan2::PartitioningProblem<InputAdapterType> ProblemType;
238 
239  RCP<const Teuchos::MpiComm<int> > dupMpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(rowMap->getComm()->duplicate());
240  RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > zoltanComm = dupMpiComm->getRawMpiComm();
241 
242  InputAdapterType adapter(A->getCrsGraph());
243  RCP<ProblemType> problem(new ProblemType(&adapter, &Zoltan2Params, (*zoltanComm)()));
244 
245  {
246  SubFactoryMonitor m1(*this, "Zoltan2 " + toString(algo), level);
247  problem->solve();
248  }
249 
251  ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
252 
253  const typename InputAdapterType::part_t * parts = problem->getSolution().getPartListView();
254 
255  for (GO i = 0; i < numElements; i++) {
256  int partNum = parts[i];
257 
258  for (LO j = 0; j < blkSize; j++)
259  decompEntries[i*blkSize + j] = partNum;
260  }
261 
262  Set(level, "Partition", decomposition);
263  }
264  }
265 
266 } //namespace MueLu
267 
268 #endif //if defined(HAVE_MUELU_ZOLTAN2) && defined(HAVE_MPI)
269 
270 #endif // MUELU_ZOLTAN2INTERFACE_DEF_HPP
const std::string & name() const
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
GlobalOrdinal GO
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
LocalOrdinal LO
One-liner description of what is happening.
T * get() const
void Build(Level &currentLevel) const
Build an object with this factory.
Exception throws to report incompatible objects (like maps).
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
params_t::ConstIterator ConstIterator
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
any & getAny(bool activeQry=true)
bool isType(const std::string &name) const
Exception throws to report errors in the internal logical of the program.
ParameterEntry & getEntry(const std::string &name)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.