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 #undef SET_VALID_ENTRY
107 
108  // general variables needed in AggregationFactory
109  validParamList->set< RCP<const FactoryBase> >("Graph", null, "Generating factory of the graph");
110  validParamList->set< RCP<const FactoryBase> >("DofsPerNode", null, "Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
111 
112  // special variables necessary for OnePtAggregationAlgorithm
113  validParamList->set< std::string > ("OnePt aggregate map name", "", "Name of input map for single node aggregates. (default='')");
114  validParamList->set< std::string > ("OnePt aggregate map factory", "", "Generating factory of (DOF) map for single node aggregates.");
115  //validParamList->set< RCP<const FactoryBase> >("OnePt aggregate map factory", NoFactory::getRCP(), "Generating factory of (DOF) map for single node aggregates.");
116 
117  // InterfaceAggregation parameters
118  //validParamList->set< bool > ("aggregation: use interface aggregation", "false", "Flag to trigger aggregation along an interface using specified aggregate seeds.");
119  validParamList->set< std::string > ("Interface aggregate map name", "", "Name of input map for interface aggregates. (default='')");
120  validParamList->set< std::string > ("Interface aggregate map factory", "", "Generating factory of (DOF) map for interface aggregates.");
121  validParamList->set<RCP<const FactoryBase> > ("nodeOnInterface", Teuchos::null, "Array specifying whether or not a node is on the interface (1 or 0).");
122 
123  return validParamList;
124  }
125 
126  template <class LocalOrdinal, class GlobalOrdinal, class Node>
127  void UncoupledAggregationFactory<LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
128  Input(currentLevel, "Graph");
129  Input(currentLevel, "DofsPerNode");
130 
131  const ParameterList& pL = GetParameterList();
132 
133  // request special data necessary for OnePtAggregationAlgorithm
134  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
135  if (mapOnePtName.length() > 0) {
136  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
137  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
138  currentLevel.DeclareInput(mapOnePtName, NoFactory::get());
139  } else {
140  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
141  currentLevel.DeclareInput(mapOnePtName, mapOnePtFact.get());
142  }
143  }
144 
145  // request special data necessary for InterfaceAggregation
146  if (pL.get<bool>("aggregation: use interface aggregation") == true){
147  if(currentLevel.GetLevelID() == 0) {
148  if(currentLevel.IsAvailable("nodeOnInterface", NoFactory::get())) {
149  currentLevel.DeclareInput("nodeOnInterface", NoFactory::get(), this);
150  } else {
151  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("nodeOnInterface", NoFactory::get()),
152  Exceptions::RuntimeError,
153  "nodeOnInterface was not provided by the user on level0!");
154  }
155  } else {
156  Input(currentLevel, "nodeOnInterface");
157  }
158  }
159  }
160 
161  template <class LocalOrdinal, class GlobalOrdinal, class Node>
162  void UncoupledAggregationFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(Level &currentLevel) const {
163  FactoryMonitor m(*this, "Build", currentLevel);
164 
165  ParameterList pL = GetParameterList();
166  bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
167 
168  if (pL.get<int>("aggregation: max agg size") == -1)
169  pL.set("aggregation: max agg size", INT_MAX);
170 
171  // define aggregation algorithms
172  RCP<const FactoryBase> graphFact = GetFactory("Graph");
173 
174  // TODO Can we keep different aggregation algorithms over more Build calls?
175  algos_.clear();
176  algos_.push_back(rcp(new PreserveDirichletAggregationAlgorithm(graphFact)));
177  if (pL.get<bool>("aggregation: use interface aggregation") == true) algos_.push_back(rcp(new InterfaceAggregationAlgorithm (graphFact)));
178  if (pL.get<bool>("aggregation: allow user-specified singletons") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm (graphFact)));
179  if (pL.get<bool>("aggregation: enable phase 1" ) == true) algos_.push_back(rcp(new AggregationPhase1Algorithm (graphFact)));
180  if (pL.get<bool>("aggregation: enable phase 2a") == true) algos_.push_back(rcp(new AggregationPhase2aAlgorithm (graphFact)));
181  if (pL.get<bool>("aggregation: enable phase 2b") == true) algos_.push_back(rcp(new AggregationPhase2bAlgorithm (graphFact)));
182  if (pL.get<bool>("aggregation: enable phase 3" ) == true) algos_.push_back(rcp(new AggregationPhase3Algorithm (graphFact)));
183 
184  // TODO: remove old aggregation mode
185  //if (pL.get<bool>("UseOnePtAggregationAlgorithm") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm (graphFact)));
186  //if (pL.get<bool>("UseUncoupledAggregationAlgorithm") == true) algos_.push_back(rcp(new AggregationPhase1Algorithm (graphFact)));
187  //if (pL.get<bool>("UseMaxLinkAggregationAlgorithm") == true) algos_.push_back(rcp(new MaxLinkAggregationAlgorithm (graphFact)));
188  //if (pL.get<bool>("UseEmergencyAggregationAlgorithm") == true) algos_.push_back(rcp(new EmergencyAggregationAlgorithm (graphFact)));
189 
190  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
191  RCP<Map> OnePtMap = Teuchos::null;
192  if (mapOnePtName.length()) {
193  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
194  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
195  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, NoFactory::get());
196  } else {
197  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
198  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
199  }
200  }
201 
202  // Set map for interface aggregates
203  std::string mapInterfaceName = pL.get<std::string>("Interface aggregate map name");
204  RCP<Map> InterfaceMap = Teuchos::null;
205 
206  RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel, "Graph");
207 
208  // Build
209  RCP<Aggregates> aggregates = rcp(new Aggregates(*graph));
210  aggregates->setObjectLabel("UC");
211 
212  const LO numRows = graph->GetNodeNumVertices();
213 
214  // construct aggStat information
215  std::vector<unsigned> aggStat(numRows, READY);
216 
217  // interface
218  if (pL.get<bool>("aggregation: use interface aggregation") == true){
219  Teuchos::Array<LO> nodeOnInterface = Get<Array<LO>>(currentLevel,"nodeOnInterface");
220  for (LO i = 0; i < numRows; i++) {
221  if (nodeOnInterface[i])
222  aggStat[i] = INTERFACE;
223  }
224  }
225 
226  ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
227  if (dirichletBoundaryMap != Teuchos::null)
228  for (LO i = 0; i < numRows; i++)
229  if (dirichletBoundaryMap[i] == true)
230  aggStat[i] = BOUNDARY;
231 
232  LO nDofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
233  GO indexBase = graph->GetDomainMap()->getIndexBase();
234  if (OnePtMap != Teuchos::null) {
235  for (LO i = 0; i < numRows; i++) {
236  // reconstruct global row id (FIXME only works for contiguous maps)
237  GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
238 
239  for (LO kr = 0; kr < nDofsPerNode; kr++)
240  if (OnePtMap->isNodeGlobalElement(grid + kr))
241  aggStat[i] = ONEPT;
242  }
243  }
244 
245 
246 
247  const RCP<const Teuchos::Comm<int> > comm = graph->GetComm();
248  GO numGlobalRows = 0;
249  if (IsPrint(Statistics1))
250  MueLu_sumAll(comm, as<GO>(numRows), numGlobalRows);
251 
252  LO numNonAggregatedNodes = numRows;
253  GO numGlobalAggregatedPrev = 0, numGlobalAggsPrev = 0;
254  for (size_t a = 0; a < algos_.size(); a++) {
255  std::string phase = algos_[a]->description();
256  SubFactoryMonitor sfm(*this, "Algo \"" + phase + "\"", currentLevel);
257 
258  int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
259  algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
260  algos_[a]->SetProcRankVerbose(oldRank);
261 
262  if (IsPrint(Statistics1)) {
263  GO numLocalAggregated = numRows - numNonAggregatedNodes, numGlobalAggregated = 0;
264  GO numLocalAggs = aggregates->GetNumAggregates(), numGlobalAggs = 0;
265  MueLu_sumAll(comm, numLocalAggregated, numGlobalAggregated);
266  MueLu_sumAll(comm, numLocalAggs, numGlobalAggs);
267 
268  double aggPercent = 100*as<double>(numGlobalAggregated)/as<double>(numGlobalRows);
269  if (aggPercent > 99.99 && aggPercent < 100.00) {
270  // Due to round off (for instance, for 140465733/140466897), we could
271  // get 100.00% display even if there are some remaining nodes. This
272  // is bad from the users point of view. It is much better to change
273  // it to display 99.99%.
274  aggPercent = 99.99;
275  }
276  GetOStream(Statistics1) << " aggregated : " << (numGlobalAggregated - numGlobalAggregatedPrev) << " (phase), " << std::fixed
277  << std::setprecision(2) << numGlobalAggregated << "/" << numGlobalRows << " [" << aggPercent << "%] (total)\n"
278  << " remaining : " << numGlobalRows - numGlobalAggregated << "\n"
279  << " aggregates : " << numGlobalAggs-numGlobalAggsPrev << " (phase), " << numGlobalAggs << " (total)" << std::endl;
280  numGlobalAggregatedPrev = numGlobalAggregated;
281  numGlobalAggsPrev = numGlobalAggs;
282  }
283  }
284 
285  TEUCHOS_TEST_FOR_EXCEPTION(numNonAggregatedNodes, Exceptions::RuntimeError, "MueLu::UncoupledAggregationFactory::Build: Leftover nodes found! Error!");
286 
287  aggregates->AggregatesCrossProcessors(false);
288  aggregates->ComputeAggregateSizes(true/*forceRecompute*/);
289 
290  Set(currentLevel, "Aggregates", aggregates);
291 
292  GetOStream(Statistics1) << aggregates->description() << std::endl;
293  }
294 
295 } //namespace MueLu
296 
297 
298 #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:127
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:228
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
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:85
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:301
MueLu::OnePtAggregationAlgorithm
Algorithm for coarsening a graph with uncoupled aggregation. keep special marked nodes as singleton n...
Definition: MueLu_OnePtAggregationAlgorithm_decl.hpp:87
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:332
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:191
MueLu::Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:71
MueLu::Aggregates
Container class for aggregation information.
Definition: MueLu_Aggregates_decl.hpp:100
MueLu::ONEPT
Definition: MueLu_Types.hpp:76
MueLu::AggregationPhase1Algorithm
Algorithm for coarsening a graph with uncoupled aggregation.
Definition: MueLu_AggregationPhase1Algorithm_decl.hpp:85
MueLu::UncoupledAggregationFactory::Build
void Build(Level &currentLevel) const
Build aggregates.
Definition: MueLu_UncoupledAggregationFactory_def.hpp:162
MueLu::AggregationPhase3Algorithm
Handle leftover nodes. Try to avoid singleton nodes.
Definition: MueLu_AggregationPhase3Algorithm_decl.hpp:78
MueLu::AggregationPhase2aAlgorithm
Among unaggregated points, see if we can make a reasonable size aggregate out of it.
Definition: MueLu_AggregationPhase2aAlgorithm_decl.hpp:84
MueLu::AggregationPhase2bAlgorithm
Add leftovers to existing aggregates.
Definition: MueLu_AggregationPhase2bAlgorithm_decl.hpp:81
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:84
MueLu_GraphBase.hpp
MueLu_MasterList.hpp
MueLu::READY
Definition: MueLu_Types.hpp:65
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:59
MueLu::BOUNDARY
Definition: MueLu_Types.hpp:82
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99