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_VectorFactory.hpp>
54 
56 
57 #include "MueLu_OnePtAggregationAlgorithm.hpp"
58 #include "MueLu_PreserveDirichletAggregationAlgorithm.hpp"
59 #include "MueLu_IsolatedNodeAggregationAlgorithm.hpp"
60 
61 #include "MueLu_AggregationPhase1Algorithm.hpp"
62 #include "MueLu_AggregationPhase2aAlgorithm.hpp"
63 #include "MueLu_AggregationPhase2bAlgorithm.hpp"
64 #include "MueLu_AggregationPhase3Algorithm.hpp"
65 
66 #include "MueLu_Level.hpp"
67 #include "MueLu_GraphBase.hpp"
68 #include "MueLu_Aggregates.hpp"
69 #include "MueLu_MasterList.hpp"
70 #include "MueLu_Monitor.hpp"
71 #include "MueLu_AmalgamationInfo.hpp"
72 #include "MueLu_Utilities.hpp"
73 
74 namespace MueLu {
75 
76  template <class LocalOrdinal, class GlobalOrdinal, class Node>
78  : bDefinitionPhase_(true)
79  { }
80 
81  template <class LocalOrdinal, class GlobalOrdinal, class Node>
83  RCP<ParameterList> validParamList = rcp(new ParameterList());
84 
85  // Aggregation parameters (used in aggregation algorithms)
86  // TODO introduce local member function for each aggregation algorithm such that each aggregation algorithm can define its own parameters
87 
89 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
90  SET_VALID_ENTRY("aggregation: max agg size");
91  SET_VALID_ENTRY("aggregation: min agg size");
92  SET_VALID_ENTRY("aggregation: max selected neighbors");
93  SET_VALID_ENTRY("aggregation: ordering");
94  validParamList->getEntry("aggregation: ordering").setValidator(
95  rcp(new validatorType(Teuchos::tuple<std::string>("natural", "graph", "random"), "aggregation: ordering")));
96  SET_VALID_ENTRY("aggregation: enable phase 1");
97  SET_VALID_ENTRY("aggregation: enable phase 2a");
98  SET_VALID_ENTRY("aggregation: enable phase 2b");
99  SET_VALID_ENTRY("aggregation: enable phase 3");
100  SET_VALID_ENTRY("aggregation: preserve Dirichlet points");
101  SET_VALID_ENTRY("aggregation: allow user-specified singletons");
102 #undef SET_VALID_ENTRY
103 
104  // general variables needed in AggregationFactory
105  validParamList->set< RCP<const FactoryBase> >("Graph", null, "Generating factory of the graph");
106  validParamList->set< RCP<const FactoryBase> >("DofsPerNode", null, "Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
107 
108  // special variables necessary for OnePtAggregationAlgorithm
109  validParamList->set< std::string > ("OnePt aggregate map name", "", "Name of input map for single node aggregates. (default='')");
110  validParamList->set< std::string > ("OnePt aggregate map factory", "", "Generating factory of (DOF) map for single node aggregates.");
111  //validParamList->set< RCP<const FactoryBase> >("OnePt aggregate map factory", NoFactory::getRCP(), "Generating factory of (DOF) map for single node aggregates.");
112 
113  return validParamList;
114  }
115 
116  template <class LocalOrdinal, class GlobalOrdinal, class Node>
118  Input(currentLevel, "Graph");
119  Input(currentLevel, "DofsPerNode");
120 
121  const ParameterList& pL = GetParameterList();
122 
123  // request special data necessary for OnePtAggregationAlgorithm
124  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
125  if (mapOnePtName.length() > 0) {
126  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
127  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
128  currentLevel.DeclareInput(mapOnePtName, NoFactory::get());
129  } else {
130  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
131  currentLevel.DeclareInput(mapOnePtName, mapOnePtFact.get());
132  }
133  }
134  }
135 
136  template <class LocalOrdinal, class GlobalOrdinal, class Node>
138  FactoryMonitor m(*this, "Build", currentLevel);
139 
140  ParameterList pL = GetParameterList();
141  bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
142 
143  if (pL.get<int>("aggregation: max agg size") == -1)
144  pL.set("aggregation: max agg size", INT_MAX);
145 
146  // define aggregation algorithms
147  RCP<const FactoryBase> graphFact = GetFactory("Graph");
148 
149  // TODO Can we keep different aggregation algorithms over more Build calls?
150  algos_.clear();
151  algos_.push_back(rcp(new PreserveDirichletAggregationAlgorithm(graphFact)));
152  if (pL.get<bool>("aggregation: allow user-specified singletons") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm (graphFact)));
153  if (pL.get<bool>("aggregation: enable phase 1" ) == true) algos_.push_back(rcp(new AggregationPhase1Algorithm (graphFact)));
154  if (pL.get<bool>("aggregation: enable phase 2a") == true) algos_.push_back(rcp(new AggregationPhase2aAlgorithm (graphFact)));
155  if (pL.get<bool>("aggregation: enable phase 2b") == true) algos_.push_back(rcp(new AggregationPhase2bAlgorithm (graphFact)));
156  if (pL.get<bool>("aggregation: enable phase 3" ) == true) algos_.push_back(rcp(new AggregationPhase3Algorithm (graphFact)));
157 
158  // TODO: remove old aggregation mode
159  //if (pL.get<bool>("UseOnePtAggregationAlgorithm") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm (graphFact)));
160  //if (pL.get<bool>("UseUncoupledAggregationAlgorithm") == true) algos_.push_back(rcp(new AggregationPhase1Algorithm (graphFact)));
161  //if (pL.get<bool>("UseMaxLinkAggregationAlgorithm") == true) algos_.push_back(rcp(new MaxLinkAggregationAlgorithm (graphFact)));
162  //if (pL.get<bool>("UseEmergencyAggregationAlgorithm") == true) algos_.push_back(rcp(new EmergencyAggregationAlgorithm (graphFact)));
163 
164  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
165  RCP<Map> OnePtMap = Teuchos::null;
166  if (mapOnePtName.length()) {
167  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
168  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
169  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, NoFactory::get());
170  } else {
171  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
172  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
173  }
174  }
175 
176  RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel, "Graph");
177 
178  // Build
179  RCP<Aggregates> aggregates = rcp(new Aggregates(*graph));
180  aggregates->setObjectLabel("UC");
181 
182  const LO numRows = graph->GetNodeNumVertices();
183 
184  // construct aggStat information
185  std::vector<unsigned> aggStat(numRows, READY);
186 
187  ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
188  if (dirichletBoundaryMap != Teuchos::null)
189  for (LO i = 0; i < numRows; i++)
190  if (dirichletBoundaryMap[i] == true)
191  aggStat[i] = BOUNDARY;
192 
193  LO nDofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
194  GO indexBase = graph->GetDomainMap()->getIndexBase();
195  if (OnePtMap != Teuchos::null) {
196  for (LO i = 0; i < numRows; i++) {
197  // reconstruct global row id (FIXME only works for contiguous maps)
198  GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
199 
200  for (LO kr = 0; kr < nDofsPerNode; kr++)
201  if (OnePtMap->isNodeGlobalElement(grid + kr))
202  aggStat[i] = ONEPT;
203  }
204  }
205 
206 
207  const RCP<const Teuchos::Comm<int> > comm = graph->GetComm();
208  GO numGlobalRows = 0;
209  if (IsPrint(Statistics1))
210  MueLu_sumAll(comm, as<GO>(numRows), numGlobalRows);
211 
212  LO numNonAggregatedNodes = numRows;
213  GO numGlobalAggregatedPrev = 0, numGlobalAggsPrev = 0;
214  for (size_t a = 0; a < algos_.size(); a++) {
215  std::string phase = algos_[a]->description();
216  SubFactoryMonitor sfm(*this, "Algo \"" + phase + "\"", currentLevel);
217 
218  int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
219  algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
220  algos_[a]->SetProcRankVerbose(oldRank);
221 
222  if (IsPrint(Statistics1)) {
223  GO numLocalAggregated = numRows - numNonAggregatedNodes, numGlobalAggregated = 0;
224  GO numLocalAggs = aggregates->GetNumAggregates(), numGlobalAggs = 0;
225  MueLu_sumAll(comm, numLocalAggregated, numGlobalAggregated);
226  MueLu_sumAll(comm, numLocalAggs, numGlobalAggs);
227 
228  double aggPercent = 100*as<double>(numGlobalAggregated)/as<double>(numGlobalRows);
229  if (aggPercent > 99.99 && aggPercent < 100.00) {
230  // Due to round off (for instance, for 140465733/140466897), we could
231  // get 100.00% display even if there are some remaining nodes. This
232  // is bad from the users point of view. It is much better to change
233  // it to display 99.99%.
234  aggPercent = 99.99;
235  }
236  GetOStream(Statistics1) << " aggregated : " << (numGlobalAggregated - numGlobalAggregatedPrev) << " (phase), " << std::fixed
237  << std::setprecision(2) << numGlobalAggregated << "/" << numGlobalRows << " [" << aggPercent << "%] (total)\n"
238  << " remaining : " << numGlobalRows - numGlobalAggregated << "\n"
239  << " aggregates : " << numGlobalAggs-numGlobalAggsPrev << " (phase), " << numGlobalAggs << " (total)" << std::endl;
240  numGlobalAggregatedPrev = numGlobalAggregated;
241  numGlobalAggsPrev = numGlobalAggs;
242  }
243  }
244 
245  TEUCHOS_TEST_FOR_EXCEPTION(numNonAggregatedNodes, Exceptions::RuntimeError, "MueLu::UncoupledAggregationFactory::Build: Leftover nodes found! Error!");
246 
247  aggregates->AggregatesCrossProcessors(false);
248 
249  Set(currentLevel, "Aggregates", aggregates);
250 
251  GetOStream(Statistics1) << aggregates->description() << std::endl;
252  }
253 
254 } //namespace MueLu
255 
256 
257 #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.
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
void DeclareInput(Level &currentLevel) const
Input.
Namespace for MueLu classes and methods.
void Build(Level &currentLevel) const
Build aggregates.
T * get() const
static const NoFactory * get()
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.
Add leftovers to existing aggregatesIdeaIn phase 2b non-aggregated nodes are added to existing aggreg...
Algorithm for coarsening a graph with uncoupled aggregation.
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()