MueLu Version of the Day
Loading...
Searching...
No Matches
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
22namespace 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->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!");
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->getLocalNumElements(); 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.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
void DeclareInput(Level &level) const
Specifies the data that this class needs, and the factories that generate that data.
void Build(Level &level) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Namespace for MueLu classes and methods.
@ Statistics0
Print statistics that do not involve significant additional computation.