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 
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>
80  : bDefinitionPhase_(true)
81  { }
82 
83  template <class LocalOrdinal, class GlobalOrdinal, class Node>
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: phase2a include root");
103  SET_VALID_ENTRY("aggregation: phase2a agg factor");
104  SET_VALID_ENTRY("aggregation: preserve Dirichlet points");
105  SET_VALID_ENTRY("aggregation: allow user-specified singletons");
106  SET_VALID_ENTRY("aggregation: use interface aggregation");
107  SET_VALID_ENTRY("aggregation: error on nodes with no on-rank neighbors");
108  SET_VALID_ENTRY("aggregation: phase3 avoid singletons");
109  SET_VALID_ENTRY("aggregation: compute aggregate qualities");
110 #undef SET_VALID_ENTRY
111 
112  // general variables needed in AggregationFactory
113  validParamList->set< RCP<const FactoryBase> >("Graph", null, "Generating factory of the graph");
114  validParamList->set< RCP<const FactoryBase> >("DofsPerNode", null, "Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
115  validParamList->set< RCP<const FactoryBase> >("AggregateQualities", null, "Generating factory for variable \'AggregateQualities\'");
116 
117  // special variables necessary for OnePtAggregationAlgorithm
118  validParamList->set< std::string > ("OnePt aggregate map name", "", "Name of input map for single node aggregates. (default='')");
119  validParamList->set< std::string > ("OnePt aggregate map factory", "", "Generating factory of (DOF) map for single node aggregates.");
120  //validParamList->set< RCP<const FactoryBase> >("OnePt aggregate map factory", NoFactory::getRCP(), "Generating factory of (DOF) map for single node aggregates.");
121 
122  // InterfaceAggregation parameters
123  //validParamList->set< bool > ("aggregation: use interface aggregation", "false", "Flag to trigger aggregation along an interface using specified aggregate seeds.");
124  validParamList->set< std::string > ("Interface aggregate map name", "", "Name of input map for interface aggregates. (default='')");
125  validParamList->set< std::string > ("Interface aggregate map factory", "", "Generating factory of (DOF) map for interface aggregates.");
126  validParamList->set<RCP<const FactoryBase> > ("nodeOnInterface", Teuchos::null, "Array specifying whether or not a node is on the interface (1 or 0).");
127 
128  return validParamList;
129  }
130 
131  template <class LocalOrdinal, class GlobalOrdinal, class Node>
133  Input(currentLevel, "Graph");
134  Input(currentLevel, "DofsPerNode");
135 
136  const ParameterList& pL = GetParameterList();
137 
138  // request special data necessary for OnePtAggregationAlgorithm
139  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
140  if (mapOnePtName.length() > 0) {
141  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
142  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
143  currentLevel.DeclareInput(mapOnePtName, NoFactory::get());
144  } else {
145  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
146  currentLevel.DeclareInput(mapOnePtName, mapOnePtFact.get());
147  }
148  }
149 
150  // request special data necessary for InterfaceAggregation
151  if (pL.get<bool>("aggregation: use interface aggregation") == true){
152  if(currentLevel.GetLevelID() == 0) {
153  if(currentLevel.IsAvailable("nodeOnInterface", NoFactory::get())) {
154  currentLevel.DeclareInput("nodeOnInterface", NoFactory::get(), this);
155  } else {
156  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("nodeOnInterface", NoFactory::get()),
158  "nodeOnInterface was not provided by the user on level0!");
159  }
160  } else {
161  Input(currentLevel, "nodeOnInterface");
162  }
163  }
164 
165  if (pL.get<bool>("aggregation: compute aggregate qualities")) {
166  Input(currentLevel, "AggregateQualities");
167  }
168  }
169 
170  template <class LocalOrdinal, class GlobalOrdinal, class Node>
172  FactoryMonitor m(*this, "Build", currentLevel);
173 
174  ParameterList pL = GetParameterList();
175  bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
176 
177  if (pL.get<int>("aggregation: max agg size") == -1)
178  pL.set("aggregation: max agg size", INT_MAX);
179 
180  // define aggregation algorithms
181  RCP<const FactoryBase> graphFact = GetFactory("Graph");
182 
183  // TODO Can we keep different aggregation algorithms over more Build calls?
184  algos_.clear();
185  algos_.push_back(rcp(new PreserveDirichletAggregationAlgorithm(graphFact)));
186  if (pL.get<bool>("aggregation: use interface aggregation") == true) algos_.push_back(rcp(new InterfaceAggregationAlgorithm (graphFact)));
187  if (pL.get<bool>("aggregation: allow user-specified singletons") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm (graphFact)));
188  if (pL.get<bool>("aggregation: enable phase 1" ) == true) algos_.push_back(rcp(new AggregationPhase1Algorithm (graphFact)));
189  if (pL.get<bool>("aggregation: enable phase 2a") == true) algos_.push_back(rcp(new AggregationPhase2aAlgorithm (graphFact)));
190  if (pL.get<bool>("aggregation: enable phase 2b") == true) algos_.push_back(rcp(new AggregationPhase2bAlgorithm (graphFact)));
191  if (pL.get<bool>("aggregation: enable phase 3" ) == true) algos_.push_back(rcp(new AggregationPhase3Algorithm (graphFact)));
192 
193  // TODO: remove old aggregation mode
194  //if (pL.get<bool>("UseOnePtAggregationAlgorithm") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm (graphFact)));
195  //if (pL.get<bool>("UseUncoupledAggregationAlgorithm") == true) algos_.push_back(rcp(new AggregationPhase1Algorithm (graphFact)));
196  //if (pL.get<bool>("UseMaxLinkAggregationAlgorithm") == true) algos_.push_back(rcp(new MaxLinkAggregationAlgorithm (graphFact)));
197  //if (pL.get<bool>("UseEmergencyAggregationAlgorithm") == true) algos_.push_back(rcp(new EmergencyAggregationAlgorithm (graphFact)));
198 
199  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
200  RCP<Map> OnePtMap = Teuchos::null;
201  if (mapOnePtName.length()) {
202  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
203  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
204  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, NoFactory::get());
205  } else {
206  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
207  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
208  }
209  }
210 
211  // Set map for interface aggregates
212  std::string mapInterfaceName = pL.get<std::string>("Interface aggregate map name");
213  RCP<Map> InterfaceMap = Teuchos::null;
214 
215  RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel, "Graph");
216 
217  // Build
218  RCP<Aggregates> aggregates = rcp(new Aggregates(*graph));
219  aggregates->setObjectLabel("UC");
220 
221  const LO numRows = graph->GetNodeNumVertices();
222 
223  // construct aggStat information
224  std::vector<unsigned> aggStat(numRows, READY);
225 
226  // interface
227  if (pL.get<bool>("aggregation: use interface aggregation") == true){
228  Teuchos::Array<LO> nodeOnInterface = Get<Array<LO>>(currentLevel,"nodeOnInterface");
229  for (LO i = 0; i < numRows; i++) {
230  if (nodeOnInterface[i])
231  aggStat[i] = INTERFACE;
232  }
233  }
234 
235  ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
236  if (dirichletBoundaryMap != Teuchos::null)
237  for (LO i = 0; i < numRows; i++)
238  if (dirichletBoundaryMap[i] == true)
239  aggStat[i] = BOUNDARY;
240 
241  LO nDofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
242  GO indexBase = graph->GetDomainMap()->getIndexBase();
243  if (OnePtMap != Teuchos::null) {
244  for (LO i = 0; i < numRows; i++) {
245  // reconstruct global row id (FIXME only works for contiguous maps)
246  GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
247 
248  for (LO kr = 0; kr < nDofsPerNode; kr++)
249  if (OnePtMap->isNodeGlobalElement(grid + kr))
250  aggStat[i] = ONEPT;
251  }
252  }
253 
254 
255 
256  const RCP<const Teuchos::Comm<int> > comm = graph->GetComm();
257  GO numGlobalRows = 0;
258  if (IsPrint(Statistics1))
259  MueLu_sumAll(comm, as<GO>(numRows), numGlobalRows);
260 
261  LO numNonAggregatedNodes = numRows;
262  GO numGlobalAggregatedPrev = 0, numGlobalAggsPrev = 0;
263  for (size_t a = 0; a < algos_.size(); a++) {
264  std::string phase = algos_[a]->description();
265  SubFactoryMonitor sfm(*this, "Algo " + phase, currentLevel);
266 
267  int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
268  algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
269  algos_[a]->SetProcRankVerbose(oldRank);
270 
271  if (IsPrint(Statistics1)) {
272  GO numLocalAggregated = numRows - numNonAggregatedNodes, numGlobalAggregated = 0;
273  GO numLocalAggs = aggregates->GetNumAggregates(), numGlobalAggs = 0;
274  MueLu_sumAll(comm, numLocalAggregated, numGlobalAggregated);
275  MueLu_sumAll(comm, numLocalAggs, numGlobalAggs);
276 
277  double aggPercent = 100*as<double>(numGlobalAggregated)/as<double>(numGlobalRows);
278  if (aggPercent > 99.99 && aggPercent < 100.00) {
279  // Due to round off (for instance, for 140465733/140466897), we could
280  // get 100.00% display even if there are some remaining nodes. This
281  // is bad from the users point of view. It is much better to change
282  // it to display 99.99%.
283  aggPercent = 99.99;
284  }
285  GetOStream(Statistics1) << " aggregated : " << (numGlobalAggregated - numGlobalAggregatedPrev) << " (phase), " << std::fixed
286  << std::setprecision(2) << numGlobalAggregated << "/" << numGlobalRows << " [" << aggPercent << "%] (total)\n"
287  << " remaining : " << numGlobalRows - numGlobalAggregated << "\n"
288  << " aggregates : " << numGlobalAggs-numGlobalAggsPrev << " (phase), " << numGlobalAggs << " (total)" << std::endl;
289  numGlobalAggregatedPrev = numGlobalAggregated;
290  numGlobalAggsPrev = numGlobalAggs;
291  }
292  }
293 
294  TEUCHOS_TEST_FOR_EXCEPTION(numNonAggregatedNodes, Exceptions::RuntimeError, "MueLu::UncoupledAggregationFactory::Build: Leftover nodes found! Error!");
295 
296  aggregates->AggregatesCrossProcessors(false);
297  aggregates->ComputeAggregateSizes(true/*forceRecompute*/);
298 
299  Set(currentLevel, "Aggregates", aggregates);
300 
301  if (pL.get<bool>("aggregation: compute aggregate qualities")) {
302  RCP<Xpetra::MultiVector<double,LO,GO,Node>> aggQualities = Get<RCP<Xpetra::MultiVector<double,LO,GO,Node>>>(currentLevel, "AggregateQualities");
303  }
304 
305  if (IsPrint(Statistics1))
306  GetOStream(Statistics1) << aggregates->description() << std::endl;
307  }
308 
309 } //namespace MueLu
310 
311 
312 #endif /* MUELU_UNCOUPLEDAGGREGATIONFACTORY_DEF_HPP_ */
Algorithm for coarsening a graph with uncoupled aggregation. keep special marked nodes as singleton n...
#define MueLu_sumAll(rcpComm, in, out)
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
Container class for aggregation information.
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
void DeclareInput(Level &currentLevel) const
Input.
Namespace for MueLu classes and methods.
void Build(Level &currentLevel) const
Build aggregates.
static const NoFactory * get()
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
Algorithm for coarsening a graph with uncoupled aggregation. creates aggregates along an interface us...
Builds one-to-one aggregates for all Dirichlet boundary nodes. For some applications this might be ne...
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
#define SET_VALID_ENTRY(name)
Among unaggregated points, see if we can make a reasonable size aggregate out of it.IdeaAmong unaggregated points, see if we can make a reasonable size aggregate out of it. We do this by looking at neighbors and seeing how many are unaggregated and on my processor. Loosely, base the number of new aggregates created on the percentage of unaggregated nodes.
Add leftovers to existing aggregatesIdeaIn phase 2b non-aggregated nodes are added to existing aggreg...
Algorithm for coarsening a graph with uncoupled aggregation.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Exception throws to report errors in the internal logical of the program.
Handle leftover nodes. Try to avoid singleton nodesIdeaIn phase 3 we try to stick unaggregated nodes ...
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()