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