MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_Zoltan2Interface_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_ZOLTAN2INTERFACE_DEF_HPP
47#define MUELU_ZOLTAN2INTERFACE_DEF_HPP
48
49#include <sstream>
50#include <set>
51
53#if defined(HAVE_MUELU_ZOLTAN2) && defined(HAVE_MPI)
54
55#include <Zoltan2_XpetraMultiVectorAdapter.hpp>
56#include <Zoltan2_XpetraCrsGraphAdapter.hpp>
57#include <Zoltan2_PartitioningProblem.hpp>
58
59#include <Teuchos_Utils.hpp>
60#include <Teuchos_DefaultMpiComm.hpp>
61#include <Teuchos_OpaqueWrapper.hpp>
62
63#include "MueLu_Level.hpp"
64#include "MueLu_Exceptions.hpp"
65#include "MueLu_Monitor.hpp"
66#include "MueLu_Utilities.hpp"
67
68namespace MueLu {
69
70 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
72 defaultZoltan2Params = rcp(new ParameterList());
73 defaultZoltan2Params->set("algorithm", "multijagged");
74 defaultZoltan2Params->set("partitioning_approach", "partition");
75
76 // Improve scaling for communication bound algorithms by premigrating
77 // coordinates to a subset of processors.
78 // For more information, see Github issue #1538
79 defaultZoltan2Params->set("mj_premigration_option", 1);
80 }
81
82 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84 RCP<ParameterList> validParamList = rcp(new ParameterList());
85
86 validParamList->set< RCP<const FactoryBase> > ("A", Teuchos::null, "Factory of the matrix A");
87 validParamList->set< RCP<const FactoryBase> > ("number of partitions", Teuchos::null, "Instance of RepartitionHeuristicFactory.");
88 validParamList->set< RCP<const FactoryBase> > ("Coordinates", Teuchos::null, "Factory of the coordinates");
89 validParamList->set< RCP<const ParameterList> > ("ParameterList", Teuchos::null, "Zoltan2 parameters");
90 validParamList->set< RCP<const FactoryBase> > ("repartition: heuristic target rows per process", Teuchos::null, "Factory for number of rows per process to use with MultiJagged");
91
92 return validParamList;
93 }
94
95
96 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
98 Input(currentLevel, "A");
99 Input(currentLevel, "number of partitions");
100 const ParameterList& pL = GetParameterList();
101 // We do this dance, because we don't want "ParameterList" to be marked as used.
102 // Is there a better way?
103 Teuchos::ParameterEntry entry = pL.getEntry("ParameterList");
104 RCP<const Teuchos::ParameterList> providedList = Teuchos::any_cast<RCP<const Teuchos::ParameterList> >(entry.getAny(false));
105 if (providedList != Teuchos::null && providedList->isType<std::string>("algorithm")) {
106 const std::string algo = providedList->get<std::string>("algorithm");
107 if (algo == "multijagged") {
108 Input(currentLevel, "Coordinates");
109 Input(currentLevel, "repartition: heuristic target rows per process");
110 } else if (algo == "rcb") {
111 Input(currentLevel, "Coordinates");
112 }
113 } else {
114 Input(currentLevel, "repartition: heuristic target rows per process");
115 Input(currentLevel, "Coordinates");
116 }
117 }
118
119 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
121 FactoryMonitor m(*this, "Build", level);
122
123 typedef typename Teuchos::ScalarTraits<SC>::coordinateType real_type;
124 typedef typename Xpetra::MultiVector<real_type,LO,GO,NO> RealValuedMultiVector;
125
126 RCP<Matrix> A = Get<RCP<Matrix> >(level, "A");
127 RCP<const Map> rowMap = A->getRowMap();
128 LO blkSize = A->GetFixedBlockSize();
129
130 int numParts = Get<int>(level, "number of partitions");
131 if (numParts == 1 || numParts == -1) {
132 // Single processor, decomposition is trivial: all zeros
133 RCP<Xpetra::Vector<GO,LO,GO,NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
134 Set(level, "Partition", decomposition);
135 return;
136 }/* else if (numParts == -1) {
137 // No repartitioning
138 RCP<Xpetra::Vector<GO,LO,GO,NO> > decomposition = Teuchos::null; //Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
139 //decomposition->putScalar(Teuchos::as<Scalar>(rowMap->getComm()->getRank()));
140 Set(level, "Partition", decomposition);
141 return;
142 }*/
143
144 const ParameterList& pL = GetParameterList();
145
146 RCP<const ParameterList> providedList = pL.get<RCP<const ParameterList> >("ParameterList");
147 ParameterList Zoltan2Params;
148 if (providedList != Teuchos::null)
149 Zoltan2Params = *providedList;
150
151 // Merge defalt Zoltan2 parameters with user provided
152 // If default and user parameters contain the same parameter name, user one is always preferred
153 for (ParameterList::ConstIterator param = defaultZoltan2Params->begin(); param != defaultZoltan2Params->end(); param++) {
154 const std::string& pName = defaultZoltan2Params->name(param);
155 if (!Zoltan2Params.isParameter(pName))
156 Zoltan2Params.setEntry(pName, defaultZoltan2Params->getEntry(pName));
157 }
158 Zoltan2Params.set("num_global_parts", Teuchos::as<int>(numParts));
159
160 GetOStream(Runtime0) << "Zoltan2 parameters:\n----------\n" << Zoltan2Params << "----------" << std::endl;
161
162 const std::string& algo = Zoltan2Params.get<std::string>("algorithm");
163
164 if (algo == "multijagged" || algo == "rcb") {
165
166 RCP<RealValuedMultiVector> coords = Get<RCP<RealValuedMultiVector> >(level, "Coordinates");
167 RCP<const Map> map = coords->getMap();
168 GO numElements = map->getLocalNumElements();
169
170 // Check that the number of local coordinates is consistent with the #rows in A
171 TEUCHOS_TEST_FOR_EXCEPTION(rowMap->getLocalNumElements()/blkSize != coords->getLocalLength(), Exceptions::Incompatible,
172 "Coordinate vector length (" + toString(coords->getLocalLength()) << " is incompatible with number of block rows in A ("
173 + toString(rowMap->getLocalNumElements()/blkSize) + "The vector length should be the same as the number of mesh points.");
174#ifdef HAVE_MUELU_DEBUG
175 GO indexBase = rowMap->getIndexBase();
176 GetOStream(Runtime0) << "Checking consistence of row and coordinates maps" << std::endl;
177 // Make sure that logical blocks in row map coincide with logical nodes in coordinates map
178 ArrayView<const GO> rowElements = rowMap->getLocalElementList();
179 ArrayView<const GO> coordsElements = map ->getLocalElementList();
180 for (LO i = 0; i < Teuchos::as<LO>(numElements); i++)
181 TEUCHOS_TEST_FOR_EXCEPTION((coordsElements[i]-indexBase)*blkSize + indexBase != rowElements[i*blkSize],
182 Exceptions::RuntimeError, "i = " << i << ", coords GID = " << coordsElements[i]
183 << ", row GID = " << rowElements[i*blkSize] << ", blkSize = " << blkSize << std::endl);
184#endif
185
186 typedef Zoltan2::XpetraMultiVectorAdapter<RealValuedMultiVector> InputAdapterType;
187 typedef Zoltan2::PartitioningProblem<InputAdapterType> ProblemType;
188
189 Array<real_type> weightsPerRow(numElements);
190 for (LO i = 0; i < numElements; i++) {
191 weightsPerRow[i] = 0.0;
192
193 for (LO j = 0; j < blkSize; j++) {
194 weightsPerRow[i] += A->getNumEntriesInLocalRow(i*blkSize+j);
195 }
196 }
197
198 // MultiJagged: Grab the target rows per process from the Heuristic to use unless the Zoltan2 list says otherwise
199 if(algo == "multijagged" && !Zoltan2Params.isParameter("mj_premigration_coordinate_count")) {
200 LO heuristicTargetRowsPerProcess = Get<LO>(level,"repartition: heuristic target rows per process");
201 Zoltan2Params.set("mj_premigration_coordinate_count", heuristicTargetRowsPerProcess);
202 }
203 const bool writeZoltan2DebuggingFiles = Zoltan2Params.get("mj_debug",false);
204 Zoltan2Params.remove("mj_debug");
205
206 std::vector<int> strides;
207 std::vector<const real_type*> weights(1, weightsPerRow.getRawPtr());
208
209 RCP<const Teuchos::MpiComm<int> > dupMpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(rowMap->getComm()->duplicate());
210 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > zoltanComm = dupMpiComm->getRawMpiComm();
211
212 InputAdapterType adapter(coords, weights, strides);
213 RCP<ProblemType> problem(new ProblemType(&adapter, &Zoltan2Params, (*zoltanComm)()));
214
215 {
216 SubFactoryMonitor m1(*this, "Zoltan2 " + toString(algo), level);
217 if (writeZoltan2DebuggingFiles)
218 adapter.generateFiles(("mj_debug.lvl_"+std::to_string(level.GetLevelID())).c_str(), *(rowMap->getComm()));
219 problem->solve();
220 }
221
222 RCP<Xpetra::Vector<GO,LO,GO,NO> > decomposition = Xpetra::VectorFactory<GO,LO,GO,NO>::Build(rowMap, false);
223 ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
224
225 const typename InputAdapterType::part_t * parts = problem->getSolution().getPartListView();
226
227 for (GO i = 0; i < numElements; i++) {
228 int partNum = parts[i];
229
230 for (LO j = 0; j < blkSize; j++)
231 decompEntries[i*blkSize + j] = partNum;
232 }
233
234 Set(level, "Partition", decomposition);
235
236 } else {
237
238 GO numElements = rowMap->getLocalNumElements();
239
240 typedef Zoltan2::XpetraCrsGraphAdapter<CrsGraph> InputAdapterType;
241 typedef Zoltan2::PartitioningProblem<InputAdapterType> ProblemType;
242
243 RCP<const Teuchos::MpiComm<int> > dupMpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(rowMap->getComm()->duplicate());
244 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > zoltanComm = dupMpiComm->getRawMpiComm();
245
246 InputAdapterType adapter(A->getCrsGraph());
247 RCP<ProblemType> problem(new ProblemType(&adapter, &Zoltan2Params, (*zoltanComm)()));
248
249 {
250 SubFactoryMonitor m1(*this, "Zoltan2 " + toString(algo), level);
251 problem->solve();
252 }
253
254 RCP<Xpetra::Vector<GO,LO,GO,NO> > decomposition = Xpetra::VectorFactory<GO,LO,GO,NO>::Build(rowMap, false);
255 ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
256
257 const typename InputAdapterType::part_t * parts = problem->getSolution().getPartListView();
258
259 // For blkSize > 1, ignore solution for every row but the first ones in a block.
260 for (GO i = 0; i < numElements/blkSize; i++) {
261 int partNum = parts[i*blkSize];
262
263 for (LO j = 0; j < blkSize; j++)
264 decompEntries[i*blkSize + j] = partNum;
265 }
266
267 Set(level, "Partition", decomposition);
268 }
269 }
270
271} //namespace MueLu
272
273#endif //if defined(HAVE_MUELU_ZOLTAN2) && defined(HAVE_MPI)
274
275#endif // MUELU_ZOLTAN2INTERFACE_DEF_HPP
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.
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void Build(Level &currentLevel) const
Build an object with this factory.
Namespace for MueLu classes and methods.
@ Runtime0
One-liner description of what is happening.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.