MueLu  Version of the Day
MueLu_UncoupledAggregationFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_UNCOUPLEDAGGREGATIONFACTORY_DEF_HPP_
47 #define MUELU_UNCOUPLEDAGGREGATIONFACTORY_DEF_HPP_
48 
49 #include <climits>
50 
51 #include <Xpetra_Map.hpp>
52 #include <Xpetra_Vector.hpp>
53 #include <Xpetra_MultiVectorFactory.hpp>
54 #include <Xpetra_VectorFactory.hpp>
55 
56 #include "MueLu_UncoupledAggregationFactory_decl.hpp"
57 
58 #include "MueLu_InterfaceAggregationAlgorithm.hpp"
59 #include "MueLu_OnePtAggregationAlgorithm.hpp"
60 #include "MueLu_PreserveDirichletAggregationAlgorithm.hpp"
61 #include "MueLu_IsolatedNodeAggregationAlgorithm.hpp"
62 
63 #include "MueLu_AggregationPhase1Algorithm.hpp"
64 #include "MueLu_AggregationPhase2aAlgorithm.hpp"
65 #include "MueLu_AggregationPhase2bAlgorithm.hpp"
66 #include "MueLu_AggregationPhase3Algorithm.hpp"
67 
68 #include "MueLu_Level.hpp"
69 #include "MueLu_GraphBase.hpp"
70 #include "MueLu_Aggregates.hpp"
71 #include "MueLu_MasterList.hpp"
72 #include "MueLu_Monitor.hpp"
73 #include "MueLu_AmalgamationInfo.hpp"
74 #include "MueLu_Utilities.hpp"
75 
76 namespace MueLu {
77 
78  template <class LocalOrdinal, class GlobalOrdinal, class Node>
79  UncoupledAggregationFactory<LocalOrdinal, GlobalOrdinal, Node>::UncoupledAggregationFactory()
80  : bDefinitionPhase_(true)
81  { }
82 
83  template <class LocalOrdinal, class GlobalOrdinal, class Node>
84  RCP<const ParameterList> UncoupledAggregationFactory<LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
85  RCP<ParameterList> validParamList = rcp(new ParameterList());
86 
87  // Aggregation parameters (used in aggregation algorithms)
88  // TODO introduce local member function for each aggregation algorithm such that each aggregation algorithm can define its own parameters
89 
90  typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
91 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
92  SET_VALID_ENTRY("aggregation: max agg size");
93  SET_VALID_ENTRY("aggregation: min agg size");
94  SET_VALID_ENTRY("aggregation: max selected neighbors");
95  SET_VALID_ENTRY("aggregation: ordering");
96  validParamList->getEntry("aggregation: ordering").setValidator(
97  rcp(new validatorType(Teuchos::tuple<std::string>("natural", "graph", "random"), "aggregation: ordering")));
98  SET_VALID_ENTRY("aggregation: enable phase 1");
99  SET_VALID_ENTRY("aggregation: enable phase 2a");
100  SET_VALID_ENTRY("aggregation: enable phase 2b");
101  SET_VALID_ENTRY("aggregation: enable phase 3");
102  SET_VALID_ENTRY("aggregation: preserve Dirichlet points");
103  SET_VALID_ENTRY("aggregation: allow user-specified singletons");
104  SET_VALID_ENTRY("aggregation: use interface aggregation");
105  SET_VALID_ENTRY("aggregation: error on nodes with no on-rank neighbors");
106  SET_VALID_ENTRY("aggregation: phase3 avoid singletons");
107 #undef SET_VALID_ENTRY
108 
109  // general variables needed in AggregationFactory
110  validParamList->set< RCP<const FactoryBase> >("Graph", null, "Generating factory of the graph");
111  validParamList->set< RCP<const FactoryBase> >("DofsPerNode", null, "Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
112 
113  // special variables necessary for OnePtAggregationAlgorithm
114  validParamList->set< std::string > ("OnePt aggregate map name", "", "Name of input map for single node aggregates. (default='')");
115  validParamList->set< std::string > ("OnePt aggregate map factory", "", "Generating factory of (DOF) map for single node aggregates.");
116  //validParamList->set< RCP<const FactoryBase> >("OnePt aggregate map factory", NoFactory::getRCP(), "Generating factory of (DOF) map for single node aggregates.");
117 
118  // InterfaceAggregation parameters
119  //validParamList->set< bool > ("aggregation: use interface aggregation", "false", "Flag to trigger aggregation along an interface using specified aggregate seeds.");
120  validParamList->set< std::string > ("Interface aggregate map name", "", "Name of input map for interface aggregates. (default='')");
121  validParamList->set< std::string > ("Interface aggregate map factory", "", "Generating factory of (DOF) map for interface aggregates.");
122  validParamList->set<RCP<const FactoryBase> > ("nodeOnInterface", Teuchos::null, "Array specifying whether or not a node is on the interface (1 or 0).");
123 
124  return validParamList;
125  }
126 
127  template <class LocalOrdinal, class GlobalOrdinal, class Node>
128  void UncoupledAggregationFactory<LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
129  Input(currentLevel, "Graph");
130  Input(currentLevel, "DofsPerNode");
131 
132  const ParameterList& pL = GetParameterList();
133 
134  // request special data necessary for OnePtAggregationAlgorithm
135  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
136  if (mapOnePtName.length() > 0) {
137  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
138  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
139  currentLevel.DeclareInput(mapOnePtName, NoFactory::get());
140  } else {
141  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
142  currentLevel.DeclareInput(mapOnePtName, mapOnePtFact.get());
143  }
144  }
145 
146  // request special data necessary for InterfaceAggregation
147  if (pL.get<bool>("aggregation: use interface aggregation") == true){
148  if(currentLevel.GetLevelID() == 0) {
149  if(currentLevel.IsAvailable("nodeOnInterface", NoFactory::get())) {
150  currentLevel.DeclareInput("nodeOnInterface", NoFactory::get(), this);
151  } else {
152  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("nodeOnInterface", NoFactory::get()),
153  Exceptions::RuntimeError,
154  "nodeOnInterface was not provided by the user on level0!");
155  }
156  } else {
157  Input(currentLevel, "nodeOnInterface");
158  }
159  }
160  }
161 
162  template <class LocalOrdinal, class GlobalOrdinal, class Node>
163  void UncoupledAggregationFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(Level &currentLevel) const {
164  FactoryMonitor m(*this, "Build", currentLevel);
165 
166  ParameterList pL = GetParameterList();
167  bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
168 
169  if (pL.get<int>("aggregation: max agg size") == -1)
170  pL.set("aggregation: max agg size", INT_MAX);
171 
172  // define aggregation algorithms
173  RCP<const FactoryBase> graphFact = GetFactory("Graph");
174 
175  // TODO Can we keep different aggregation algorithms over more Build calls?
176  algos_.clear();
177  algos_.push_back(rcp(new PreserveDirichletAggregationAlgorithm(graphFact)));
178  if (pL.get<bool>("aggregation: use interface aggregation") == true) algos_.push_back(rcp(new InterfaceAggregationAlgorithm (graphFact)));
179  if (pL.get<bool>("aggregation: allow user-specified singletons") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm (graphFact)));
180  if (pL.get<bool>("aggregation: enable phase 1" ) == true) algos_.push_back(rcp(new AggregationPhase1Algorithm (graphFact)));
181  if (pL.get<bool>("aggregation: enable phase 2a") == true) algos_.push_back(rcp(new AggregationPhase2aAlgorithm (graphFact)));
182  if (pL.get<bool>("aggregation: enable phase 2b") == true) algos_.push_back(rcp(new AggregationPhase2bAlgorithm (graphFact)));
183  if (pL.get<bool>("aggregation: enable phase 3" ) == true) algos_.push_back(rcp(new AggregationPhase3Algorithm (graphFact)));
184 
185  // TODO: remove old aggregation mode
186  //if (pL.get<bool>("UseOnePtAggregationAlgorithm") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm (graphFact)));
187  //if (pL.get<bool>("UseUncoupledAggregationAlgorithm") == true) algos_.push_back(rcp(new AggregationPhase1Algorithm (graphFact)));
188  //if (pL.get<bool>("UseMaxLinkAggregationAlgorithm") == true) algos_.push_back(rcp(new MaxLinkAggregationAlgorithm (graphFact)));
189  //if (pL.get<bool>("UseEmergencyAggregationAlgorithm") == true) algos_.push_back(rcp(new EmergencyAggregationAlgorithm (graphFact)));
190 
191  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
192  RCP<Map> OnePtMap = Teuchos::null;
193  if (mapOnePtName.length()) {
194  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
195  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
196  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, NoFactory::get());
197  } else {
198  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
199  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
200  }
201  }
202 
203  // Set map for interface aggregates
204  std::string mapInterfaceName = pL.get<std::string>("Interface aggregate map name");
205  RCP<Map> InterfaceMap = Teuchos::null;
206 
207  RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel, "Graph");
208 
209  // Build
210  RCP<Aggregates> aggregates = rcp(new Aggregates(*graph));
211  aggregates->setObjectLabel("UC");
212 
213  const LO numRows = graph->GetNodeNumVertices();
214 
215  // construct aggStat information
216  std::vector<unsigned> aggStat(numRows, READY);
217 
218  // interface
219  if (pL.get<bool>("aggregation: use interface aggregation") == true){
220  Teuchos::Array<LO> nodeOnInterface = Get<Array<LO>>(currentLevel,"nodeOnInterface");
221  for (LO i = 0; i < numRows; i++) {
222  if (nodeOnInterface[i])
223  aggStat[i] = INTERFACE;
224  }
225  }
226 
227  ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
228  if (dirichletBoundaryMap != Teuchos::null)
229  for (LO i = 0; i < numRows; i++)
230  if (dirichletBoundaryMap[i] == true)
231  aggStat[i] = BOUNDARY;
232 
233  LO nDofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
234  GO indexBase = graph->GetDomainMap()->getIndexBase();
235  if (OnePtMap != Teuchos::null) {
236  for (LO i = 0; i < numRows; i++) {
237  // reconstruct global row id (FIXME only works for contiguous maps)
238  GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
239 
240  for (LO kr = 0; kr < nDofsPerNode; kr++)
241  if (OnePtMap->isNodeGlobalElement(grid + kr))
242  aggStat[i] = ONEPT;
243  }
244  }
245 
246 
247 
248  const RCP<const Teuchos::Comm<int> > comm = graph->GetComm();
249  GO numGlobalRows = 0;
250  if (IsPrint(Statistics1))
251  MueLu_sumAll(comm, as<GO>(numRows), numGlobalRows);
252 
253  LO numNonAggregatedNodes = numRows;
254  GO numGlobalAggregatedPrev = 0, numGlobalAggsPrev = 0;
255  for (size_t a = 0; a < algos_.size(); a++) {
256  std::string phase = algos_[a]->description();
257  SubFactoryMonitor sfm(*this, "Algo \"" + phase + "\"", currentLevel);
258 
259  int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
260  algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
261  algos_[a]->SetProcRankVerbose(oldRank);
262 
263  if (IsPrint(Statistics1)) {
264  GO numLocalAggregated = numRows - numNonAggregatedNodes, numGlobalAggregated = 0;
265  GO numLocalAggs = aggregates->GetNumAggregates(), numGlobalAggs = 0;
266  MueLu_sumAll(comm, numLocalAggregated, numGlobalAggregated);
267  MueLu_sumAll(comm, numLocalAggs, numGlobalAggs);
268 
269  double aggPercent = 100*as<double>(numGlobalAggregated)/as<double>(numGlobalRows);
270  if (aggPercent > 99.99 && aggPercent < 100.00) {
271  // Due to round off (for instance, for 140465733/140466897), we could
272  // get 100.00% display even if there are some remaining nodes. This
273  // is bad from the users point of view. It is much better to change
274  // it to display 99.99%.
275  aggPercent = 99.99;
276  }
277  GetOStream(Statistics1) << " aggregated : " << (numGlobalAggregated - numGlobalAggregatedPrev) << " (phase), " << std::fixed
278  << std::setprecision(2) << numGlobalAggregated << "/" << numGlobalRows << " [" << aggPercent << "%] (total)\n"
279  << " remaining : " << numGlobalRows - numGlobalAggregated << "\n"
280  << " aggregates : " << numGlobalAggs-numGlobalAggsPrev << " (phase), " << numGlobalAggs << " (total)" << std::endl;
281  numGlobalAggregatedPrev = numGlobalAggregated;
282  numGlobalAggsPrev = numGlobalAggs;
283  }
284  }
285 
286  TEUCHOS_TEST_FOR_EXCEPTION(numNonAggregatedNodes, Exceptions::RuntimeError, "MueLu::UncoupledAggregationFactory::Build: Leftover nodes found! Error!");
287 
288  aggregates->AggregatesCrossProcessors(false);
289  aggregates->ComputeAggregateSizes(true/*forceRecompute*/);
290 
291  Set(currentLevel, "Aggregates", aggregates);
292 
293  GetOStream(Statistics1) << aggregates->description() << std::endl;
294  }
295 
296 } //namespace MueLu
297 
298 
299 #endif /* MUELU_UNCOUPLEDAGGREGATIONFACTORY_DEF_HPP_ */
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
MueLu_sumAll
#define MueLu_sumAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:61
MueLu::UncoupledAggregationFactory::DeclareInput
void DeclareInput(Level &currentLevel) const
Input.
Definition: MueLu_UncoupledAggregationFactory_def.hpp:128
string
DASHPOS string(SUBSTRING ${CLASS} 0 ${DASHPOS} CLASS_NAME) set(CONDITION_NAME "") set(CONDITION_NAME_END "") string(SUBSTRING $
Definition: src/CMakeLists.txt:134
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:209
MueLu::UncoupledAggregationFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_UncoupledAggregationFactory_def.hpp:84
MueLu_Monitor.hpp
MueLu::Level::DeclareInput
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
Definition: MueLu_Level.cpp:150
MueLu::INTERFACE
@ INTERFACE
Definition: MueLu_Types.hpp:85
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
MueLu::UncoupledAggregationFactory::UncoupledAggregationFactory
UncoupledAggregationFactory()
Constructor.
Definition: MueLu_UncoupledAggregationFactory_def.hpp:79
MueLu_UncoupledAggregationFactory_decl.hpp
MueLu::PreserveDirichletAggregationAlgorithm
Builds one-to-one aggregates for all Dirichlet boundary nodes. For some applications this might be ne...
Definition: MueLu_PreserveDirichletAggregationAlgorithm_decl.hpp:86
MueLu_Level.hpp
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition: MueLu_Monitor.hpp:261
MueLu::OnePtAggregationAlgorithm
Algorithm for coarsening a graph with uncoupled aggregation. keep special marked nodes as singleton n...
Definition: MueLu_OnePtAggregationAlgorithm_decl.hpp:88
MueLu::Level::IsAvailable
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
Definition: MueLu_Level.hpp:333
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::Level::Get
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
Definition: MueLu_Level.hpp:192
MueLu::Statistics1
@ Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:71
MueLu::Aggregates
Container class for aggregation information.
Definition: MueLu_Aggregates_decl.hpp:102
MueLu::ONEPT
@ ONEPT
Definition: MueLu_Types.hpp:76
MueLu::AggregationPhase1Algorithm
Algorithm for coarsening a graph with uncoupled aggregation.
Definition: MueLu_AggregationPhase1Algorithm_decl.hpp:86
MueLu::UncoupledAggregationFactory::Build
void Build(Level &currentLevel) const
Build aggregates.
Definition: MueLu_UncoupledAggregationFactory_def.hpp:163
MueLu::AggregationPhase3Algorithm
Handle leftover nodes. Try to avoid singleton nodes.
Definition: MueLu_AggregationPhase3Algorithm_decl.hpp:79
MueLu::AggregationPhase2aAlgorithm
Among unaggregated points, see if we can make a reasonable size aggregate out of it.
Definition: MueLu_AggregationPhase2aAlgorithm_decl.hpp:85
MueLu::AggregationPhase2bAlgorithm
Add leftovers to existing aggregates.
Definition: MueLu_AggregationPhase2bAlgorithm_decl.hpp:84
SET_VALID_ENTRY
#define SET_VALID_ENTRY(name)
MueLu::InterfaceAggregationAlgorithm
Algorithm for coarsening a graph with uncoupled aggregation. creates aggregates along an interface us...
Definition: MueLu_InterfaceAggregationAlgorithm_decl.hpp:85
MueLu_GraphBase.hpp
MueLu_MasterList.hpp
MueLu::READY
@ READY
Definition: MueLu_Types.hpp:65
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:59
MueLu::BOUNDARY
@ BOUNDARY
Definition: MueLu_Types.hpp:82
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99