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> >("AmalgamatedPartition", Teuchos::null, "(advanced) Factory generating the AmalgamatedPartition (e.g. an IsorropiaInterface)");
29 
30  return validParamList;
31  }
32 
33 
34  template <class LocalOrdinal, class GlobalOrdinal, class Node>
36  Input(currentLevel, "A");
37  Input(currentLevel, "AmalgamatedPartition");
38  } //DeclareInput()
39 
40  template <class LocalOrdinal, class GlobalOrdinal, class Node>
42  FactoryMonitor m(*this, "Build", level);
43 
44  RCP<Matrix> A = Get< RCP<Matrix> > (level, "A");
45  RCP<Xpetra::Vector<GO, LO, GO, NO> > amalgPartition = Get< RCP<Xpetra::Vector<GO, LO, GO, NO> > >(level, "AmalgamatedPartition");
46  //RCP<AmalgamationInfo> amalgInfo = Get< RCP<AmalgamationInfo> >(level, "UnAmalgamationInfo");
47 
48  RCP<const Teuchos::Comm< int > > comm = A->getRowMap()->getComm();
49 
50  ArrayRCP<GO> amalgPartitionData = amalgPartition->getDataNonConst(0);
51 
52  RCP<const Map> rowMap = A->getRowMap();
53  RCP<const Map> nodeMap = amalgPartition->getMap();
54 
55  // extract amalgamation information from matrix A
56  LO blockdim = 1; // block dim for fixed size blocks
57  GO offset = 0;
58  LO blockid = -1; // block id in strided map
59  LO nStridedOffset = 0; // DOF offset for strided block id "blockid" (default = 0)
60  LO stridedblocksize = blockdim; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
61 
62  // 1) check for blocking/striding information
63  // fill above variables
64  if(A->IsView("stridedMaps") &&
65  Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
66  Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
67  RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
68  TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,Exceptions::BadCast,"MueLu::RepartitionInterface::Build: cast to strided row map failed.");
69  blockdim = strMap->getFixedBlockSize();
70  offset = strMap->getOffset();
71  blockid = strMap->getStridedBlockId();
72  if (blockid > -1) {
73  std::vector<size_t> stridingInfo = strMap->getStridingData();
74  for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
75  nStridedOffset += stridingInfo[j];
76  stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
77 
78  } else {
79  stridedblocksize = blockdim;
80  }
81  oldView = A->SwitchToView(oldView);
82  GetOStream(Statistics0, -1) << "RepartitionInterface::Build():" << " found blockdim=" << blockdim << " from strided maps (blockid=" << blockid << ", strided block size=" << stridedblocksize << "). offset=" << offset << std::endl;
83  } else GetOStream(Statistics0, -1) << "RepartitionInterface::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
84 
85  // vector which stores final (unamalgamated) repartitioning
86  RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, false);
87  ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
88 
89  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!");
90 
91  //RCP<std::map<GO,std::vector<GO> > > nodegid2dofgids = amalgInfo->GetGlobalAmalgamationParams();
92 
93  // fill vector with information about partitioning
94  // TODO: we assume simple block maps here
95  // TODO: adapt this to usage of nodegid2dofgids
96  for(size_t i = 0; i < nodeMap->getNodeNumElements(); i++) {
97  // not fully sure about this. We're filling local ids in the decomposition vector with
98  // the results stored in array. The decomposition vector is created using the rowMap of A
99 
100  // transform local node id to global node id.
101  //GO gNodeId = nodeMap->getGlobalElement(i);
102 
103  // extract global DOF ids that belong to gNodeId
104  /*std::vector<GlobalOrdinal> DOFs = (*nodegid2dofgids)[gNodeId];
105  for(size_t j=0; j<stridedblocksize; j++) {
106  decompEntries[i*stridedblocksize + j] = myRank;
107  }*/
108  for (LO j = 0; j < stridedblocksize/*DOFs.size()*/; j++) {
109  // transform global DOF ids to local DOF ids using rowMap
110  // note: The vector decomposition is based on rowMap
111  //LO lDofId = rowMap->getLocalElement(DOFs[j]); // -> i doubt that we need this!
112 
113  // put the same domain id to all DOFs of the same node
114  decompEntries[i*stridedblocksize + j] = amalgPartitionData[i];
115  //decompEntries[lDofId] = amalgPartitionData[i];
116  }
117 
118  }
119 
120  Set(level, "Partition", decomposition);
121 
122  } //Build()
123 
124 
125 
126 } //namespace MueLu
127 
128 
129 #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.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Build(Level &level) const
Build an object with this factory.
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.