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 
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.
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 
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.
LocalOrdinal LO
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.
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
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)