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