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 
17 namespace MueLu {
18 
19 template <class LocalOrdinal, class GlobalOrdinal, class Node>
21  RCP<ParameterList> validParamList = rcp(new ParameterList());
22  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Factory of the matrix A");
23  validParamList->set<RCP<const FactoryBase> >("number of partitions", Teuchos::null, "Instance of RepartitionHeuristicFactory.");
24  validParamList->set<RCP<const FactoryBase> >("AmalgamatedPartition", Teuchos::null, "(advanced) Factory generating the AmalgamatedPartition (e.g. an IsorropiaInterface)");
25 
26  return validParamList;
27 }
28 
29 template <class LocalOrdinal, class GlobalOrdinal, class Node>
31  Input(currentLevel, "A");
32  Input(currentLevel, "number of partitions");
33  Input(currentLevel, "AmalgamatedPartition");
34 } // DeclareInput()
35 
36 template <class LocalOrdinal, class GlobalOrdinal, class Node>
38  FactoryMonitor m(*this, "Build", level);
39 
40  RCP<Matrix> A = Get<RCP<Matrix> >(level, "A");
41  RCP<Xpetra::Vector<GO, LO, GO, NO> > amalgPartition = Get<RCP<Xpetra::Vector<GO, LO, GO, NO> > >(level, "AmalgamatedPartition");
42  int numParts = Get<int>(level, "number of partitions");
43 
44  RCP<const Map> rowMap = A->getRowMap();
45 
46  // standard case: use matrix info and amalgamated rebalancing info to create "Partition" vector
47  RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
48 
49  // Short cut: if we only need one partition, then create a dummy partition vector
50  if (numParts == 1 || numParts == -1) {
51  // Single processor, decomposition is trivial: all zeros
53  Set(level, "Partition", decomposition);
54  return;
55  } /* else if (numParts == -1) {
56  // No repartitioning
57  RCP<Xpetra::Vector<GO,LO,GO,NO> > decomposition = Teuchos::null; //Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
58  //decomposition->putScalar(Teuchos::as<Scalar>(comm->getRank()));
59  Set(level, "Partition", decomposition);
60  return;
61  }*/
62 
63  ArrayRCP<GO> amalgPartitionData = amalgPartition->getDataNonConst(0);
64  RCP<const Map> nodeMap = amalgPartition->getMap();
65 
66  // extract amalgamation information from matrix A
67  LO blockdim = 1; // block dim for fixed size blocks
68  LO blockid = -1; // block id in strided map
69  LO nStridedOffset = 0; // DOF offset for strided block id "blockid" (default = 0)
70  LO stridedblocksize = blockdim; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
71 
72  // 1) check for blocking/striding information
73  // fill above variables
74  if (A->IsView("stridedMaps") &&
75  Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
76  Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
77  RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
78  TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null, Exceptions::BadCast, "MueLu::RepartitionInterface::Build: cast to strided row map failed.");
79  blockdim = strMap->getFixedBlockSize();
80  blockid = strMap->getStridedBlockId();
81  if (blockid > -1) {
82  std::vector<size_t> stridingInfo = strMap->getStridingData();
83  for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
84  nStridedOffset += stridingInfo[j];
85  stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
86 
87  } else {
88  stridedblocksize = blockdim;
89  }
90  oldView = A->SwitchToView(oldView);
91  // GetOStream(Statistics0, -1) << "RepartitionInterface::Build():" << " found blockdim=" << blockdim << " from strided maps (blockid=" << blockid << ", strided block size=" << stridedblocksize << "). offset=" << offset << std::endl;
92  } else
93  GetOStream(Statistics0, -1) << "RepartitionInterface::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
94 
95  // vector which stores final (unamalgamated) repartitioning
97  ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
98 
99  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<int>(nodeMap->getLocalNumElements()) * stridedblocksize != Teuchos::as<int>(rowMap->getLocalNumElements()), Exceptions::RuntimeError, "Inconsistency between nodeMap and dofMap: we are supporting block maps only. No support for general strided maps, yet!");
100 
101  // RCP<std::map<GO,std::vector<GO> > > nodegid2dofgids = amalgInfo->GetGlobalAmalgamationParams();
102 
103  // fill vector with information about partitioning
104  // TODO: we assume simple block maps here
105  // TODO: adapt this to usage of nodegid2dofgids
106  for (size_t i = 0; i < nodeMap->getLocalNumElements(); i++) {
107  // not fully sure about this. We're filling local ids in the decomposition vector with
108  // the results stored in array. The decomposition vector is created using the rowMap of A
109 
110  // transform local node id to global node id.
111  // GO gNodeId = nodeMap->getGlobalElement(i);
112 
113  // extract global DOF ids that belong to gNodeId
114  /*std::vector<GlobalOrdinal> DOFs = (*nodegid2dofgids)[gNodeId];
115  for(size_t j=0; j<stridedblocksize; j++) {
116  decompEntries[i*stridedblocksize + j] = myRank;
117  }*/
118  for (LO j = 0; j < stridedblocksize /*DOFs.size()*/; j++) {
119  // transform global DOF ids to local DOF ids using rowMap
120  // note: The vector decomposition is based on rowMap
121  // LO lDofId = rowMap->getLocalElement(DOFs[j]); // -> i doubt that we need this!
122 
123  // put the same domain id to all DOFs of the same node
124  decompEntries[i * stridedblocksize + j] = amalgPartitionData[i];
125  // decompEntries[lDofId] = amalgPartitionData[i];
126  }
127  }
128 
129  Set(level, "Partition", decomposition);
130 
131 } // Build()
132 
133 } // namespace MueLu
134 
135 #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.