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  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
101  RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, false);
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.
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.