MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_IsorropiaInterface_def.hpp
Go to the documentation of this file.
1/*
2 * MueLu_IsorropiaInterface_def.hpp
3 *
4 * Created on: Jun 10, 2013
5 * Author: tobias
6 */
7
8#ifndef MUELU_ISORROPIAINTERFACE_DEF_HPP_
9#define MUELU_ISORROPIAINTERFACE_DEF_HPP_
10
12
13#include <Teuchos_Utils.hpp>
14//#include <Teuchos_DefaultMpiComm.hpp> //TODO: fwd decl.
15//#include <Teuchos_OpaqueWrapper.hpp> //TODO: fwd decl.
16
17#include <Xpetra_MapFactory.hpp>
18#include <Xpetra_CrsGraphFactory.hpp>
19#include <Xpetra_BlockedMap.hpp>
20#include <Xpetra_BlockedCrsMatrix.hpp>
21
22#ifdef HAVE_MUELU_ISORROPIA
23#include <Isorropia_Exception.hpp>
24
25
26#ifdef HAVE_MUELU_EPETRA
27#include <Xpetra_EpetraCrsGraph.hpp>
28#include <Epetra_CrsGraph.h>
29#include <Isorropia_EpetraPartitioner.hpp>
30#endif
31
32#ifdef HAVE_MUELU_TPETRA
33#include <Xpetra_TpetraCrsGraph.hpp>
34#endif
35#endif // ENDIF HAVE_MUELU_ISORROPIA
36
37#include "MueLu_Level.hpp"
38#include "MueLu_Exceptions.hpp"
39#include "MueLu_Monitor.hpp"
40#include "MueLu_Graph.hpp"
41#include "MueLu_AmalgamationInfo.hpp"
42#include "MueLu_Utilities.hpp"
43
44namespace MueLu {
45
46 template <class LocalOrdinal, class GlobalOrdinal, class Node>
48 RCP<ParameterList> validParamList = rcp(new ParameterList());
49
50 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Factory of the matrix A");
51 validParamList->set< RCP<const FactoryBase> >("number of partitions", Teuchos::null, "Instance of RepartitionHeuristicFactory.");
52 validParamList->set< RCP<const FactoryBase> >("UnAmalgamationInfo", Teuchos::null, "Generating factory of UnAmalgamationInfo");
53
54 return validParamList;
55 }
56
57
58 template <class LocalOrdinal, class GlobalOrdinal, class Node>
60 Input(currentLevel, "A");
61 Input(currentLevel, "number of partitions");
62 Input(currentLevel, "UnAmalgamationInfo");
63 }
64
65 template <class LocalOrdinal, class GlobalOrdinal, class Node>
67 FactoryMonitor m(*this, "Build", level);
68 typedef Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node> BlockMap;
69
70 RCP<Matrix> A = Get< RCP<Matrix> >(level, "A");
71 RCP<AmalgamationInfo> amalInfo = Get< RCP<AmalgamationInfo> >(level, "UnAmalgamationInfo");
72 GO numParts = Get< int >(level, "number of partitions");
73
74 RCP<const Map> rowMap = A->getRowMap();
75 RCP<const Map> colMap = A->getColMap();
76
77 if (numParts == 1 || numParts == -1) {
78 // Running on one processor, so decomposition is the trivial one, all zeros.
79 RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
80 Set(level, "AmalgamatedPartition", decomposition);
81 return;
82 }
83
84 // ok, reconstruct graph information of matrix A
85 // Note, this is the non-rebalanced matrix A and we need the graph
86 // of the non-rebalanced matrix A. We cannot make use of the "Graph"
87 // that is/will be built for the aggregates later for several reasons
88 // 1) the "Graph" for the aggregates is meant to be based on the rebalanced matrix A
89 // 2) we cannot call a CoalesceDropFactory::Build here since this is very complicated and
90 // completely messes up the whole factory chain
91 // 3) CoalesceDropFactory is meant to provide some minimal Graph information for the aggregation
92 // (LWGraph), but here we need the full CrsGraph for Isorropia as input
93
94 // That is, why this code is somewhat redundant to CoalesceDropFactory
95
96 LO blockdim = 1; // block dim for fixed size blocks
97 GO indexBase = rowMap->getIndexBase(); // index base of maps
98 GO offset = 0;
99 LO blockid = -1; // block id in strided map
100 LO nStridedOffset = 0; // DOF offset for strided block id "blockid" (default = 0)
101 //LO stridedblocksize = blockdim; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
102
103 // 1) check for blocking/striding information
104 // fill above variables
105 if(A->IsView("stridedMaps") &&
106 Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
107 Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
108 RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
109 TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,Exceptions::BadCast,"MueLu::IsorropiaInterface::Build: cast to strided row map failed.");
110 blockdim = strMap->getFixedBlockSize();
111 offset = strMap->getOffset();
112 blockid = strMap->getStridedBlockId();
113 if (blockid > -1) {
114 std::vector<size_t> stridingInfo = strMap->getStridingData();
115 for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
116 nStridedOffset += stridingInfo[j];
117 //stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
118
119 }// else {
120 // stridedblocksize = blockdim;
121 //}
122 oldView = A->SwitchToView(oldView);
123 //GetOStream(Statistics0) << "IsorropiaInterface::Build():" << " found blockdim=" << blockdim << " from strided maps (blockid=" << blockid << ", strided block size=" << stridedblocksize << "). offset=" << offset << std::endl;
124 } else GetOStream(Statistics0) << "IsorropiaInterface::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
125
126 // 2) get row map for amalgamated matrix (graph of A)
127 // with same distribution over all procs as row map of A
128 RCP<const Map> nodeMap= amalInfo->getNodeRowMap();
129 RCP<const BlockedMap> bnodeMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(nodeMap);
130 if(!bnodeMap.is_null()) nodeMap=bnodeMap->getMap();
131
132 GetOStream(Statistics0) << "IsorropiaInterface:Build(): nodeMap " << nodeMap->getLocalNumElements() << "/" << nodeMap->getGlobalNumElements() << " elements" << std::endl;
133
134
135 // 3) create graph of amalgamated matrix
136 RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, A->getLocalMaxNumRowEntries()*blockdim);
137
138 // 4) do amalgamation. generate graph of amalgamated matrix
139 for(LO row=0; row<Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); row++) {
140 // get global DOF id
141 GO grid = rowMap->getGlobalElement(row);
142
143 // translate grid to nodeid
144 // JHU 2019-20-May this is identical to AmalgamationFactory::DOFGid2NodeId(), and is done
145 // to break a circular dependence when libraries are built statically
146 GO nodeId = (grid - offset - indexBase) / blockdim + indexBase;
147
148 size_t nnz = A->getNumEntriesInLocalRow(row);
149 Teuchos::ArrayView<const LO> indices;
150 Teuchos::ArrayView<const SC> vals;
151 A->getLocalRowView(row, indices, vals);
152
153 RCP<std::vector<GO> > cnodeIds = Teuchos::rcp(new std::vector<GO>); // global column block ids
154 LO realnnz = 0;
155 for(LO col=0; col<Teuchos::as<LO>(nnz); col++) {
156 GO gcid = colMap->getGlobalElement(indices[col]); // global column id
157
158 if(vals[col]!=0.0) {
159 // JHU 2019-20-May this is identical to AmalgamationFactory::DOFGid2NodeId(), and is done
160 // to break a circular dependence when libraries are built statically
161 GO cnodeId = (gcid - offset - indexBase) / blockdim + indexBase;
162 cnodeIds->push_back(cnodeId);
163 realnnz++; // increment number of nnz in matrix row
164 }
165 }
166
167 Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp( cnodeIds );
168
169 if(arr_cnodeIds.size() > 0 )
170 crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
171 }
172 // fill matrix graph
173 crsGraph->fillComplete(nodeMap,nodeMap);
174
175#ifdef HAVE_MPI
176#ifdef HAVE_MUELU_ISORROPIA
177
178 // prepare parameter list for Isorropia
179 Teuchos::ParameterList paramlist;
180 paramlist.set("NUM PARTS", toString(numParts));
181
182 /*STRUCTURALLY SYMMETRIC [NO/yes] (is symmetrization required?)
183 PARTITIONING METHOD [block/cyclic/random/rcb/rib/hsfc/graph/HYPERGRAPH]
184 NUM PARTS [int k] (global number of parts)
185 IMBALANCE TOL [float tol] (1.0 is perfect balance)
186 BALANCE OBJECTIVE [ROWS/nonzeros]
187 */
188 Teuchos::ParameterList& sublist = paramlist.sublist("Zoltan");
189 sublist.set("LB_APPROACH", "PARTITION");
190
191
192
193#ifdef HAVE_MUELU_EPETRA
194 RCP< Xpetra::EpetraCrsGraphT<GO, Node> > epCrsGraph = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsGraphT<GO, Node> >(crsGraph);
195 if(epCrsGraph != Teuchos::null) {
196 RCP< const Epetra_CrsGraph> epetraCrsGraph = epCrsGraph->getEpetra_CrsGraph();
197
198 RCP<Isorropia::Epetra::Partitioner> isoPart = Teuchos::rcp(new Isorropia::Epetra::Partitioner(epetraCrsGraph,paramlist));
199
200 int size = 0;
201 const int* array = NULL;
202 isoPart->extractPartsView(size,array);
203
204 TEUCHOS_TEST_FOR_EXCEPTION(size != Teuchos::as<int>(nodeMap->getLocalNumElements()), Exceptions::RuntimeError, "length of array returned from extractPartsView does not match local length of rowMap");
205
206 RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(nodeMap, false);
207 ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
208
209 // fill vector with amalgamated information about partitioning
210 for(int i = 0; i<size; i++) {
211 decompEntries[i] = Teuchos::as<GO>(array[i]);
212 }
213
214 Set(level, "AmalgamatedPartition", decomposition);
215
216 }
217#endif // ENDIF HAVE_MUELU_EPETRA
218
219#ifdef HAVE_MUELU_TPETRA
220#ifdef HAVE_MUELU_INST_DOUBLE_INT_INT
221 RCP< Xpetra::TpetraCrsGraph<LO, GO, Node> > tpCrsGraph = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsGraph<LO, GO, Node> >(crsGraph);
222 TEUCHOS_TEST_FOR_EXCEPTION(tpCrsGraph != Teuchos::null, Exceptions::RuntimeError, "Tpetra is not supported with Isorropia.");
223#else
224 TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "Isorropia is an interface to Zoltan which only has support for LO=GO=int and SC=double.");
225#endif // ENDIF HAVE_MUELU_INST_DOUBLE_INT_INT
226#endif // ENDIF HAVE_MUELU_TPETRA
227#endif // HAVE_MUELU_ISORROPIA
228#else // if we don't have MPI
229
230
231 // Running on one processor, so decomposition is the trivial one, all zeros.
232 RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
233 Set(level, "AmalgamatedPartition", decomposition);
234
235#endif // HAVE_MPI
236 // throw a more helpful error message if something failed
237 //TEUCHOS_TEST_FOR_EXCEPTION(!level.IsAvailable("AmalgamatedPartition"), Exceptions::RuntimeError, "IsorropiaInterface::Build : no \'Partition\' vector available on level. Isorropia failed to build a partition of the non-repartitioned graph of A. Please make sure, that Isorropia is correctly compiled (Epetra/Tpetra).");
238
239 } //Build()
240
241
242
243} //namespace MueLu
244
245//#endif //if defined(HAVE_MUELU_ISORROPIA) && defined(HAVE_MPI)
246
247
248#endif /* MUELU_ISORROPIAINTERFACE_DEF_HPP_ */
Exception indicating invalid cast attempted.
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.
void DeclareInput(Level &level) const
Specifies the data that this class needs, and the factories that generate that data.
void Build(Level &level) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Namespace for MueLu classes and methods.
@ Statistics0
Print statistics that do not involve significant additional computation.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.