MueLu  Version of the Day
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
58  RCP<Xpetra::Vector<GO,LO,GO,NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
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  GO offset = 0;
75  LO blockid = -1; // block id in strided map
76  LO nStridedOffset = 0; // DOF offset for strided block id "blockid" (default = 0)
77  LO stridedblocksize = blockdim; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
78 
79  // 1) check for blocking/striding information
80  // fill above variables
81  if(A->IsView("stridedMaps") &&
82  Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
83  Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
84  RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
85  TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,Exceptions::BadCast,"MueLu::RepartitionInterface::Build: cast to strided row map failed.");
86  blockdim = strMap->getFixedBlockSize();
87  offset = strMap->getOffset();
88  blockid = strMap->getStridedBlockId();
89  if (blockid > -1) {
90  std::vector<size_t> stridingInfo = strMap->getStridingData();
91  for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
92  nStridedOffset += stridingInfo[j];
93  stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
94 
95  } else {
96  stridedblocksize = blockdim;
97  }
98  oldView = A->SwitchToView(oldView);
99  //GetOStream(Statistics0, -1) << "RepartitionInterface::Build():" << " found blockdim=" << blockdim << " from strided maps (blockid=" << blockid << ", strided block size=" << stridedblocksize << "). offset=" << offset << std::endl;
100  } else GetOStream(Statistics0, -1) << "RepartitionInterface::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
101 
102  // vector which stores final (unamalgamated) repartitioning
103  RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, false);
104  ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
105 
106  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!");
107 
108  //RCP<std::map<GO,std::vector<GO> > > nodegid2dofgids = amalgInfo->GetGlobalAmalgamationParams();
109 
110  // fill vector with information about partitioning
111  // TODO: we assume simple block maps here
112  // TODO: adapt this to usage of nodegid2dofgids
113  for(size_t i = 0; i < nodeMap->getNodeNumElements(); i++) {
114  // not fully sure about this. We're filling local ids in the decomposition vector with
115  // the results stored in array. The decomposition vector is created using the rowMap of A
116 
117  // transform local node id to global node id.
118  //GO gNodeId = nodeMap->getGlobalElement(i);
119 
120  // extract global DOF ids that belong to gNodeId
121  /*std::vector<GlobalOrdinal> DOFs = (*nodegid2dofgids)[gNodeId];
122  for(size_t j=0; j<stridedblocksize; j++) {
123  decompEntries[i*stridedblocksize + j] = myRank;
124  }*/
125  for (LO j = 0; j < stridedblocksize/*DOFs.size()*/; j++) {
126  // transform global DOF ids to local DOF ids using rowMap
127  // note: The vector decomposition is based on rowMap
128  //LO lDofId = rowMap->getLocalElement(DOFs[j]); // -> i doubt that we need this!
129 
130  // put the same domain id to all DOFs of the same node
131  decompEntries[i*stridedblocksize + j] = amalgPartitionData[i];
132  //decompEntries[lDofId] = amalgPartitionData[i];
133  }
134 
135  }
136 
137  Set(level, "Partition", decomposition);
138 
139  } //Build()
140 
141 
142 
143 } //namespace MueLu
144 
145 
146 #endif /* MUELU_REPARTITIONINTERFACE_DEF_HPP_ */
Exception indicating invalid cast attempted.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
Print statistics that do not involve significant additional computation.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Build(Level &level) const
Build an object with this factory.
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.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.