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> >("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);
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.
size_type size() const
void DeclareInput(Level &level) const
Specifies the data that this class needs, and the factories that generate that data.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Build(Level &level) const
Build an object with this factory.
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Exception throws to report errors in the internal logical of the program.