46#ifndef MUELU_HYBRIDAGGREGATIONFACTORY_DEF_HPP_
47#define MUELU_HYBRIDAGGREGATIONFACTORY_DEF_HPP_
49#include <Xpetra_Matrix.hpp>
50#include <Xpetra_Map.hpp>
51#include <Xpetra_Vector.hpp>
52#include <Xpetra_MultiVectorFactory.hpp>
53#include <Xpetra_VectorFactory.hpp>
58#include "MueLu_InterfaceAggregationAlgorithm.hpp"
59#include "MueLu_OnePtAggregationAlgorithm.hpp"
60#include "MueLu_PreserveDirichletAggregationAlgorithm.hpp"
61#include "MueLu_IsolatedNodeAggregationAlgorithm.hpp"
63#include "MueLu_AggregationPhase1Algorithm.hpp"
64#include "MueLu_AggregationPhase2aAlgorithm.hpp"
65#include "MueLu_AggregationPhase2bAlgorithm.hpp"
66#include "MueLu_AggregationPhase3Algorithm.hpp"
69#include "MueLu_AggregationStructuredAlgorithm.hpp"
70#include "MueLu_UncoupledIndexManager.hpp"
77#include "MueLu_Aggregates.hpp"
80#include "MueLu_Utilities.hpp"
81#include "MueLu_AmalgamationInfo.hpp"
86 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
91 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 RCP<ParameterList> validParamList = rcp(
new ParameterList());
96 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
97#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
103 validParamList->getEntry(
"aggregation: ordering").setValidator(
104 rcp(
new validatorType(Teuchos::tuple<std::string>(
"natural",
"graph",
"random"),
"aggregation: ordering")));
111 SET_VALID_ENTRY(
"aggregation: error on nodes with no on-rank neighbors");
123#undef SET_VALID_ENTRY
127 validParamList->set< RCP<const FactoryBase> >(
"Graph", null,
"Generating factory of the graph");
128 validParamList->set< RCP<const FactoryBase> >(
"DofsPerNode", null,
"Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
130 validParamList->set<std::string> (
"OnePt aggregate map name",
"",
131 "Name of input map for single node aggregates. (default='')");
132 validParamList->set<std::string> (
"OnePt aggregate map factory",
"",
133 "Generating factory of (DOF) map for single node aggregates.");
136 validParamList->set<std::string> (
"Interface aggregate map name",
"",
137 "Name of input map for interface aggregates. (default='')");
138 validParamList->set<std::string> (
"Interface aggregate map factory",
"",
139 "Generating factory of (DOF) map for interface aggregates.");
140 validParamList->set<RCP<const FactoryBase> > (
"interfacesDimensions", Teuchos::null,
141 "Describes the dimensions of all the interfaces on this rank.");
142 validParamList->set<RCP<const FactoryBase> > (
"nodeOnInterface", Teuchos::null,
143 "List the LIDs of the nodes on any interface.");
147 validParamList->set<RCP<const FactoryBase> >(
"numDimensions", Teuchos::null,
148 "Number of spatial dimension provided by CoordinatesTransferFactory.");
149 validParamList->set<RCP<const FactoryBase> >(
"lNodesPerDim", Teuchos::null,
150 "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
154 validParamList->set<RCP<const FactoryBase> > (
"aggregationRegionType", Teuchos::null,
155 "Type of aggregation to use on the region (\"structured\" or \"uncoupled\")");
157 return validParamList;
160 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
163 Input(currentLevel,
"Graph");
165 ParameterList pL = GetParameterList();
178 "Aggregation region type was not provided by the user!");
185 "numDimensions was not provided by the user on level0!");
192 "lNodesPerDim was not provided by the user on level0!");
195 Input(currentLevel,
"aggregationRegionType");
196 Input(currentLevel,
"numDimensions");
197 Input(currentLevel,
"lNodesPerDim");
203 Input(currentLevel,
"DofsPerNode");
206 if (pL.get<
bool>(
"aggregation: use interface aggregation") ==
true){
213 "interfacesDimensions was not provided by the user on level0!");
220 "nodeOnInterface was not provided by the user on level0!");
223 Input(currentLevel,
"interfacesDimensions");
224 Input(currentLevel,
"nodeOnInterface");
229 std::string mapOnePtName = pL.get<std::string>(
"OnePt aggregate map name");
230 if (mapOnePtName.length() > 0) {
231 std::string mapOnePtFactName = pL.get<std::string>(
"OnePt aggregate map factory");
232 if (mapOnePtFactName ==
"" || mapOnePtFactName ==
"NoFactory") {
235 RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
236 currentLevel.
DeclareInput(mapOnePtName, mapOnePtFact.get());
241 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
246 RCP<Teuchos::FancyOStream> out;
247 if(
const char* dbg = std::getenv(
"MUELU_HYBRIDAGGREGATION_DEBUG")) {
248 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
249 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
251 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
254 *out <<
"Entering hybrid aggregation" << std::endl;
256 ParameterList pL = GetParameterList();
257 bDefinitionPhase_ =
false;
259 if (pL.get<
int>(
"aggregation: max agg size") == -1)
260 pL.set(
"aggregation: max agg size", INT_MAX);
263 RCP<const FactoryBase> graphFact = GetFactory(
"Graph");
266 RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel,
"Graph");
267 RCP<const Map> fineMap = graph->GetDomainMap();
268 const int myRank = fineMap->getComm()->getRank();
269 const int numRanks = fineMap->getComm()->getSize();
271 out->setProcRankAndSize(graph->GetImportMap()->getComm()->getRank(),
272 graph->GetImportMap()->getComm()->getSize());
275 RCP<Aggregates> aggregates = rcp(
new Aggregates(*graph));
276 aggregates->setObjectLabel(
"HB");
279 const LO numRows = graph->GetNodeNumVertices();
280 std::vector<unsigned> aggStat(numRows,
READY);
283 std::string regionType;
286 regionType = currentLevel.
Get<std::string>(
"aggregationRegionType",
NoFactory::get());
289 regionType = Get< std::string >(currentLevel,
"aggregationRegionType");
292 int numDimensions = 0;
298 numDimensions = Get<int>(currentLevel,
"numDimensions");
302 std::string coarseningRate = pL.get<std::string>(
"aggregation: coarsening rate");
303 Teuchos::Array<LO> coarseRate;
305 coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
306 }
catch(
const Teuchos::InvalidArrayStringRepresentation& e) {
307 GetOStream(
Errors,-1) <<
" *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** "
311 TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
313 "\"aggregation: coarsening rate\" must have at least as many"
314 " components as the number of spatial dimensions in the problem.");
317 LO numNonAggregatedNodes = numRows;
318 if (regionType ==
"structured") {
324 const int interpolationOrder = pL.get<
int>(
"aggregation: coarsening order");
325 Array<LO> lFineNodesPerDir(3);
328 lFineNodesPerDir = currentLevel.
Get<Array<LO> >(
"lNodesPerDim",
NoFactory::get());
331 lFineNodesPerDir = Get<Array<LO> >(currentLevel,
"lNodesPerDim");
335 for(
int dim = numDimensions; dim < 3; ++dim) {
336 lFineNodesPerDir[dim] = 1;
340 RCP<MueLu::IndexManager<LO,GO,NO> > geoData;
351 TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getLocalNumElements()
352 !=
static_cast<size_t>(geoData->getNumLocalFineNodes()),
354 "The local number of elements in the graph's map is not equal to "
355 "the number of nodes given by: lNodesPerDim!");
357 aggregates->SetIndexManager(geoData);
358 aggregates->SetNumAggregates(geoData->getNumLocalCoarseNodes());
360 Set(currentLevel,
"lCoarseNodesPerDim", geoData->getLocalCoarseNodesPerDir());
364 if (regionType ==
"uncoupled"){
368 if (pL.get<
bool>(
"aggregation: allow user-specified singletons") ==
true) algos_.push_back(rcp(
new OnePtAggregationAlgorithm (graphFact)));
374 *out <<
" Build interface aggregates" << std::endl;
376 if (pL.get<
bool>(
"aggregation: use interface aggregation") ==
true) {
377 BuildInterfaceAggregates(currentLevel, aggregates, aggStat, numNonAggregatedNodes,
381 *out <<
"Treat Dirichlet BC" << std::endl;
383 ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
384 if (dirichletBoundaryMap != Teuchos::null)
385 for (LO i = 0; i < numRows; i++)
386 if (dirichletBoundaryMap[i] ==
true)
390 std::string mapOnePtName = pL.get<std::string>(
"OnePt aggregate map name");
391 RCP<Map> OnePtMap = Teuchos::null;
392 if (mapOnePtName.length()) {
393 std::string mapOnePtFactName = pL.get<std::string>(
"OnePt aggregate map factory");
394 if (mapOnePtFactName ==
"" || mapOnePtFactName ==
"NoFactory") {
397 RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
398 OnePtMap = currentLevel.
Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
402 LO nDofsPerNode = Get<LO>(currentLevel,
"DofsPerNode");
403 GO indexBase = graph->GetDomainMap()->getIndexBase();
404 if (OnePtMap != Teuchos::null) {
405 for (LO i = 0; i < numRows; i++) {
407 GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
408 for (LO kr = 0; kr < nDofsPerNode; kr++)
409 if (OnePtMap->isNodeGlobalElement(grid + kr))
415 Array<LO> lCoarseNodesPerDir(3,-1);
416 Set(currentLevel,
"lCoarseNodesPerDim", lCoarseNodesPerDir);
419 aggregates->AggregatesCrossProcessors(
false);
421 *out <<
"Run all the algorithms on the local rank" << std::endl;
422 for (
size_t a = 0; a < algos_.size(); a++) {
423 std::string phase = algos_[a]->description();
425 *out << regionType <<
" | Executing phase " << a << std::endl;
427 int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
428 algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
429 algos_[a]->SetProcRankVerbose(oldRank);
430 *out << regionType <<
" | Done Executing phase " << a << std::endl;
433 *out <<
"Compute statistics on aggregates" << std::endl;
434 aggregates->ComputeAggregateSizes(
true);
436 Set(currentLevel,
"Aggregates", aggregates);
437 Set(currentLevel,
"numDimensions", numDimensions);
438 Set(currentLevel,
"aggregationRegionTypeCoarse", regionType);
440 GetOStream(
Statistics1) << aggregates->description() << std::endl;
441 *out <<
"HybridAggregation done!" << std::endl;
444 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
447 std::vector<unsigned>& aggStat, LO& numNonAggregatedNodes,
448 Array<LO> coarseRate)
const {
449 FactoryMonitor m(*
this,
"BuildInterfaceAggregates", currentLevel);
451 RCP<Teuchos::FancyOStream> out;
452 if(
const char* dbg = std::getenv(
"MUELU_HYBRIDAGGREGATION_DEBUG")) {
453 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
454 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
456 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
460 if(coarseRate.size() == 1) {coarseRate.resize(3, coarseRate[0]);}
461 ArrayRCP<LO> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
462 ArrayRCP<LO> procWinner = aggregates->GetProcWinner() ->getDataNonConst(0);
463 Array<LO> interfacesDimensions = Get<Array<LO> >(currentLevel,
"interfacesDimensions");
464 Array<LO> nodesOnInterfaces = Get<Array<LO> >(currentLevel,
"nodeOnInterface");
465 const int numInterfaces = interfacesDimensions.size() / 3;
466 const int myRank = aggregates->GetMap()->getComm()->getRank();
469 Array<LO> coarseInterfacesDimensions(interfacesDimensions.size());
470 Array<LO> nodesOnCoarseInterfaces;
472 LO endRate, totalNumCoarseNodes = 0, numCoarseNodes;
473 for(
int interfaceIdx = 0; interfaceIdx < numInterfaces; ++interfaceIdx) {
475 for(
int dim = 0; dim < 3; ++dim) {
476 endRate = (interfacesDimensions[3*interfaceIdx + dim] - 1) % coarseRate[dim];
477 if(interfacesDimensions[3*interfaceIdx + dim] == 1) {
478 coarseInterfacesDimensions[3*interfaceIdx + dim] = 1;
480 coarseInterfacesDimensions[3*interfaceIdx + dim]
481 = (interfacesDimensions[3*interfaceIdx+dim]-1) / coarseRate[dim] + 2;
482 if(endRate==0){ coarseInterfacesDimensions[3*interfaceIdx + dim]--;}
484 numCoarseNodes *= coarseInterfacesDimensions[3*interfaceIdx + dim];
486 totalNumCoarseNodes += numCoarseNodes;
488 nodesOnCoarseInterfaces.resize(totalNumCoarseNodes, -1);
491 Array<LO> endRate(3);
492 LO interfaceOffset = 0, aggregateCount = 0, coarseNodeCount = 0;
493 for(
int interfaceIdx = 0; interfaceIdx < numInterfaces; ++interfaceIdx) {
494 ArrayView<LO> fineNodesPerDim = interfacesDimensions(3*interfaceIdx, 3);
495 ArrayView<LO> coarseNodesPerDim = coarseInterfacesDimensions(3*interfaceIdx, 3);
496 LO numInterfaceNodes = 1, numCoarseNodes = 1;
497 for(
int dim = 0; dim < 3; ++dim) {
498 numInterfaceNodes *= fineNodesPerDim[dim];
499 numCoarseNodes *= coarseNodesPerDim[dim];
500 endRate[dim] = (fineNodesPerDim[dim]-1) % coarseRate[dim];
502 ArrayView<LO> interfaceNodes = nodesOnInterfaces(interfaceOffset, numInterfaceNodes);
504 interfaceOffset += numInterfaceNodes;
506 LO rem, rate, fineNodeIdx;
507 Array<LO> nodeIJK(3), coarseIJK(3), rootIJK(3);
510 for(LO coarseNodeIdx = 0; coarseNodeIdx < numCoarseNodes; ++coarseNodeIdx) {
511 coarseIJK[2] = coarseNodeIdx / (coarseNodesPerDim[0]*coarseNodesPerDim[1]);
512 rem = coarseNodeIdx % (coarseNodesPerDim[0]*coarseNodesPerDim[1]);
513 coarseIJK[1] = rem / coarseNodesPerDim[0];
514 coarseIJK[0] = rem % coarseNodesPerDim[0];
516 for(LO dim = 0; dim < 3; ++dim) {
517 if(coarseIJK[dim] == coarseNodesPerDim[dim] - 1) {
518 nodeIJK[dim] = fineNodesPerDim[dim] - 1;
520 nodeIJK[dim] = coarseIJK[dim]*coarseRate[dim];
523 fineNodeIdx = (nodeIJK[2]*fineNodesPerDim[1] + nodeIJK[1])*fineNodesPerDim[0] + nodeIJK[0];
525 if(aggStat[interfaceNodes[fineNodeIdx]] ==
READY) {
526 vertex2AggId[interfaceNodes[fineNodeIdx]] = aggregateCount;
527 procWinner[interfaceNodes[fineNodeIdx]] = myRank;
528 aggStat[interfaceNodes[fineNodeIdx]] =
AGGREGATED;
530 --numNonAggregatedNodes;
532 nodesOnCoarseInterfaces[coarseNodeCount] = vertex2AggId[interfaceNodes[fineNodeIdx]];
539 for(LO nodeIdx = 0; nodeIdx < numInterfaceNodes; ++nodeIdx) {
542 if(aggStat[interfaceNodes[nodeIdx]] ==
AGGREGATED) {
continue;}
544 nodeIJK[2] = nodeIdx / (fineNodesPerDim[0]*fineNodesPerDim[1]);
545 rem = nodeIdx % (fineNodesPerDim[0]*fineNodesPerDim[1]);
546 nodeIJK[1] = rem / fineNodesPerDim[0];
547 nodeIJK[0] = rem % fineNodesPerDim[0];
549 for(
int dim = 0; dim < 3; ++dim) {
550 coarseIJK[dim] = nodeIJK[dim] / coarseRate[dim];
551 rem = nodeIJK[dim] % coarseRate[dim];
552 if(nodeIJK[dim] < fineNodesPerDim[dim] - endRate[dim]) {
553 rate = coarseRate[dim];
557 if(rem > (rate / 2)) {++coarseIJK[dim];}
560 for(LO dim = 0; dim < 3; ++dim) {
561 if(coarseIJK[dim] == coarseNodesPerDim[dim] - 1) {
562 nodeIJK[dim] = fineNodesPerDim[dim] - 1;
564 nodeIJK[dim] = coarseIJK[dim]*coarseRate[dim];
567 fineNodeIdx = (nodeIJK[2]*fineNodesPerDim[1] + nodeIJK[1])*fineNodesPerDim[0] + nodeIJK[0];
569 vertex2AggId[interfaceNodes[nodeIdx]] = vertex2AggId[interfaceNodes[fineNodeIdx]];
570 procWinner[interfaceNodes[nodeIdx]] = myRank;
571 aggStat[interfaceNodes[nodeIdx]] =
AGGREGATED;
572 --numNonAggregatedNodes;
577 aggregates->SetNumAggregates(aggregateCount);
580 Set(currentLevel,
"coarseInterfacesDimensions", coarseInterfacesDimensions);
581 Set(currentLevel,
"nodeOnCoarseInterface", nodesOnCoarseInterfaces);
#define SET_VALID_ENTRY(name)
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.
Add leftovers to existing aggregates.
Handle leftover nodes. Try to avoid singleton nodes.
Algorithm for coarsening a graph with structured aggregation.
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.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level ¤tLevel) const
Input.
HybridAggregationFactory()
Constructor.
void Build(Level ¤tLevel) const
Build aggregates.
void BuildInterfaceAggregates(Level ¤tLevel, RCP< Aggregates > aggregates, std::vector< unsigned > &aggStat, LO &numNonAggregatedNodes, Array< LO > coarseRate) const
Specifically build aggregates along interfaces.
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.
Namespace for MueLu classes and methods.
@ Statistics1
Print more statistics.