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 #include <Xpetra_BlockedMap.hpp>
21 
22 #ifdef HAVE_MUELU_ISORROPIA
23 #include <Isorropia_Exception.hpp>
24 
25 
26 #ifdef HAVE_MUELU_EPETRA
28 #include <Epetra_CrsGraph.h>
29 #include <Isorropia_EpetraPartitioner.hpp>
30 #endif
31 
32 #ifdef HAVE_MUELU_TPETRA
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_AmalgamationFactory.hpp"
42 #include "MueLu_AmalgamationInfo.hpp"
43 #include "MueLu_Utilities.hpp"
44 
45 namespace MueLu {
46 
47  template <class LocalOrdinal, class GlobalOrdinal, class Node>
49  RCP<ParameterList> validParamList = rcp(new ParameterList());
50 
51  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Factory of the matrix A");
52  validParamList->set< RCP<const FactoryBase> >("number of partitions", Teuchos::null, "Instance of RepartitionHeuristicFactory.");
53  validParamList->set< RCP<const FactoryBase> >("UnAmalgamationInfo", Teuchos::null, "Generating factory of UnAmalgamationInfo");
54 
55  return validParamList;
56  }
57 
58 
59  template <class LocalOrdinal, class GlobalOrdinal, class Node>
61  Input(currentLevel, "A");
62  Input(currentLevel, "number of partitions");
63  Input(currentLevel, "UnAmalgamationInfo");
64  }
65 
66  template <class LocalOrdinal, class GlobalOrdinal, class Node>
68  FactoryMonitor m(*this, "Build", level);
70 
71  RCP<Matrix> A = Get< RCP<Matrix> >(level, "A");
72  RCP<AmalgamationInfo> amalInfo = Get< RCP<AmalgamationInfo> >(level, "UnAmalgamationInfo");
73  GO numParts = Get< int >(level, "number of partitions");
74 
75  RCP<const Map> rowMap = A->getRowMap();
76  RCP<const Map> colMap = A->getColMap();
77 
78  if (numParts == 1 || numParts == -1) {
79  // Running on one processor, so decomposition is the trivial one, all zeros.
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) << "IsorropiaInterface::Build():" << " found blockdim=" << blockdim << " from strided maps (blockid=" << blockid << ", strided block size=" << stridedblocksize << "). offset=" << offset << std::endl;
125  } else GetOStream(Statistics0) << "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  RCP<const BlockedMap> bnodeMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(nodeMap);
131  if(!bnodeMap.is_null()) nodeMap=bnodeMap->getMap();
132 
133  GetOStream(Statistics0) << "IsorropiaInterface:Build(): nodeMap " << nodeMap->getNodeNumElements() << "/" << nodeMap->getGlobalNumElements() << " elements" << std::endl;
134 
135 
136  // 3) create graph of amalgamated matrix
137  RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, 10, Xpetra::DynamicProfile);
138 
139  // 4) do amalgamation. generate graph of amalgamated matrix
140  for(LO row=0; row<Teuchos::as<LO>(A->getRowMap()->getNodeNumElements()); row++) {
141  // get global DOF id
142  GO grid = rowMap->getGlobalElement(row);
143 
144  // translate grid to nodeid
145  GO nodeId = AmalgamationFactory::DOFGid2NodeId(grid, blockdim, offset, indexBase);
146 
147  size_t nnz = A->getNumEntriesInLocalRow(row);
150  A->getLocalRowView(row, indices, vals);
151 
152  RCP<std::vector<GO> > cnodeIds = Teuchos::rcp(new std::vector<GO>); // global column block ids
153  LO realnnz = 0;
154  for(LO col=0; col<Teuchos::as<LO>(nnz); col++) {
155  GO gcid = colMap->getGlobalElement(indices[col]); // global column id
156 
157  if(vals[col]!=0.0) {
158  GO cnodeId = AmalgamationFactory::DOFGid2NodeId(gcid, blockdim, offset, indexBase);
159  cnodeIds->push_back(cnodeId);
160  realnnz++; // increment number of nnz in matrix row
161  }
162  }
163 
164  Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp( cnodeIds );
165 
166  if(arr_cnodeIds.size() > 0 )
167  crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
168  }
169  // fill matrix graph
170  crsGraph->fillComplete(nodeMap,nodeMap);
171 
172 #ifdef HAVE_MPI
173 #ifdef HAVE_MUELU_ISORROPIA
174 
175  // prepare parameter list for Isorropia
176  Teuchos::ParameterList paramlist;
177  paramlist.set("NUM PARTS", toString(numParts));
178 
179  /*STRUCTURALLY SYMMETRIC [NO/yes] (is symmetrization required?)
180  PARTITIONING METHOD [block/cyclic/random/rcb/rib/hsfc/graph/HYPERGRAPH]
181  NUM PARTS [int k] (global number of parts)
182  IMBALANCE TOL [float tol] (1.0 is perfect balance)
183  BALANCE OBJECTIVE [ROWS/nonzeros]
184  */
185  Teuchos::ParameterList& sublist = paramlist.sublist("Zoltan");
186  sublist.set("LB_APPROACH", "PARTITION");
187 
188 
189 
190 #ifdef HAVE_MUELU_EPETRA
191  RCP< Xpetra::EpetraCrsGraphT<GO, Node> > epCrsGraph = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsGraphT<GO, Node> >(crsGraph);
192  if(epCrsGraph != Teuchos::null) {
193  RCP< const Epetra_CrsGraph> epetraCrsGraph = epCrsGraph->getEpetra_CrsGraph();
194 
195  RCP<Isorropia::Epetra::Partitioner> isoPart = Teuchos::rcp(new Isorropia::Epetra::Partitioner(epetraCrsGraph,paramlist));
196 
197  int size = 0;
198  const int* array = NULL;
199  isoPart->extractPartsView(size,array);
200 
201  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");
202 
204  ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
205 
206  // fill vector with amalgamated information about partitioning
207  for(int i = 0; i<size; i++) {
208  decompEntries[i] = Teuchos::as<GO>(array[i]);
209  }
210 
211  Set(level, "AmalgamatedPartition", decomposition);
212 
213  }
214 #endif // ENDIF HAVE_MUELU_EPETRA
215 
216 #ifdef HAVE_MUELU_TPETRA
217 #ifdef HAVE_MUELU_INST_DOUBLE_INT_INT
218  RCP< Xpetra::TpetraCrsGraph<LO, GO, Node> > tpCrsGraph = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsGraph<LO, GO, Node> >(crsGraph);
219  TEUCHOS_TEST_FOR_EXCEPTION(tpCrsGraph != Teuchos::null, Exceptions::RuntimeError, "Tpetra is not supported with Isorropia.");
220 #else
221  TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "Isorropia is an interface to Zoltan which only has support for LO=GO=int and SC=double.");
222 #endif // ENDIF HAVE_MUELU_INST_DOUBLE_INT_INT
223 #endif // ENDIF HAVE_MUELU_TPETRA
224 #endif // HAVE_MUELU_ISORROPIA
225 #else // if we don't have MPI
226 
227 
228  // Running on one processor, so decomposition is the trivial one, all zeros.
230  Set(level, "AmalgamatedPartition", decomposition);
231 
232 #endif // HAVE_MPI
233  // throw a more helpful error message if something failed
234  //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).");
235 
236  } //Build()
237 
238 
239 
240 } //namespace MueLu
241 
242 //#endif //if defined(HAVE_MUELU_ISORROPIA) && defined(HAVE_MPI)
243 
244 
245 #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
GlobalOrdinal GO
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.
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
bool is_null() const
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
std::string viewLabel_t
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
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.