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