MueLu Version of the Day
Loading...
Searching...
No Matches
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
61namespace 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.
105 RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
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()->getLocalNumElements() / 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->getLocalNumElements();
240 ArrayView<const GO> mapGIDs = map->getLocalElementList();
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
302 ArrayRCP<ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > CoordsData(dim);
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
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
Exception throws to report incompatible objects (like maps).
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
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)
static int GetLocalNumberOfRows(void *data, int *ierr)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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 void GetProblemGeometry(void *data, int numGIDEntries, int numLIDEntries, int numObjectIDs, ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids, int dim, double *coordinates, int *ierr)
static int GetProblemDimension(void *data, int *ierr)
Namespace for MueLu classes and methods.
@ Statistics1
Print more statistics.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.