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: 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>
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()),
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>
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_ */
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()