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<CrsGraph>& myGraph, RCP<const Map>& coarseCoordinatesFineMap,
134 RCP<const Map>& coarseCoordinatesMap)
const {
135 Monitor m(*
this,
"BuildGraphP");
137 RCP<Teuchos::FancyOStream> out;
138 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
139 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
140 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
142 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
145 const bool coupled = geoData->isAggregationCoupled();
148 int numInterpolationPoints = 0;
149 if(geoData->getInterpolationOrder() == 0) {
150 numInterpolationPoints = 1;
151 }
else if(geoData->getInterpolationOrder() == 1) {
153 numInterpolationPoints = 1 << geoData->getNumDimensions();
155 *out <<
"numInterpolationPoints=" << numInterpolationPoints << std::endl;
157 Array<LO> colIndex((geoData->getNumLocalCoarseNodes() + numInterpolationPoints*
158 (geoData->getNumLocalFineNodes() - geoData->getNumLocalCoarseNodes()))*dofsPerNode);
159 Array<size_t> rowPtr(geoData->getNumLocalFineNodes()*dofsPerNode + 1);
161 ArrayRCP<size_t> nnzOnRow(geoData->getNumLocalFineNodes()*dofsPerNode);
163 *out <<
"Compute prolongatorGraph data" << std::endl;
164 if(geoData->getInterpolationOrder() == 0) {
165 ComputeGraphDataConstant(graph, geoData, dofsPerNode, numInterpolationPoints,
166 nnzOnRow, rowPtr, colIndex);
167 }
else if(geoData->getInterpolationOrder() == 1) {
168 ComputeGraphDataLinear(graph, geoData, dofsPerNode, numInterpolationPoints,
169 nnzOnRow, rowPtr, colIndex);
173 RCP<Map> rowMap = MapFactory::Build(graph.
GetDomainMap(), dofsPerNode);
174 RCP<Map> colMap, domainMap;
175 *out <<
"Compute domain and column maps of the CrsGraph" << std::endl;
177 *out <<
"Extract data for ghosted nodes" << std::endl;
178 Array<LO> ghostedCoarseNodeCoarseLIDs;
179 Array<int> ghostedCoarseNodeCoarsePIDs;
180 Array<GO> ghostedCoarseNodeCoarseGIDs;
181 geoData->getGhostedNodesData(graph.
GetDomainMap(), ghostedCoarseNodeCoarseLIDs,
182 ghostedCoarseNodeCoarsePIDs, ghostedCoarseNodeCoarseGIDs);
187 geoData->getNumGlobalCoarseNodes(),
188 ghostedCoarseNodeCoarseGIDs(),
192 LO coarseNodeIdx = 0;
193 Array<GO> coarseNodeCoarseGIDs, coarseNodeFineGIDs;
194 geoData->getCoarseNodesData(graph.
GetDomainMap(), coarseNodeCoarseGIDs, coarseNodeFineGIDs);
195 for(LO nodeIdx = 0; nodeIdx < ghostedCoarseNodeCoarseGIDs.size(); ++nodeIdx) {
196 if(ghostedCoarseNodeCoarsePIDs[nodeIdx] == colMap->getComm()->getRank()) {
197 coarseNodeCoarseGIDs[coarseNodeIdx] = ghostedCoarseNodeCoarseGIDs[nodeIdx];
201 domainMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
202 geoData->getNumGlobalCoarseNodes(),
203 coarseNodeCoarseGIDs(),
206 coarseCoordinatesMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
207 geoData->getNumGlobalCoarseNodes(),
208 coarseNodeCoarseGIDs(),
211 coarseCoordinatesFineMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
212 geoData->getNumGlobalCoarseNodes(),
213 coarseNodeFineGIDs(),
220 Teuchos::OrdinalTraits<GO>::invalid(),
221 geoData->getNumLocalCoarseNodes()*dofsPerNode,
226 Array<GO> coarseNodeCoarseGIDs(geoData->getNumLocalCoarseNodes());
227 Array<GO> coarseNodeFineGIDs(geoData->getNumLocalCoarseNodes());
228 geoData->getCoarseNodesData(graph.
GetDomainMap(), coarseNodeCoarseGIDs, coarseNodeFineGIDs);
229 coarseCoordinatesMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
230 Teuchos::OrdinalTraits<GO>::invalid(),
231 geoData->getNumLocalCoarseNodes(),
234 coarseCoordinatesFineMap = MapFactory::Build(graph.
GetDomainMap()->lib(),
235 Teuchos::OrdinalTraits<GO>::invalid(),
236 coarseNodeFineGIDs(),
241 *out <<
"Call constructor of CrsGraph" << std::endl;
242 myGraph = CrsGraphFactory::Build(rowMap,
246 *out <<
"Fill CrsGraph" << std::endl;
248 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
249 for(LO dof = 0; dof < dofsPerNode; ++dof) {
250 rowIdx = nodeIdx*dofsPerNode + dof;
251 myGraph->insertLocalIndices(rowIdx, colIndex(rowPtr[rowIdx], nnzOnRow[rowIdx]) );
255 *out <<
"Call fillComplete on CrsGraph" << std::endl;
256 myGraph->fillComplete(domainMap, rowMap);
257 *out <<
"Prolongator CrsGraph computed" << std::endl;
262 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
265 const LO dofsPerNode,
const int ,
266 ArrayRCP<size_t>& nnzOnRow, Array<size_t>& rowPtr,
267 Array<LO>& colIndex)
const {
269 RCP<Teuchos::FancyOStream> out;
270 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
271 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
272 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
274 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
277 Array<LO> ghostedCoarseNodeCoarseLIDs;
278 Array<int> ghostedCoarseNodeCoarsePIDs;
279 Array<GO> ghostedCoarseNodeCoarseGIDs;
280 geoData->getGhostedNodesData(graph.
GetDomainMap(), ghostedCoarseNodeCoarseLIDs,
281 ghostedCoarseNodeCoarsePIDs, ghostedCoarseNodeCoarseGIDs);
283 LO ghostedCoarseNodeCoarseLID, rem, rate;
284 Array<LO> ghostedIdx(3), coarseIdx(3);
285 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
288 geoData->getFineNodeGhostedTuple(nodeIdx, ghostedIdx[0], ghostedIdx[1], ghostedIdx[2]);
290 for(
int dim = 0; dim < 3; ++dim) {
291 coarseIdx[dim] = ghostedIdx[dim] / geoData->getCoarseningRate(dim);
292 rem = ghostedIdx[dim] % geoData->getCoarseningRate(dim);
293 if(ghostedIdx[dim] - geoData->getOffset(dim)
294 < geoData->getLocalFineNodesInDir(dim) - geoData->getCoarseningEndRate(dim)) {
295 rate = geoData->getCoarseningRate(dim);
297 rate = geoData->getCoarseningEndRate(dim);
299 if(rem > (rate / 2)) {++coarseIdx[dim];}
300 if( (geoData->getStartGhostedCoarseNode(dim)*geoData->getCoarseningRate(dim)
301 > geoData->getStartIndex(dim)) && geoData->isAggregationCoupled() ) {--coarseIdx[dim];}
304 geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2],
305 ghostedCoarseNodeCoarseLID);
307 for(LO dof = 0; dof < dofsPerNode; ++dof) {
308 nnzOnRow[nodeIdx*dofsPerNode + dof] = 1;
309 rowPtr[nodeIdx*dofsPerNode + dof + 1] = rowPtr[nodeIdx*dofsPerNode + dof] + 1;
310 colIndex[rowPtr[nodeIdx*dofsPerNode + dof]] =
311 ghostedCoarseNodeCoarseLIDs[ghostedCoarseNodeCoarseLID]*dofsPerNode + dof;
318 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
321 const LO dofsPerNode,
const int numInterpolationPoints,
322 ArrayRCP<size_t>& nnzOnRow, Array<size_t>& rowPtr,
323 Array<LO>& colIndex)
const {
325 RCP<Teuchos::FancyOStream> out;
326 if(
const char* dbg = std::getenv(
"MUELU_STRUCTUREDALGORITHM_DEBUG")) {
327 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
328 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
330 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
333 const bool coupled = geoData->isAggregationCoupled();
334 const int numDimensions = geoData->getNumDimensions();
335 Array<LO> ghostedIdx(3,0);
336 Array<LO> coarseIdx(3,0);
337 Array<LO> ijkRem(3,0);
338 const LO coarsePointOffset[8][3] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0},
339 {0, 0, 1}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1}};
341 for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
344 geoData->getFineNodeGhostedTuple(nodeIdx, ghostedIdx[0], ghostedIdx[1], ghostedIdx[2]);
345 for(
int dim=0; dim < numDimensions; dim++){
346 coarseIdx[dim] = ghostedIdx[dim] / geoData->getCoarseningRate(dim);
347 ijkRem[dim] = ghostedIdx[dim] % geoData->getCoarseningRate(dim);
349 if (geoData->getStartGhostedCoarseNode(dim)*geoData->getCoarseningRate(dim)
350 > geoData->getStartIndex(dim)) {
354 if(ghostedIdx[dim] == geoData->getLocalFineNodesInDir(dim) - 1) {
355 coarseIdx[dim] = geoData->getLocalCoarseNodesInDir(dim) - 1;
362 bool allCoarse =
true;
363 Array<bool> isCoarse(numDimensions);
364 for(
int dim = 0; dim < numDimensions; ++dim) {
365 isCoarse[dim] =
false;
367 isCoarse[dim] =
true;
370 if( ghostedIdx[dim]-geoData->getOffset(dim) == geoData->getLocalFineNodesInDir(dim)-1 &&
371 geoData->getMeshEdge(dim*2+1) )
372 isCoarse[dim] =
true;
374 if( ghostedIdx[dim]-geoData->getOffset(dim) == geoData->getLocalFineNodesInDir(dim)-1)
375 isCoarse[dim] =
true;
382 LO rowIdx = 0, colIdx = 0;
384 for(LO dof = 0; dof < dofsPerNode; ++dof) {
385 rowIdx = nodeIdx*dofsPerNode + dof;
386 nnzOnRow[rowIdx] = 1;
387 rowPtr[rowIdx + 1] = rowPtr[rowIdx] + 1;
390 geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2], colIdx);
391 colIndex[rowPtr[rowIdx]] = colIdx*dofsPerNode + dof;
395 for(
int dim = 0; dim < numDimensions; ++dim) {
396 if(coarseIdx[dim] == geoData->getGhostedNodesInDir(dim) - 1)
400 for(LO dof = 0; dof < dofsPerNode; ++dof) {
402 rowIdx = nodeIdx*dofsPerNode + dof;
403 nnzOnRow[rowIdx] = Teuchos::as<size_t>( numInterpolationPoints );
404 rowPtr[rowIdx + 1] = rowPtr[rowIdx] + Teuchos::as<LO>(numInterpolationPoints);
406 for(LO interpIdx = 0; interpIdx < numInterpolationPoints; ++interpIdx) {
407 geoData->getCoarseNodeGhostedLID(coarseIdx[0] + coarsePointOffset[interpIdx][0],
408 coarseIdx[1] + coarsePointOffset[interpIdx][1],
409 coarseIdx[2] + coarsePointOffset[interpIdx][2],
411 colIndex[rowPtr[rowIdx] + interpIdx] = colIdx*dofsPerNode + dof;
Container class for aggregation information.
RCP< IndexManager > & GetIndexManager()
Get the index manager used by structured aggregation algorithms.
virtual const RCP< const Map > GetDomainMap() const =0
Namespace for MueLu classes and methods.
void ComputeGraphDataLinear(const GraphBase &graph, RCP< IndexManager > &geoData, const LO dofsPerNode, const int numInterpolationPoints, ArrayRCP< size_t > &nnzOnRow, Array< size_t > &rowPtr, Array< LO > &colIndex) const
void BuildAggregates(const Teuchos::ParameterList ¶ms, const GraphBase &graph, Aggregates &aggregates, std::vector< unsigned > &aggStat, LO &numNonAggregatedNodes) const
Local aggregation.
void ComputeGraphDataConstant(const GraphBase &graph, RCP< IndexManager > &geoData, const LO dofsPerNode, const int numInterpolationPoints, ArrayRCP< size_t > &nnzOnRow, Array< size_t > &rowPtr, Array< LO > &colIndex) const
MueLu representation of a graph.
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, const LO dofsPerNode, RCP< CrsGraph > &myGraph, RCP< const Map > &coarseCoordinatesFineMap, RCP< const Map > &coarseCoordinatesMap) const
Local aggregation.