MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_ZoltanInterface_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_ZOLTANINTERFACE_DEF_HPP
47 #define MUELU_ZOLTANINTERFACE_DEF_HPP
48 
50 #if defined(HAVE_MUELU_ZOLTAN) && defined(HAVE_MPI)
51 
52 #include <Teuchos_Utils.hpp>
53 #include <Teuchos_DefaultMpiComm.hpp> //TODO: fwd decl.
54 #include <Teuchos_OpaqueWrapper.hpp> //TODO: fwd decl.
55 
56 #include "MueLu_Level.hpp"
57 #include "MueLu_Exceptions.hpp"
58 #include "MueLu_Monitor.hpp"
59 #include "MueLu_Utilities.hpp"
60 
61 namespace MueLu {
62 
63  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65  RCP<ParameterList> validParamList = rcp(new ParameterList());
66 
67  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Factory of the matrix A");
68  validParamList->set< RCP<const FactoryBase> >("number of partitions", Teuchos::null, "Instance of RepartitionHeuristicFactory.");
69  validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Factory of the coordinates");
70 
71  return validParamList;
72  }
73 
74 
75  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
77  Input(currentLevel, "A");
78  Input(currentLevel, "number of partitions");
79  Input(currentLevel, "Coordinates");
80  }
81 
82  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84  FactoryMonitor m(*this, "Build", level);
85 
86  RCP<Matrix> A = Get< RCP<Matrix> > (level, "A");
87  RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A);
88  RCP<const Map> rowMap;
89  if(bA != Teuchos::null) {
90  // Extracting the full the row map here...
91  RCP<const Map> bArowMap = bA->getRowMap();
92  RCP<const BlockedMap> bRowMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(bArowMap);
93  rowMap = bRowMap->getFullMap();
94  } else {
95  rowMap = A->getRowMap();
96  }
97 
98  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LocalOrdinal, GlobalOrdinal, Node> double_multivector_type;
99  RCP<double_multivector_type> Coords = Get< RCP<double_multivector_type> >(level, "Coordinates");
100  size_t dim = Coords->getNumVectors();
101  int numParts = Get<int>(level, "number of partitions");
102 
103  if (numParts == 1 || numParts == -1) {
104  // Running on one processor, so decomposition is the trivial one, all zeros.
106  Set(level, "Partition", decomposition);
107  return;
108  } else if (numParts == -1) {
109  // No repartitioning
110  RCP<Xpetra::Vector<GO,LO,GO,NO> > decomposition = Teuchos::null;
111  Set(level, "Partition", decomposition);
112  return;
113  }
114 
115  float zoltanVersion_;
116  Zoltan_Initialize(0, NULL, &zoltanVersion_);
117 
118  RCP<const Teuchos::MpiComm<int> > dupMpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(rowMap->getComm()->duplicate());
119  RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > zoltanComm = dupMpiComm->getRawMpiComm();
120 
121  RCP<Zoltan> zoltanObj_ = rcp(new Zoltan((*zoltanComm)())); //extract the underlying MPI_Comm handle and create a Zoltan object
122  if (zoltanObj_ == Teuchos::null)
123  throw Exceptions::RuntimeError("MueLu::Zoltan : Unable to create Zoltan data structure");
124 
125  // Tell Zoltan what kind of local/global IDs we will use.
126  // In our case, each GID is two ints and there are no local ids.
127  // One can skip this step if the IDs are just single ints.
128  int rv;
129  if ((rv = zoltanObj_->Set_Param("num_gid_entries", "1")) != ZOLTAN_OK)
130  throw Exceptions::RuntimeError("MueLu::Zoltan::Setup : setting parameter 'num_gid_entries' returned error code " + Teuchos::toString(rv));
131  if ((rv = zoltanObj_->Set_Param("num_lid_entries", "0") ) != ZOLTAN_OK)
132  throw Exceptions::RuntimeError("MueLu::Zoltan::Setup : setting parameter 'num_lid_entries' returned error code " + Teuchos::toString(rv));
133  if ((rv = zoltanObj_->Set_Param("obj_weight_dim", "1") ) != ZOLTAN_OK)
134  throw Exceptions::RuntimeError("MueLu::Zoltan::Setup : setting parameter 'obj_weight_dim' returned error code " + Teuchos::toString(rv));
135 
136  if (GetVerbLevel() & Statistics1) zoltanObj_->Set_Param("debug_level", "1");
137  else zoltanObj_->Set_Param("debug_level", "0");
138 
139  zoltanObj_->Set_Param("num_global_partitions", toString(numParts));
140 
141  zoltanObj_->Set_Num_Obj_Fn(GetLocalNumberOfRows, (void *) A.getRawPtr());
142  zoltanObj_->Set_Obj_List_Fn(GetLocalNumberOfNonzeros, (void *) A.getRawPtr());
143  zoltanObj_->Set_Num_Geom_Fn(GetProblemDimension, (void *) &dim);
144  zoltanObj_->Set_Geom_Multi_Fn(GetProblemGeometry, (void *) Coords.get());
145 
146  // Data pointers that Zoltan requires.
147  ZOLTAN_ID_PTR import_gids = NULL; // Global nums of objs to be imported
148  ZOLTAN_ID_PTR import_lids = NULL; // Local indices to objs to be imported
149  int *import_procs = NULL; // Proc IDs of procs owning objs to be imported.
150  int *import_to_part = NULL; // Partition #s to which imported objs should be assigned.
151  ZOLTAN_ID_PTR export_gids = NULL; // Global nums of objs to be exported
152  ZOLTAN_ID_PTR export_lids = NULL; // local indices to objs to be exported
153  int *export_procs = NULL; // Proc IDs of destination procs for objs to be exported.
154  int *export_to_part = NULL; // Partition #s for objs to be exported.
155  int num_imported; // Number of objs to be imported.
156  int num_exported; // Number of objs to be exported.
157  int newDecomp; // Flag indicating whether the decomposition has changed
158  int num_gid_entries; // Number of array entries in a global ID.
159  int num_lid_entries;
160 
161  {
162  SubFactoryMonitor m1(*this, "Zoltan RCB", level);
163  rv = zoltanObj_->LB_Partition(newDecomp, num_gid_entries, num_lid_entries,
164  num_imported, import_gids, import_lids, import_procs, import_to_part,
165  num_exported, export_gids, export_lids, export_procs, export_to_part);
166  if (rv == ZOLTAN_FATAL)
167  throw Exceptions::RuntimeError("Zoltan::LB_Partition() returned error code");
168  }
169 
170  // TODO check that A's row map is 1-1. Zoltan requires this.
171 
172  RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition;
173  if (newDecomp) {
174  decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, false); // Don't initialize, will be overwritten
175  ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
176 
177  int mypid = rowMap->getComm()->getRank();
178  for (typename ArrayRCP<GO>::iterator i = decompEntries.begin(); i != decompEntries.end(); ++i)
179  *i = mypid;
180 
181  LO blockSize = A->GetFixedBlockSize();
182  for (int i = 0; i < num_exported; ++i) {
183  // We have assigned Zoltan gids to first row GID in the block
184  // NOTE: Zoltan GIDs are different from GIDs in the Coordinates vector
185  LO localEl = rowMap->getLocalElement(export_gids[i]);
186  int partNum = export_to_part[i];
187  for (LO j = 0; j < blockSize; ++j)
188  decompEntries[localEl + j] = partNum;
189  }
190  }
191 
192  Set(level, "Partition", decomposition);
193 
194  zoltanObj_->LB_Free_Part(&import_gids, &import_lids, &import_procs, &import_to_part);
195  zoltanObj_->LB_Free_Part(&export_gids, &export_lids, &export_procs, &export_to_part);
196 
197  } //Build()
198 
199  //-------------------------------------------------------------------------------------------------------------
200  // GetLocalNumberOfRows
201  //-------------------------------------------------------------------------------------------------------------
202 
203  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
205  if (data == NULL) {
206  *ierr = ZOLTAN_FATAL;
207  return -1;
208  }
209  Matrix *A = (Matrix*) data;
210  *ierr = ZOLTAN_OK;
211 
212  LO blockSize = A->GetFixedBlockSize();
213  TEUCHOS_TEST_FOR_EXCEPTION(blockSize == 0, Exceptions::RuntimeError, "MueLu::Zoltan : Matrix has block size 0.");
214 
215  return A->getRowMap()->getNodeNumElements() / blockSize;
216  } //GetLocalNumberOfRows()
217 
218  //-------------------------------------------------------------------------------------------------------------
219  // GetLocalNumberOfNonzeros
220  //-------------------------------------------------------------------------------------------------------------
221 
222  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
224  GetLocalNumberOfNonzeros(void *data, int NumGidEntries, int /* NumLidEntries */, ZOLTAN_ID_PTR gids,
225  ZOLTAN_ID_PTR /* lids */, int /* wgtDim */, float *weights, int *ierr) {
226  if (data == NULL || NumGidEntries < 1) {
227  *ierr = ZOLTAN_FATAL;
228  return;
229  } else {
230  *ierr = ZOLTAN_OK;
231  }
232 
233  Matrix *A = (Matrix*) data;
234  RCP<const Map> map = A->getRowMap();
235 
236  LO blockSize = A->GetFixedBlockSize();
237  TEUCHOS_TEST_FOR_EXCEPTION(blockSize == 0, Exceptions::RuntimeError, "MueLu::Zoltan : Matrix has block size 0.");
238 
239  size_t numElements = map->getNodeNumElements();
240  ArrayView<const GO> mapGIDs = map->getNodeElementList();
241 
242  if (blockSize == 1) {
243  for (size_t i = 0; i < numElements; i++) {
244  gids[i] = as<ZOLTAN_ID_TYPE>(mapGIDs[i]);
245  weights[i] = A->getNumEntriesInLocalRow(i);
246  }
247 
248  } else {
249  LO numBlockElements = numElements / blockSize;
250 
251  for (LO i = 0; i < numBlockElements; i++) {
252  // Assign zoltan GID to the first row GID in the block
253  // NOTE: Zoltan GIDs are different from GIDs in the Coordinates vector
254  gids[i] = as<ZOLTAN_ID_TYPE>(mapGIDs[i*blockSize]);
255  weights[i] = 0.0;
256  for (LO j = 0; j < blockSize; j++)
257  weights[i] += A->getNumEntriesInLocalRow(i*blockSize+j);
258  }
259  }
260 
261  }
262 
263  //-------------------------------------------------------------------------------------------------------------
264  // GetProblemDimension
265  //-------------------------------------------------------------------------------------------------------------
266 
267  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
269  GetProblemDimension(void *data, int *ierr)
270  {
271  int dim = *((int*)data);
272  *ierr = ZOLTAN_OK;
273 
274  return dim;
275  } //GetProblemDimension
276 
277  //-------------------------------------------------------------------------------------------------------------
278  // GetProblemGeometry
279  //-------------------------------------------------------------------------------------------------------------
280 
281  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
283  GetProblemGeometry(void *data, int /* numGIDEntries */, int /* numLIDEntries */, int numObjectIDs,
284  ZOLTAN_ID_PTR /* gids */, ZOLTAN_ID_PTR /* lids */, int dim, double *coordinates, int *ierr)
285  {
286  if (data == NULL) {
287  *ierr = ZOLTAN_FATAL;
288  return;
289  }
290 
291  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> double_multivector_type;
292  double_multivector_type *Coords = (double_multivector_type*) data;
293 
294  if (dim != Teuchos::as<int>(Coords->getNumVectors())) {
295  //FIXME I'm assuming dim should be 1, 2, or 3 coming in?!
296  *ierr = ZOLTAN_FATAL;
297  return;
298  }
299 
300  TEUCHOS_TEST_FOR_EXCEPTION(numObjectIDs != Teuchos::as<int>(Coords->getLocalLength()), Exceptions::Incompatible, "Length of coordinates must be the same as the number of objects");
301 
303  for (int j = 0; j < dim; ++j)
304  CoordsData[j] = Coords->getData(j);
305 
306  size_t numElements = Coords->getLocalLength();
307  for (size_t i = 0; i < numElements; ++i)
308  for (int j = 0; j < dim; ++j)
309  coordinates[i*dim+j] = (double) CoordsData[j][i];
310 
311  *ierr = ZOLTAN_OK;
312 
313  } //GetProblemGeometry
314 
315 } //namespace MueLu
316 
317 #endif //if defined(HAVE_MUELU_ZOLTAN) && defined(HAVE_MPI)
318 
319 #endif // MUELU_ZOLTANINTERFACE_DEF_HPP
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
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)
Print more statistics.
LocalOrdinal LO
T * get() const
static int GetProblemDimension(void *data, int *ierr)
static void GetProblemGeometry(void *data, int numGIDEntries, int numLIDEntries, int numObjectIDs, ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids, int dim, double *coordinates, int *ierr)
Exception throws to report incompatible objects (like maps).
T * getRawPtr() const
void Build(Level &level) const
Build an object with this factory.
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.
static void GetLocalNumberOfNonzeros(void *data, int NumGidEntries, int NumLidEntries, ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids, int wgtDim, float *weights, int *ierr)
void DeclareInput(Level &level) const
Specifies the data that this class needs, and the factories that generate that data.
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
Exception throws to report errors in the internal logical of the program.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
iterator end() const
iterator begin() const
std::string toString(const T &t)
static int GetLocalNumberOfRows(void *data, int *ierr)