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