MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_RepartitionInterface_def.hpp
Go to the documentation of this file.
1 /*
2  * MueLu_RepartitionInterface_def.hpp
3  *
4  * Created on: 5 Sep 2013
5  * Author: wiesner
6  */
7 
8 #ifndef MUELU_REPARTITIONINTERFACE_DEF_HPP_
9 #define MUELU_REPARTITIONINTERFACE_DEF_HPP_
10 
12 
13 #include "MueLu_Level.hpp"
14 #include "MueLu_Exceptions.hpp"
15 #include "MueLu_Monitor.hpp"
16 #include "MueLu_Graph.hpp"
17 #include "MueLu_AmalgamationFactory.hpp"
18 #include "MueLu_AmalgamationInfo.hpp"
19 #include "MueLu_Utilities.hpp"
20 
21 
22 namespace MueLu {
23 
24  template <class LocalOrdinal, class GlobalOrdinal, class Node>
26  RCP<ParameterList> validParamList = rcp(new ParameterList());
27  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Factory of the matrix A");
28  validParamList->set< RCP<const FactoryBase> >("number of partitions", Teuchos::null, "Instance of RepartitionHeuristicFactory.");
29  validParamList->set< RCP<const FactoryBase> >("AmalgamatedPartition", Teuchos::null, "(advanced) Factory generating the AmalgamatedPartition (e.g. an IsorropiaInterface)");
30 
31  return validParamList;
32  }
33 
34 
35  template <class LocalOrdinal, class GlobalOrdinal, class Node>
37  Input(currentLevel, "A");
38  Input(currentLevel, "number of partitions");
39  Input(currentLevel, "AmalgamatedPartition");
40  } //DeclareInput()
41 
42  template <class LocalOrdinal, class GlobalOrdinal, class Node>
44  FactoryMonitor m(*this, "Build", level);
45 
46  RCP<Matrix> A = Get< RCP<Matrix> > (level, "A");
47  RCP<Xpetra::Vector<GO, LO, GO, NO> > amalgPartition = Get< RCP<Xpetra::Vector<GO, LO, GO, NO> > >(level, "AmalgamatedPartition");
48  int numParts = Get<int>(level, "number of partitions");
49 
50  RCP<const Map> rowMap = A->getRowMap();
51 
52  // standard case: use matrix info and amalgamated rebalancing info to create "Partition" vector
53  RCP<const Teuchos::Comm< int > > comm = A->getRowMap()->getComm();
54 
55  // Short cut: if we only need one partition, then create a dummy partition vector
56  if (numParts == 1 || numParts == -1) {
57  // Single processor, decomposition is trivial: all zeros
59  Set(level, "Partition", decomposition);
60  return;
61  }/* else if (numParts == -1) {
62  // No repartitioning
63  RCP<Xpetra::Vector<GO,LO,GO,NO> > decomposition = Teuchos::null; //Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
64  //decomposition->putScalar(Teuchos::as<Scalar>(comm->getRank()));
65  Set(level, "Partition", decomposition);
66  return;
67  }*/
68 
69  ArrayRCP<GO> amalgPartitionData = amalgPartition->getDataNonConst(0);
70  RCP<const Map> nodeMap = amalgPartition->getMap();
71 
72  // extract amalgamation information from matrix A
73  LO blockdim = 1; // block dim for fixed size blocks
74  LO blockid = -1; // block id in strided map
75  LO nStridedOffset = 0; // DOF offset for strided block id "blockid" (default = 0)
76  LO stridedblocksize = blockdim; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
77 
78  // 1) check for blocking/striding information
79  // fill above variables
80  if(A->IsView("stridedMaps") &&
81  Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
82  Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
83  RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
84  TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,Exceptions::BadCast,"MueLu::RepartitionInterface::Build: cast to strided row map failed.");
85  blockdim = strMap->getFixedBlockSize();
86  blockid = strMap->getStridedBlockId();
87  if (blockid > -1) {
88  std::vector<size_t> stridingInfo = strMap->getStridingData();
89  for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
90  nStridedOffset += stridingInfo[j];
91  stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
92 
93  } else {
94  stridedblocksize = blockdim;
95  }
96  oldView = A->SwitchToView(oldView);
97  //GetOStream(Statistics0, -1) << "RepartitionInterface::Build():" << " found blockdim=" << blockdim << " from strided maps (blockid=" << blockid << ", strided block size=" << stridedblocksize << "). offset=" << offset << std::endl;
98  } else GetOStream(Statistics0, -1) << "RepartitionInterface::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
99 
100  // vector which stores final (unamalgamated) repartitioning
102  ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
103 
104  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<int>(nodeMap->getNodeNumElements())*stridedblocksize != Teuchos::as<int>(rowMap->getNodeNumElements()), Exceptions::RuntimeError, "Inconsistency between nodeMap and dofMap: we are supporting block maps only. No support for general strided maps, yet!");
105 
106  //RCP<std::map<GO,std::vector<GO> > > nodegid2dofgids = amalgInfo->GetGlobalAmalgamationParams();
107 
108  // fill vector with information about partitioning
109  // TODO: we assume simple block maps here
110  // TODO: adapt this to usage of nodegid2dofgids
111  for(size_t i = 0; i < nodeMap->getNodeNumElements(); i++) {
112  // not fully sure about this. We're filling local ids in the decomposition vector with
113  // the results stored in array. The decomposition vector is created using the rowMap of A
114 
115  // transform local node id to global node id.
116  //GO gNodeId = nodeMap->getGlobalElement(i);
117 
118  // extract global DOF ids that belong to gNodeId
119  /*std::vector<GlobalOrdinal> DOFs = (*nodegid2dofgids)[gNodeId];
120  for(size_t j=0; j<stridedblocksize; j++) {
121  decompEntries[i*stridedblocksize + j] = myRank;
122  }*/
123  for (LO j = 0; j < stridedblocksize/*DOFs.size()*/; j++) {
124  // transform global DOF ids to local DOF ids using rowMap
125  // note: The vector decomposition is based on rowMap
126  //LO lDofId = rowMap->getLocalElement(DOFs[j]); // -> i doubt that we need this!
127 
128  // put the same domain id to all DOFs of the same node
129  decompEntries[i*stridedblocksize + j] = amalgPartitionData[i];
130  //decompEntries[lDofId] = amalgPartitionData[i];
131  }
132 
133  }
134 
135  Set(level, "Partition", decomposition);
136 
137  } //Build()
138 
139 
140 
141 } //namespace MueLu
142 
143 
144 #endif /* MUELU_REPARTITIONINTERFACE_DEF_HPP_ */
Exception indicating invalid cast attempted.
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)
LocalOrdinal LO
Print statistics that do not involve significant additional computation.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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
void Build(Level &level) const
Build an object with this factory.
std::string viewLabel_t
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
Exception throws to report errors in the internal logical of the program.
void DeclareInput(Level &level) const
Specifies the data that this class needs, and the factories that generate that data.