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