MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
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 
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>
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()),
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>
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_ */
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-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
Container class for aggregation information.
void setValidator(RCP< const ParameterEntryValidator > const &validator)
GlobalOrdinal GO
T & get(const std::string &name, T def_value)
virtual size_t GetNodeNumVertices() const =0
Return number of vertices owned by the calling node.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
virtual const ArrayRCP< const bool > GetBoundaryNodeMap() const =0
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Print more statistics.
virtual const RCP< const Map > GetDomainMap() const =0
LocalOrdinal LO
T * get() const
void DeclareInput(Level &currentLevel) const
Input.
static const NoFactory * get()
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...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
virtual const RCP< const Teuchos::Comm< int > > GetComm() const =0
#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.
void Build(Level &currentLevel) const
Build aggregates.
Add leftovers to existing aggregatesIdeaIn phase 2b non-aggregated nodes are added to existing aggreg...
Algorithm for coarsening a graph with uncoupled aggregation.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
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 ...
ParameterEntry & getEntry(const std::string &name)
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.