46 #ifndef MUELU_AGGREGATIONSTRUCTUREDALGORITHM_DEF_HPP_ 47 #define MUELU_AGGREGATIONSTRUCTUREDALGORITHM_DEF_HPP_ 50 #include <Teuchos_Comm.hpp> 51 #include <Teuchos_CommHelpers.hpp> 53 #include <Xpetra_MapFactory.hpp> 54 #include <Xpetra_Map.hpp> 55 #include <Xpetra_CrsGraphFactory.hpp> 56 #include <Xpetra_CrsGraph.hpp> 61 #include "MueLu_Aggregates.hpp" 62 #include "MueLu_IndexManager.hpp" 68 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 Aggregates& aggregates, std::vector<unsigned>& aggStat,
72 LO& numNonAggregatedNodes)
const {
73 Monitor m(*
this,
"BuildAggregates");
75 RCP<Teuchos::FancyOStream> out;
76 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
77 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
78 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
80 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
84 const bool coupled = geoData->isAggregationCoupled();
85 ArrayRCP<LO> vertex2AggId = aggregates.
GetVertex2AggId()->getDataNonConst(0);
86 ArrayRCP<LO> procWinner = aggregates.
GetProcWinner() ->getDataNonConst(0);
87 Array<LO> ghostedCoarseNodeCoarseLIDs;
88 Array<int> ghostedCoarseNodeCoarsePIDs;
89 Array<GO> ghostedCoarseNodeCoarseGIDs;
91 *out <<
"Extract data for ghosted nodes" << std::endl;
92 geoData->getGhostedNodesData(graph.
GetDomainMap(), ghostedCoarseNodeCoarseLIDs,
93 ghostedCoarseNodeCoarsePIDs, ghostedCoarseNodeCoarseGIDs);
96 Array<LO> ghostedIdx(3), coarseIdx(3);
97 LO ghostedCoarseNodeCoarseLID, aggId;
98 *out <<
"Loop over fine nodes and assign them to an aggregate and a rank" << std::endl;
99 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
101 geoData->getFineNodeGhostedTuple(nodeIdx, ghostedIdx[0], ghostedIdx[1], ghostedIdx[2]);
103 for(
int dim = 0; dim < 3; ++dim) {
104 coarseIdx[dim] = ghostedIdx[dim] / geoData->getCoarseningRate(dim);
105 rem = ghostedIdx[dim] % geoData->getCoarseningRate(dim);
106 if(ghostedIdx[dim] - geoData->getOffset(dim)
107 < geoData->getLocalFineNodesInDir(dim) - geoData->getCoarseningEndRate(dim)) {
108 rate = geoData->getCoarseningRate(dim);
110 rate = geoData->getCoarseningEndRate(dim);
112 if(rem > (rate / 2)) {++coarseIdx[dim];}
113 if(coupled && (geoData->getStartGhostedCoarseNode(dim)*geoData->getCoarseningRate(dim)
114 > geoData->getStartIndex(dim))) {--coarseIdx[dim];}
117 geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2],
118 ghostedCoarseNodeCoarseLID);
120 aggId = ghostedCoarseNodeCoarseLIDs[ghostedCoarseNodeCoarseLID];
121 vertex2AggId[nodeIdx] = aggId;
122 procWinner[nodeIdx] = ghostedCoarseNodeCoarsePIDs[ghostedCoarseNodeCoarseLID];
124 --numNonAggregatedNodes;
130 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
133 RCP<const Map>& coarseCoordinatesFineMap, RCP<const Map>& coarseCoordinatesMap)
const {
134 Monitor m(*
this,
"BuildAggregates");
136 RCP<Teuchos::FancyOStream> out;
137 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
138 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
139 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
141 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
144 const bool coupled = geoData->isAggregationCoupled();
147 int numInterpolationPoints = 0;
148 if(geoData->getInterpolationOrder() == 0) {
149 numInterpolationPoints = 1;
150 }
else if(geoData->getInterpolationOrder() == 1) {
152 numInterpolationPoints = 1 << geoData->getNumDimensions();
154 *out <<
"numInterpolationPoints=" << numInterpolationPoints << std::endl;
156 Array<LO> colIndex( geoData->getNumLocalCoarseNodes() + numInterpolationPoints*
157 (geoData->getNumLocalFineNodes() - geoData->getNumLocalCoarseNodes()) );
158 Array<size_t> rowPtr(geoData->getNumLocalFineNodes()+1);
160 ArrayRCP<size_t> nnzOnRow(geoData->getNumLocalFineNodes());
162 *out <<
"Compute prolongatorGraph data" << std::endl;
163 if(geoData->getInterpolationOrder() == 0) {
164 ComputeGraphDataConstant(graph, geoData, numInterpolationPoints, nnzOnRow, rowPtr, colIndex);
165 }
else if(geoData->getInterpolationOrder() == 1) {
166 ComputeGraphDataLinear(graph, geoData, numInterpolationPoints, nnzOnRow, rowPtr, colIndex);
170 RCP<Map> colMap, domainMap;
171 *out <<
"Compute domain and column maps of the CrsGraph" << std::endl;
173 *out <<
"Extract data for ghosted nodes" << std::endl;
174 Array<LO> ghostedCoarseNodeCoarseLIDs;
175 Array<int> ghostedCoarseNodeCoarsePIDs;
176 Array<GO> ghostedCoarseNodeCoarseGIDs;
177 geoData->getGhostedNodesData(graph.
GetDomainMap(), ghostedCoarseNodeCoarseLIDs,
178 ghostedCoarseNodeCoarsePIDs, ghostedCoarseNodeCoarseGIDs);
183 geoData->getNumGlobalCoarseNodes(),
184 ghostedCoarseNodeCoarseGIDs(),
189 LO coarseNodeIdx = 0;
190 Array<GO> coarseNodeCoarseGIDs, coarseNodeFineGIDs;
191 geoData->getCoarseNodesData(graph.
GetDomainMap(), coarseNodeCoarseGIDs, coarseNodeFineGIDs);
192 for(LO nodeIdx = 0; nodeIdx < ghostedCoarseNodeCoarseGIDs.size(); ++nodeIdx) {
193 if(ghostedCoarseNodeCoarsePIDs[nodeIdx] == colMap->getComm()->getRank()) {
194 coarseNodeCoarseGIDs[coarseNodeIdx] = ghostedCoarseNodeCoarseGIDs[nodeIdx];
198 domainMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
199 geoData->getNumGlobalCoarseNodes(),
200 coarseNodeCoarseGIDs(),
204 coarseCoordinatesMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
205 geoData->getNumGlobalCoarseNodes(),
206 coarseNodeCoarseGIDs(),
210 coarseCoordinatesFineMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
211 geoData->getNumGlobalCoarseNodes(),
212 coarseNodeFineGIDs(),
221 geoData->getNumGlobalCoarseNodes(),
222 geoData->getNumLocalCoarseNodes(),
228 Array<GO> coarseNodeCoarseGIDs(geoData->getNumLocalCoarseNodes());
229 Array<GO> coarseNodeFineGIDs(geoData->getNumLocalCoarseNodes());
230 geoData->getCoarseNodesData(graph.
GetDomainMap(), coarseNodeCoarseGIDs, coarseNodeFineGIDs);
231 coarseCoordinatesMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
232 geoData->getNumGlobalCoarseNodes(),
233 geoData->getNumLocalCoarseNodes(),
237 coarseCoordinatesFineMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
238 geoData->getNumGlobalCoarseNodes(),
239 coarseNodeFineGIDs(),
248 Xpetra::DynamicProfile);
249 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
250 myGraph->insertLocalIndices(nodeIdx, colIndex(rowPtr[nodeIdx], nnzOnRow[nodeIdx]) );
257 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
260 const int numInterpolationPoints, ArrayRCP<size_t>& nnzOnRow,
261 Array<size_t>& rowPtr, Array<LO>& colIndex)
const {
263 RCP<Teuchos::FancyOStream> out;
264 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
265 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
266 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
268 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
271 Array<LO> ghostedCoarseNodeCoarseLIDs;
272 Array<int> ghostedCoarseNodeCoarsePIDs;
273 Array<GO> ghostedCoarseNodeCoarseGIDs;
274 geoData->getGhostedNodesData(graph.
GetDomainMap(), ghostedCoarseNodeCoarseLIDs,
275 ghostedCoarseNodeCoarsePIDs, ghostedCoarseNodeCoarseGIDs);
277 LO ghostedCoarseNodeCoarseLID, rem, rate;
278 Array<LO> ghostedIdx(3), coarseIdx(3);
279 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
281 nnzOnRow[nodeIdx] = Teuchos::as<size_t>(1);
282 rowPtr[nodeIdx + 1] = rowPtr[nodeIdx] + 1;
285 geoData->getFineNodeGhostedTuple(nodeIdx, ghostedIdx[0], ghostedIdx[1], ghostedIdx[2]);
287 for(
int dim = 0; dim < 3; ++dim) {
288 coarseIdx[dim] = ghostedIdx[dim] / geoData->getCoarseningRate(dim);
289 rem = ghostedIdx[dim] % geoData->getCoarseningRate(dim);
290 if(ghostedIdx[dim] - geoData->getOffset(dim)
291 < geoData->getLocalFineNodesInDir(dim) - geoData->getCoarseningEndRate(dim)) {
292 rate = geoData->getCoarseningRate(dim);
294 rate = geoData->getCoarseningEndRate(dim);
296 if(rem > (rate / 2)) {++coarseIdx[dim];}
297 if( (geoData->getStartGhostedCoarseNode(dim)*geoData->getCoarseningRate(dim)
298 > geoData->getStartIndex(dim)) && geoData->isAggregationCoupled() ) {--coarseIdx[dim];}
301 geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2],
302 ghostedCoarseNodeCoarseLID);
303 colIndex[rowPtr[nodeIdx]] = ghostedCoarseNodeCoarseLIDs[ghostedCoarseNodeCoarseLID];
309 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
312 const int numInterpolationPoints, ArrayRCP<size_t>& nnzOnRow,
313 Array<size_t>& rowPtr, Array<LO>& colIndex)
const {
315 RCP<Teuchos::FancyOStream> out;
316 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
317 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
318 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
320 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
323 const bool coupled = geoData->isAggregationCoupled();
324 const int numDimensions = geoData->getNumDimensions();
325 Array<LO> ghostedIdx(3,0);
326 Array<LO> coarseIdx(3,0);
327 Array<LO> ijkRem(3,0);
330 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
333 geoData->getFineNodeGhostedTuple(nodeIdx, ghostedIdx[0], ghostedIdx[1], ghostedIdx[2]);
334 for(
int dim=0; dim < numDimensions; dim++){
335 coarseIdx[dim] = ghostedIdx[dim] / geoData->getCoarseningRate(dim);
336 ijkRem[dim] = ghostedIdx[dim] % geoData->getCoarseningRate(dim);
337 if(ghostedIdx[dim] - geoData->getOffset(dim)
338 < geoData->getLocalFineNodesInDir(dim) - geoData->getCoarseningEndRate(dim)) {
339 rate = geoData->getCoarseningRate(dim);
341 rate = geoData->getCoarseningEndRate(dim);
343 if(ijkRem[dim] > (rate / 2)) {++coarseIdx[dim];}
344 if(coupled && (geoData->getStartGhostedCoarseNode(dim)*geoData->getCoarseningRate(dim)
345 > geoData->getStartIndex(dim))) {--coarseIdx[dim];}
350 bool allCoarse =
true;
351 Array<bool> isCoarse(numDimensions);
352 for(
int dim = 0; dim < numDimensions; ++dim) {
353 isCoarse[dim] =
false;
355 isCoarse[dim] =
true;
358 if( ghostedIdx[dim]-geoData->getOffset(dim) == geoData->getLocalFineNodesInDir(dim)-1 &&
359 geoData->getMeshEdge(dim*2+1) )
360 isCoarse[dim] =
true;
362 if( ghostedIdx[dim]-geoData->getOffset(dim) == geoData->getLocalFineNodesInDir(dim)-1)
363 isCoarse[dim] =
true;
372 geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2],
373 colIndex[rowPtr[nodeIdx]]);
374 nnzOnRow[nodeIdx] = Teuchos::as<size_t>(1);
375 rowPtr[nodeIdx + 1] = rowPtr[nodeIdx] + 1;
379 nnzOnRow[nodeIdx] = Teuchos::as<size_t>( numInterpolationPoints );
380 rowPtr[nodeIdx + 1] = rowPtr[nodeIdx] + Teuchos::as<LO>( numInterpolationPoints );
382 for(
int dim = 0; dim < numDimensions; ++dim) {
383 if(coarseIdx[dim] == geoData->getGhostedNodesInDir(dim) - 1)
387 geoData->getCoarseNodeGhostedLID( coarseIdx[0], coarseIdx[1], coarseIdx[2], colIndex[ rowPtr[nodeIdx]+0]);
388 geoData->getCoarseNodeGhostedLID( coarseIdx[0]+1, coarseIdx[1], coarseIdx[2], colIndex[ rowPtr[nodeIdx]+1]);
389 if(numDimensions > 1) {
390 geoData->getCoarseNodeGhostedLID( coarseIdx[0], coarseIdx[1]+1, coarseIdx[2], colIndex[ rowPtr[nodeIdx]+2]);
391 geoData->getCoarseNodeGhostedLID( coarseIdx[0]+1, coarseIdx[1]+1, coarseIdx[2], colIndex[ rowPtr[nodeIdx]+3]);
392 if(numDimensions > 2) {
393 geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2]+1, colIndex[ rowPtr[nodeIdx]+4]);
394 geoData->getCoarseNodeGhostedLID(coarseIdx[0]+1, coarseIdx[1], coarseIdx[2]+1, colIndex[ rowPtr[nodeIdx]+5]);
395 geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1]+1, coarseIdx[2]+1, colIndex[ rowPtr[nodeIdx]+6]);
396 geoData->getCoarseNodeGhostedLID(coarseIdx[0]+1, coarseIdx[1]+1, coarseIdx[2]+1, colIndex[ rowPtr[nodeIdx]+7]);
Container class for aggregation information.
RCP< IndexManager > & GetIndexManager()
Get the index manager used by structured aggregation algorithms.
void ComputeGraphDataConstant(const GraphBase &graph, RCP< IndexManager > &geoData, const int numInterpolationPoints, ArrayRCP< size_t > &nnzOnRow, Array< size_t > &rowPtr, Array< LO > &colIndex) const
virtual const RCP< const Map > GetDomainMap() const =0
Namespace for MueLu classes and methods.
void BuildAggregates(const Teuchos::ParameterList ¶ms, const GraphBase &graph, Aggregates &aggregates, std::vector< unsigned > &aggStat, LO &numNonAggregatedNodes) const
Local aggregation.
MueLu representation of a graph.
void ComputeGraphDataLinear(const GraphBase &graph, RCP< IndexManager > &geoData, const int numInterpolationPoints, ArrayRCP< size_t > &nnzOnRow, Array< size_t > &rowPtr, Array< LO > &colIndex) const
Timer to be used in non-factories.
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
const RCP< LOVector > & GetProcWinner() const
Returns constant vector that maps local node IDs to owning processor IDs.
void BuildGraph(const GraphBase &graph, RCP< IndexManager > &geoData, RCP< CrsGraph > &myGraph, RCP< const Map > &coarseCoordinatesFineMap, RCP< const Map > &coarseCoordinatesMap) const
Local aggregation.