From b653e0e09b76507d40f3515f99567f96ad5cc0c8 Mon Sep 17 00:00:00 2001 From: Luc Berger-Vergiat Date: Fri, 11 Jan 2019 16:02:59 -0700 Subject: [PATCH] MueLu: modifying structured aggregation to compute graph of prolongator The modifies the prolongatorGraph output from structured aggregation. Instead of getting a node based graph it now returns a dofs based graph. This means that in the new format, the prolongator graph can be used to call the constructor of the prolongator matrix directly. This has the advantage of reducing some overhead and it also simplifies the implementation of the geometric interpolation factory. --- ...Lu_AggregationStructuredAlgorithm_decl.hpp | 12 +- ...eLu_AggregationStructuredAlgorithm_def.hpp | 111 +++++---- ...MueLu_StructuredAggregationFactory_def.hpp | 17 +- ...Lu_GeometricInterpolationPFactory_decl.hpp | 4 +- ...eLu_GeometricInterpolationPFactory_def.hpp | 228 +++++------------- .../muelu/test/structured/structured_2dof.xml | 5 +- .../muelu/test/structured/structured_3dof.xml | 1 + .../StructuredAggregationFactory.cpp | 10 + 8 files changed, 151 insertions(+), 237 deletions(-) diff --git a/packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_decl.hpp b/packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_decl.hpp index 87e20756b2ba..1889ce097b2c 100644 --- a/packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_decl.hpp +++ b/packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_decl.hpp @@ -106,7 +106,7 @@ namespace MueLu { /*! @brief Local aggregation. */ - void BuildGraph(const GraphBase& graph, RCP& geoData, + void BuildGraph(const GraphBase& graph, RCP& geoData, const LO dofsPerNode, RCP& myGraph, RCP& coarseCoordinatesFineMap, RCP& coarseCoordinatesMap) const; //@} @@ -116,12 +116,14 @@ namespace MueLu { private: void ComputeGraphDataConstant(const GraphBase& graph, RCP& geoData, - const int numInterpolationPoints, ArrayRCP& nnzOnRow, - Array& rowPtr, Array& colIndex) const; + const LO dofsPerNode, const int numInterpolationPoints, + ArrayRCP& nnzOnRow, Array& rowPtr, + Array& colIndex) const; void ComputeGraphDataLinear(const GraphBase& graph, RCP& geoData, - const int numInterpolationPoints, ArrayRCP& nnzOnRow, - Array& rowPtr, Array& colIndex) const; + const LO dofsPerNode, const int numInterpolationPoints, + ArrayRCP& nnzOnRow, Array& rowPtr, + Array& colIndex) const; }; diff --git a/packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_def.hpp b/packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_def.hpp index b960b192a2db..cbbb84c52c3c 100644 --- a/packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_def.hpp +++ b/packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_def.hpp @@ -129,8 +129,9 @@ namespace MueLu { template void AggregationStructuredAlgorithm:: - BuildGraph(const GraphBase& graph, RCP& geoData, RCP& myGraph, - RCP& coarseCoordinatesFineMap, RCP& coarseCoordinatesMap) const { + BuildGraph(const GraphBase& graph, RCP& geoData, const LO dofsPerNode, + RCP& myGraph, RCP& coarseCoordinatesFineMap, + RCP& coarseCoordinatesMap) const { Monitor m(*this, "BuildGraphP"); RCP out; @@ -153,20 +154,23 @@ namespace MueLu { } *out << "numInterpolationPoints=" << numInterpolationPoints << std::endl; - Array colIndex( geoData->getNumLocalCoarseNodes() + numInterpolationPoints* - (geoData->getNumLocalFineNodes() - geoData->getNumLocalCoarseNodes()) ); - Array rowPtr(geoData->getNumLocalFineNodes()+1); + Array colIndex((geoData->getNumLocalCoarseNodes() + numInterpolationPoints* + (geoData->getNumLocalFineNodes() - geoData->getNumLocalCoarseNodes()))*dofsPerNode); + Array rowPtr(geoData->getNumLocalFineNodes()*dofsPerNode + 1); rowPtr[0] = 0; - ArrayRCP nnzOnRow(geoData->getNumLocalFineNodes()); + ArrayRCP nnzOnRow(geoData->getNumLocalFineNodes()*dofsPerNode); *out << "Compute prolongatorGraph data" << std::endl; if(geoData->getInterpolationOrder() == 0) { - ComputeGraphDataConstant(graph, geoData, numInterpolationPoints, nnzOnRow, rowPtr, colIndex); + ComputeGraphDataConstant(graph, geoData, dofsPerNode, numInterpolationPoints, + nnzOnRow, rowPtr, colIndex); } else if(geoData->getInterpolationOrder() == 1) { - ComputeGraphDataLinear(graph, geoData, numInterpolationPoints, nnzOnRow, rowPtr, colIndex); + ComputeGraphDataLinear(graph, geoData, dofsPerNode, numInterpolationPoints, + nnzOnRow, rowPtr, colIndex); } - // Compute graph's colMap and domainMap + // Compute graph's rowMap, colMap and domainMap + RCP rowMap = MapFactory::Build(graph.GetDomainMap(), dofsPerNode); RCP colMap, domainMap; *out << "Compute domain and column maps of the CrsGraph" << std::endl; if(coupled){ @@ -215,11 +219,10 @@ namespace MueLu { graph.GetDomainMap()->getNode()); } else { // In this case the map will compute the global number of nodes on the coarse mesh - // since geoData->getNumGlobalCoarseNodes() == Teuchos::OrdinalTraits::invalid() // and it will assign GIDs to the local coarse nodes. colMap = MapFactory::Build(graph.GetDomainMap()->lib(), - geoData->getNumGlobalCoarseNodes(), - geoData->getNumLocalCoarseNodes(), + Teuchos::OrdinalTraits::invalid(), + geoData->getNumLocalCoarseNodes()*dofsPerNode, graph.GetDomainMap()->getIndexBase(), graph.GetDomainMap()->getComm(), graph.GetDomainMap()->getNode()); @@ -229,13 +232,13 @@ namespace MueLu { Array coarseNodeFineGIDs(geoData->getNumLocalCoarseNodes()); geoData->getCoarseNodesData(graph.GetDomainMap(), coarseNodeCoarseGIDs, coarseNodeFineGIDs); coarseCoordinatesMap = MapFactory::Build(graph.GetDomainMap()->lib(), - geoData->getNumGlobalCoarseNodes(), + Teuchos::OrdinalTraits::invalid(), geoData->getNumLocalCoarseNodes(), graph.GetDomainMap()->getIndexBase(), graph.GetDomainMap()->getComm(), graph.GetDomainMap()->getNode()); coarseCoordinatesFineMap = MapFactory::Build(graph.GetDomainMap()->lib(), - geoData->getNumGlobalCoarseNodes(), + Teuchos::OrdinalTraits::invalid(), coarseNodeFineGIDs(), graph.GetDomainMap()->getIndexBase(), graph.GetDomainMap()->getComm(), @@ -243,18 +246,22 @@ namespace MueLu { } *out << "Call constructor of CrsGraph" << std::endl; - myGraph = CrsGraphFactory::Build(graph.GetDomainMap(), + myGraph = CrsGraphFactory::Build(rowMap, colMap, nnzOnRow, Xpetra::DynamicProfile); *out << "Fill CrsGraph" << std::endl; + LO rowIdx = 0; for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) { - myGraph->insertLocalIndices(nodeIdx, colIndex(rowPtr[nodeIdx], nnzOnRow[nodeIdx]) ); + for(LO dof = 0; dof < dofsPerNode; ++dof) { + rowIdx = nodeIdx*dofsPerNode + dof; + myGraph->insertLocalIndices(rowIdx, colIndex(rowPtr[rowIdx], nnzOnRow[rowIdx]) ); + } } *out << "Call fillComplete on CrsGraph" << std::endl; - myGraph->fillComplete(domainMap, graph.GetDomainMap()); + myGraph->fillComplete(domainMap, rowMap); *out << "Prolongator CrsGraph computed" << std::endl; } // BuildAggregates() @@ -263,8 +270,9 @@ namespace MueLu { template void AggregationStructuredAlgorithm:: ComputeGraphDataConstant(const GraphBase& graph, RCP& geoData, - const int numInterpolationPoints, ArrayRCP& nnzOnRow, - Array& rowPtr, Array& colIndex) const { + const LO dofsPerNode, const int numInterpolationPoints, + ArrayRCP& nnzOnRow, Array& rowPtr, + Array& colIndex) const { RCP out; if(const char* dbg = std::getenv("MUELU_STRUCTUREDALGORITHM_DEBUG")) { @@ -283,9 +291,6 @@ namespace MueLu { LO ghostedCoarseNodeCoarseLID, rem, rate; Array ghostedIdx(3), coarseIdx(3); for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) { - // For piece-wise constant interpolation we only get one nnz per row - nnzOnRow[nodeIdx] = Teuchos::as(1); - rowPtr[nodeIdx + 1] = rowPtr[nodeIdx] + 1; // These needs to change for kokkos: rowPtr[nodeIdx + 1] = nodeIdx + 1; // Compute coarse ID associated with fine LID geoData->getFineNodeGhostedTuple(nodeIdx, ghostedIdx[0], ghostedIdx[1], ghostedIdx[2]); @@ -306,7 +311,13 @@ namespace MueLu { geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2], ghostedCoarseNodeCoarseLID); - colIndex[rowPtr[nodeIdx]] = ghostedCoarseNodeCoarseLIDs[ghostedCoarseNodeCoarseLID]; // Here too, substitute nodeIdx for rowPtr[nodeIdx] + + for(LO dof = 0; dof < dofsPerNode; ++dof) { + nnzOnRow[nodeIdx*dofsPerNode + dof] = 1; + rowPtr[nodeIdx*dofsPerNode + dof + 1] = rowPtr[nodeIdx*dofsPerNode + dof] + 1; + colIndex[rowPtr[nodeIdx*dofsPerNode + dof]] = + ghostedCoarseNodeCoarseLIDs[ghostedCoarseNodeCoarseLID]*dofsPerNode + dof; + } } // Loop over fine points } // ComputeGraphDataConstant() @@ -315,8 +326,9 @@ namespace MueLu { template void AggregationStructuredAlgorithm:: ComputeGraphDataLinear(const GraphBase& graph, RCP& geoData, - const int numInterpolationPoints, ArrayRCP& nnzOnRow, - Array& rowPtr, Array& colIndex) const { + const LO dofsPerNode, const int numInterpolationPoints, + ArrayRCP& nnzOnRow, Array& rowPtr, + Array& colIndex) const { RCP out; if(const char* dbg = std::getenv("MUELU_STRUCTUREDALGORITHM_DEBUG")) { @@ -332,6 +344,8 @@ namespace MueLu { Array coarseIdx(3,0); Array ijkRem(3,0); int rate = 0; + const LO coarsePointOffset[8][3] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0}, + {0, 0, 1}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1}}; for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) { @@ -373,35 +387,38 @@ namespace MueLu { allCoarse = false; } + LO rowIdx = 0, colIdx = 0; if(allCoarse) { - // Fine node lies on Coarse node, easy case, we only need the LID of the coarse node. - geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2], - colIndex[rowPtr[nodeIdx]]); - nnzOnRow[nodeIdx] = Teuchos::as(1); - rowPtr[nodeIdx + 1] = rowPtr[nodeIdx] + 1; + for(LO dof = 0; dof < dofsPerNode; ++dof) { + rowIdx = nodeIdx*dofsPerNode + dof; + nnzOnRow[rowIdx] = 1; + rowPtr[rowIdx + 1] = rowPtr[rowIdx] + 1; + + // Fine node lies on Coarse node, easy case, we only need the LID of the coarse node. + geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2], colIdx); + colIndex[rowPtr[rowIdx]] = colIdx*dofsPerNode + dof; + } } else { // Harder case, we need the LIDs of all the coarse nodes contributing to the interpolation - // at the current node. - nnzOnRow[nodeIdx] = Teuchos::as( numInterpolationPoints ); - rowPtr[nodeIdx + 1] = rowPtr[nodeIdx] + Teuchos::as( numInterpolationPoints ); - for(int dim = 0; dim < numDimensions; ++dim) { if(coarseIdx[dim] == geoData->getGhostedNodesInDir(dim) - 1) --coarseIdx[dim]; } - // Compute Coarse Node LID - geoData->getCoarseNodeGhostedLID( coarseIdx[0], coarseIdx[1], coarseIdx[2], colIndex[ rowPtr[nodeIdx]+0]); - geoData->getCoarseNodeGhostedLID( coarseIdx[0]+1, coarseIdx[1], coarseIdx[2], colIndex[ rowPtr[nodeIdx]+1]); - if(numDimensions > 1) { - geoData->getCoarseNodeGhostedLID( coarseIdx[0], coarseIdx[1]+1, coarseIdx[2], colIndex[ rowPtr[nodeIdx]+2]); - geoData->getCoarseNodeGhostedLID( coarseIdx[0]+1, coarseIdx[1]+1, coarseIdx[2], colIndex[ rowPtr[nodeIdx]+3]); - if(numDimensions > 2) { - geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1], coarseIdx[2]+1, colIndex[ rowPtr[nodeIdx]+4]); - geoData->getCoarseNodeGhostedLID(coarseIdx[0]+1, coarseIdx[1], coarseIdx[2]+1, colIndex[ rowPtr[nodeIdx]+5]); - geoData->getCoarseNodeGhostedLID(coarseIdx[0], coarseIdx[1]+1, coarseIdx[2]+1, colIndex[ rowPtr[nodeIdx]+6]); - geoData->getCoarseNodeGhostedLID(coarseIdx[0]+1, coarseIdx[1]+1, coarseIdx[2]+1, colIndex[ rowPtr[nodeIdx]+7]); - } - } + + for(LO dof = 0; dof < dofsPerNode; ++dof) { + // at the current node. + rowIdx = nodeIdx*dofsPerNode + dof; + nnzOnRow[rowIdx] = Teuchos::as( numInterpolationPoints ); + rowPtr[rowIdx + 1] = rowPtr[rowIdx] + Teuchos::as(numInterpolationPoints); + // Compute Coarse Node LID + for(LO interpIdx = 0; interpIdx < numInterpolationPoints; ++interpIdx) { + geoData->getCoarseNodeGhostedLID(coarseIdx[0] + coarsePointOffset[interpIdx][0], + coarseIdx[1] + coarsePointOffset[interpIdx][1], + coarseIdx[2] + coarsePointOffset[interpIdx][2], + colIdx); + colIndex[rowPtr[rowIdx] + interpIdx] = colIdx*dofsPerNode + dof; + } // Loop over numInterpolationPoints + } // Loop over dofsPerNode } } // Loop over fine points } // ComputeGraphDataLinear() diff --git a/packages/muelu/src/Graph/StructuredAggregation/MueLu_StructuredAggregationFactory_def.hpp b/packages/muelu/src/Graph/StructuredAggregation/MueLu_StructuredAggregationFactory_def.hpp index f9ea9b1e08ba..98c1373223a2 100644 --- a/packages/muelu/src/Graph/StructuredAggregation/MueLu_StructuredAggregationFactory_def.hpp +++ b/packages/muelu/src/Graph/StructuredAggregation/MueLu_StructuredAggregationFactory_def.hpp @@ -97,7 +97,6 @@ namespace MueLu { validParamList->set >("aggregation: mesh data", Teuchos::null, "Mesh ordering associated data"); - validParamList->set >("Graph", Teuchos::null, "Graph of the matrix after amalgamation but without dropping."); validParamList->set >("numDimensions", Teuchos::null, @@ -106,6 +105,8 @@ namespace MueLu { "Global number of nodes per spatial dimension provided by CoordinatesTransferFactory."); validParamList->set >("lNodesPerDim", Teuchos::null, "Local number of nodes per spatial dimension provided by CoordinatesTransferFactory."); + validParamList->set >("DofsPerNode", Teuchos::null, + "Generating factory for variable \'DofsPerNode\', usually the same as the \'Graph\' factory"); return validParamList; } // GetValidParameterList @@ -114,6 +115,7 @@ namespace MueLu { void StructuredAggregationFactory:: DeclareInput(Level& currentLevel) const { Input(currentLevel, "Graph"); + Input(currentLevel, "DofsPerNode"); ParameterList pL = GetParameterList(); std::string coupling = pL.get("aggregation: coupling"); @@ -175,10 +177,11 @@ namespace MueLu { // General problem informations are gathered from data stored in the problem matix. RCP graph = Get< RCP >(currentLevel, "Graph"); - RCP fineMap = graph->GetDomainMap(); - const int myRank = fineMap->getComm()->getRank(); - const int numRanks = fineMap->getComm()->getSize(); - const GO minGlobalIndex = fineMap->getMinGlobalIndex(); + RCP fineMap = graph->GetDomainMap(); + const int myRank = fineMap->getComm()->getRank(); + const int numRanks = fineMap->getComm()->getSize(); + const GO minGlobalIndex = fineMap->getMinGlobalIndex(); + const LO dofsPerNode = Get(currentLevel, "DofsPerNode"); // Since we want to operate on nodes and not dof, we need to modify the rowMap in order to // obtain a nodeMap. @@ -332,8 +335,8 @@ namespace MueLu { } else { // Create Coarse Data RCP myGraph; - myStructuredAlgorithm->BuildGraph(*graph, geoData, myGraph, coarseCoordinatesFineMap, - coarseCoordinatesMap); + myStructuredAlgorithm->BuildGraph(*graph, geoData, dofsPerNode, myGraph, + coarseCoordinatesFineMap, coarseCoordinatesMap); Set(currentLevel, "prolongatorGraph", myGraph); } diff --git a/packages/muelu/src/Transfers/GeneralGeometric/MueLu_GeometricInterpolationPFactory_decl.hpp b/packages/muelu/src/Transfers/GeneralGeometric/MueLu_GeometricInterpolationPFactory_decl.hpp index e71c394eb9a9..28ad8142aaca 100644 --- a/packages/muelu/src/Transfers/GeneralGeometric/MueLu_GeometricInterpolationPFactory_decl.hpp +++ b/packages/muelu/src/Transfers/GeneralGeometric/MueLu_GeometricInterpolationPFactory_decl.hpp @@ -98,8 +98,8 @@ namespace MueLu{ //@} private: - void BuildConstantP(RCP& P, RCP& prolongatorGraph, RCP& A) const; - void BuildLinearP(RCP& A, RCP& prolongatorGraph, + void BuildConstantP(RCP& P, RCP& prolongatorGraph, RCP& A) const; + void BuildLinearP(RCP& A, RCP& prolongatorGraph, RCP& fineCoordinates, RCP& ghostCoordinates, const int numDimensions, RCP& P) const; diff --git a/packages/muelu/src/Transfers/GeneralGeometric/MueLu_GeometricInterpolationPFactory_def.hpp b/packages/muelu/src/Transfers/GeneralGeometric/MueLu_GeometricInterpolationPFactory_def.hpp index a82554e5099a..b70b82217a55 100644 --- a/packages/muelu/src/Transfers/GeneralGeometric/MueLu_GeometricInterpolationPFactory_def.hpp +++ b/packages/muelu/src/Transfers/GeneralGeometric/MueLu_GeometricInterpolationPFactory_def.hpp @@ -138,7 +138,7 @@ namespace MueLu { // Declared main input/outputs to be retrieved and placed on the fine resp. coarse level RCP A = Get >(fineLevel, "A"); - RCP prolongatorGraph = Get >(fineLevel, "prolongatorGraph"); + RCP prolongatorGraph = Get >(fineLevel, "prolongatorGraph"); RCP fineCoordinates, coarseCoordinates; RCP P; const int interpolationOrder = pL.get("gmg: interpolation order"); @@ -178,7 +178,7 @@ namespace MueLu { BuildLinearP(A, prolongatorGraph, fineCoordinates, ghostCoordinates, numDimensions, P); } - *out << "The prolongator matrix has been build." << std::endl; + *out << "The prolongator matrix has been built." << std::endl; // Build the coarse nullspace RCP fineNullspace = Get< RCP > (fineLevel, "Nullspace"); @@ -201,7 +201,7 @@ namespace MueLu { template void GeometricInterpolationPFactory:: - BuildConstantP(RCP& P, RCP& prolongatorGraph, RCP& A) const { + BuildConstantP(RCP& P, RCP& prolongatorGraph, RCP& A) const { // Set debug outputs based on environment variable RCP out; @@ -214,90 +214,19 @@ namespace MueLu { *out << "BuildConstantP" << std::endl; - // Get the number of dofs per nodes from A and use that number to manually unamalgamate - // the maps of the prolongator graph. Eventually this operation will no longer be needed - // after the structured aggregation is refactored to handle this. - int dofsPerNode = A->GetFixedBlockSize(); - ArrayView initialDomainMapLIDs = - prolongatorGraph->getDomainMap()->getNodeElementList(); - Array domainMapLIDs(initialDomainMapLIDs.size()*dofsPerNode); - for(LO elementIdx = 0; elementIdx < as(initialDomainMapLIDs.size()); ++elementIdx) { - for(int dof = 0; dof < dofsPerNode; ++dof) { - domainMapLIDs[elementIdx*dofsPerNode + dof] = - initialDomainMapLIDs[elementIdx]*dofsPerNode + dof; - } - } - - *out << "Call domain map constructor" << std::endl; - - RCP domainMap = MapFactory::Build(prolongatorGraph->getRowMap()->lib(), - prolongatorGraph->getGlobalNumCols()*dofsPerNode, - domainMapLIDs(), - prolongatorGraph->getIndexBase(), - prolongatorGraph->getComm(), - prolongatorGraph->getRowMap()->getNode()); - - ArrayView initialColMapLIDs = - prolongatorGraph->getColMap()->getNodeElementList(); - Array colMapLIDs(initialColMapLIDs.size()*dofsPerNode); - for(LO elementIdx = 0; elementIdx < as(initialColMapLIDs.size()); ++elementIdx) { - for(int dof = 0; dof < dofsPerNode; ++dof) { - colMapLIDs[elementIdx*dofsPerNode + dof] = - initialColMapLIDs[elementIdx]*dofsPerNode + dof; - } - } - - *out << "Call column map constructor" << std::endl; - - RCP colMap = MapFactory::Build(prolongatorGraph->getColMap()->lib(), - prolongatorGraph->getGlobalNumCols()*dofsPerNode, - colMapLIDs(), - prolongatorGraph->getIndexBase(), - prolongatorGraph->getComm(), - prolongatorGraph->getColMap()->getNode()); - std::vector strideInfo(1); - strideInfo[0] = dofsPerNode; - RCP stridedDomainMap = StridedMapFactory::Build(domainMap, strideInfo); + strideInfo[0] = A->GetFixedBlockSize(); + RCP stridedDomainMap = + StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo); *out << "Call prolongator constructor" << std::endl; // Create the prolongator matrix and its associated objects - P = rcp(new CrsMatrixWrap(A->getDomainMap(), colMap, 0, Xpetra::StaticProfile)); - // P = rcp(new CrsMatrixWrap(graph)); + RCP dummyList = rcp(new ParameterList()); + P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList)); RCP PCrs = rcp_dynamic_cast(P)->getCrsMatrix(); - - ArrayRCP iaP; - ArrayRCP jaP; - ArrayRCP valP; - - PCrs->allocateAllValues(A->getDomainMap()->getNodeNumElements(), iaP, jaP, valP); - - ArrayView ia = iaP(); - ArrayView ja = jaP(); - ArrayView val = valP(); - ia[0] = 0; - - *out << "Fill prolongator" << std::endl; - - // Actually fill up the matrix, in this case all values are one, pretty easy! - int dofIdx; - ArrayView colIdx; - for(LO rowIdx = 0; rowIdx < static_cast(prolongatorGraph->getNodeNumRows()); ++rowIdx) { - prolongatorGraph->getLocalRowView(rowIdx, colIdx); - for(int dof = 0; dof < dofsPerNode; ++dof) { - dofIdx = rowIdx*dofsPerNode + dof; - ia[dofIdx + 1] = dofIdx + 1; - ja[dofIdx] = colIdx[0]*dofsPerNode + dof; - val[dofIdx] = 1.0; - } - } - - *out << "Set values and call fillComplete on prolongator" << std::endl; - - // call fill complete on the prolongator - PCrs->setAllValues(iaP, jaP, valP); - PCrs->expertStaticFillComplete(domainMap, A->getDomainMap()); + PCrs->setAllToScalar(1.0); + PCrs->fillComplete(); // set StridingInformation of P if (A->IsView("stridedMaps") == true) { @@ -310,7 +239,7 @@ namespace MueLu { template void GeometricInterpolationPFactory:: - BuildLinearP(RCP& A, RCP& prolongatorGraph, + BuildLinearP(RCP& A, RCP& prolongatorGraph, RCP& fineCoordinates, RCP& ghostCoordinates, const int numDimensions, RCP& P) const { @@ -327,11 +256,13 @@ namespace MueLu { *out << "Entering BuildLinearP" << std::endl; // Extract coordinates for interpolation stencil calculations + const LO numFineNodes = fineCoordinates->getLocalLength(); + const LO numGhostNodes = ghostCoordinates->getLocalLength(); Array > fineCoords(3); Array > ghostCoords(3); const real_type realZero = Teuchos::as(0.0); - ArrayRCP fineZero(fineCoordinates->getLocalLength(), realZero); - ArrayRCP ghostZero(ghostCoordinates->getLocalLength(), realZero); + ArrayRCP fineZero(numFineNodes, realZero); + ArrayRCP ghostZero(numGhostNodes, realZero); for(int dim = 0; dim < 3; ++dim) { if(dim < numDimensions) { fineCoords[dim] = fineCoordinates->getData(dim); @@ -347,117 +278,66 @@ namespace MueLu { // Compute 2^numDimensions using bit logic to avoid round-off errors const int numInterpolationPoints = 1 << numDimensions; const int dofsPerNode = A->GetFixedBlockSize(); - ArrayView initialDomainMapLIDs = - prolongatorGraph->getDomainMap()->getNodeElementList(); - Array domainMapLIDs(initialDomainMapLIDs.size()*dofsPerNode); - for(LO elementIdx = 0; elementIdx < as(initialDomainMapLIDs.size()); ++elementIdx) { - for(int dof = 0; dof < dofsPerNode; ++dof) { - domainMapLIDs[elementIdx*dofsPerNode + dof] = - initialDomainMapLIDs[elementIdx]*dofsPerNode + dof; - } - } - RCP domainMap = MapFactory::Build(prolongatorGraph->getRowMap()->lib(), - prolongatorGraph->getGlobalNumCols()*dofsPerNode, - domainMapLIDs(), - prolongatorGraph->getIndexBase(), - prolongatorGraph->getComm(), - prolongatorGraph->getRowMap()->getNode()); - - ArrayView initialColMapLIDs = - prolongatorGraph->getColMap()->getNodeElementList(); - Array colMapLIDs(initialColMapLIDs.size()*dofsPerNode); - for(LO elementIdx = 0; elementIdx < as(initialColMapLIDs.size()); ++elementIdx) { - for(int dof = 0; dof < dofsPerNode; ++dof) { - colMapLIDs[elementIdx*dofsPerNode + dof] = - initialColMapLIDs[elementIdx]*dofsPerNode + dof; - } - } - RCP colMap = MapFactory::Build(prolongatorGraph->getColMap()->lib(), - prolongatorGraph->getGlobalNumCols()*dofsPerNode, - colMapLIDs(), - prolongatorGraph->getIndexBase(), - prolongatorGraph->getComm(), - prolongatorGraph->getColMap()->getNode()); std::vector strideInfo(1); - strideInfo[0] = dofsPerNode; - RCP stridedDomainMap = StridedMapFactory::Build(domainMap, strideInfo); + strideInfo[0] = dofsPerNode; + RCP stridedDomainMap = + StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo); *out << "The maps of P have been computed" << std::endl; - P = rcp(new CrsMatrixWrap(A->getDomainMap(), colMap, 0, Xpetra::StaticProfile)); - // P = rcp(new CrsMatrixWrap(graph)); + RCP dummyList = rcp(new ParameterList()); + P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList)); RCP PCrs = rcp_dynamic_cast(P)->getCrsMatrix(); + PCrs->resumeFill(); // The Epetra matrix is considered filled at this point. - ArrayRCP iaP; - ArrayRCP jaP; - ArrayRCP valP; - - LO numNonZeroP = prolongatorGraph->getNodeNumEntries()*dofsPerNode; - PCrs->allocateAllValues(numNonZeroP, iaP, jaP, valP); - - *out << "dofsPerNode=" << dofsPerNode << std::endl; - *out << "number of non-zeroes in P: " << numNonZeroP << std::endl; - - ArrayView ia = iaP(); - ArrayView ja = jaP(); - ArrayView val = valP(); - ia[0] = 0; - - LO rowDofIdx, colDofIdx; + LO interpolationNodeIdx = 0, rowIdx = 0; ArrayView colIndices; + Array values; Array > coords(numInterpolationPoints + 1); Array stencil(numInterpolationPoints); - for(LO rowIdx = 0; rowIdx < static_cast(prolongatorGraph->getNodeNumRows()); ++rowIdx) { - prolongatorGraph->getLocalRowView(rowIdx, colIndices); - - // rowIdx and colIdx correspond to single ids of nodes on the fine resp. coarse mesh - // rowDofIdx and colDofIdx correspond to single ids of dofs on the fine resp. coarse mesh - if(colIndices.size() == 1) { - // If colIndices.size() == 1, we are handling a coarse node - // we only need to stick a one in the correct location! - for(int dof = 0; dof < dofsPerNode; ++dof) { - rowDofIdx = rowIdx*dofsPerNode + dof; - ia[rowDofIdx + 1] = ia[rowDofIdx] + 1; - - colDofIdx = ia[rowDofIdx]; - ja[colDofIdx] = colIndices[0]*dofsPerNode + dof; - val[colDofIdx] = 1.0; + for(LO nodeIdx = 0; nodeIdx < numFineNodes; ++nodeIdx) { + if(PCrs->getNumEntriesInLocalRow(nodeIdx*dofsPerNode) == 1) { + values.resize(1); + values[0] = 1.0; + for(LO dof = 0; dof < dofsPerNode; ++dof) { + rowIdx = nodeIdx*dofsPerNode + dof; + prolongatorGraph->getLocalRowView(rowIdx, colIndices); + PCrs->replaceLocalValues(rowIdx, colIndices, values()); } } else { - // If colIndices.size() > 1, we need to compute the linear interpolation coefficients - for(int dof = 0; dof < dofsPerNode; ++dof) { - rowDofIdx = rowIdx*dofsPerNode + dof; - ia[rowDofIdx + 1] = ia[rowDofIdx] + colIndices.size(); - - // Extract the coordinates associated with the current node - // and the neighboring coarse nodes - coords[0].resize(3); + // Extract the coordinates associated with the current node + // and the neighboring coarse nodes + coords[0].resize(3); + for(int dim = 0; dim < 3; ++dim) { + coords[0][dim] = fineCoords[dim][nodeIdx]; + } + prolongatorGraph->getLocalRowView(nodeIdx*dofsPerNode, colIndices); + for(int interpolationIdx=0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) { + coords[interpolationIdx + 1].resize(3); + interpolationNodeIdx = colIndices[interpolationIdx] / dofsPerNode; for(int dim = 0; dim < 3; ++dim) { - coords[0][dim] = fineCoords[dim][rowIdx]; - } - for(int interpolationIdx = 0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) { - coords[interpolationIdx + 1].resize(3); - for(int dim = 0; dim < 3; ++dim) { - coords[interpolationIdx + 1][dim] = ghostCoords[dim][colIndices[interpolationIdx]]; - } + coords[interpolationIdx + 1][dim] = ghostCoords[dim][interpolationNodeIdx]; } - ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil); + } + ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil); + values.resize(numInterpolationPoints); + for(LO valueIdx = 0; valueIdx < numInterpolationPoints; ++valueIdx) { + values[valueIdx] = stencil[valueIdx]; + } - for(int colIdx = 0; colIdx < colIndices.size(); ++colIdx) { - colDofIdx = ia[rowDofIdx] + colIdx; - ja[colDofIdx] = colIndices[colIdx]*dofsPerNode + dof; - val[colDofIdx] = stencil[colIdx]; - } + // Set values in all the rows corresponding to nodeIdx + for(LO dof = 0; dof < dofsPerNode; ++dof) { + rowIdx = nodeIdx*dofsPerNode + dof; + prolongatorGraph->getLocalRowView(rowIdx, colIndices); + PCrs->replaceLocalValues(rowIdx, colIndices, values()); } } } *out << "The calculation of the interpolation stencils has completed." << std::endl; - Xpetra::CrsMatrixUtils::sortCrsEntries(ia, ja, val, A->getDomainMap()->lib()); - PCrs->setAllValues(iaP, jaP, valP); - PCrs->expertStaticFillComplete(domainMap, A->getDomainMap()); + PCrs->fillComplete(); *out << "All values in P have been set and expertStaticFillComplete has been performed." << std::endl; diff --git a/packages/muelu/test/structured/structured_2dof.xml b/packages/muelu/test/structured/structured_2dof.xml index d33f330ceeb5..7d17dda782c0 100644 --- a/packages/muelu/test/structured/structured_2dof.xml +++ b/packages/muelu/test/structured/structured_2dof.xml @@ -18,9 +18,10 @@ - + + @@ -32,7 +33,7 @@ - + diff --git a/packages/muelu/test/structured/structured_3dof.xml b/packages/muelu/test/structured/structured_3dof.xml index e2646e5ab7bf..c1db49e01fa9 100644 --- a/packages/muelu/test/structured/structured_3dof.xml +++ b/packages/muelu/test/structured/structured_3dof.xml @@ -21,6 +21,7 @@ + diff --git a/packages/muelu/test/unit_tests/StructuredAggregationFactory.cpp b/packages/muelu/test/unit_tests/StructuredAggregationFactory.cpp index 337f28895e71..f3bfeb91a1d6 100644 --- a/packages/muelu/test/unit_tests/StructuredAggregationFactory.cpp +++ b/packages/muelu/test/unit_tests/StructuredAggregationFactory.cpp @@ -128,6 +128,7 @@ namespace MueLuTests { dropFact->SetFactory("UnAmalgamationInfo", amalgFact); RCP StructuredAggFact = rcp(new StructuredAggregationFactory()); StructuredAggFact->SetFactory("Graph", dropFact); + StructuredAggFact->SetFactory("DofsPerNode", dropFact); StructuredAggFact->SetParameter("aggregation: mesh layout", Teuchos::ParameterEntry(meshLayout)); StructuredAggFact->SetParameter("aggregation: number of spatial dimensions", @@ -232,6 +233,7 @@ namespace MueLuTests { dropFact->SetFactory("UnAmalgamationInfo", amalgFact); RCP StructuredAggFact = rcp(new StructuredAggregationFactory()); StructuredAggFact->SetFactory("Graph", dropFact); + StructuredAggFact->SetFactory("DofsPerNode", dropFact); StructuredAggFact->SetParameter("aggregation: mesh layout", Teuchos::ParameterEntry(meshLayout)); StructuredAggFact->SetParameter("aggregation: number of spatial dimensions", @@ -336,6 +338,7 @@ namespace MueLuTests { dropFact->SetFactory("UnAmalgamationInfo", amalgFact); RCP StructuredAggFact = rcp(new StructuredAggregationFactory()); StructuredAggFact->SetFactory("Graph", dropFact); + StructuredAggFact->SetFactory("DofsPerNode", dropFact); StructuredAggFact->SetParameter("aggregation: mesh layout", Teuchos::ParameterEntry(meshLayout)); StructuredAggFact->SetParameter("aggregation: number of spatial dimensions", @@ -442,6 +445,7 @@ namespace MueLuTests { dropFact->SetFactory("UnAmalgamationInfo", amalgFact); RCP StructuredAggFact = rcp(new StructuredAggregationFactory()); StructuredAggFact->SetFactory("Graph", dropFact); + StructuredAggFact->SetFactory("DofsPerNode", dropFact); StructuredAggFact->SetParameter("aggregation: mesh layout", Teuchos::ParameterEntry(meshLayout)); StructuredAggFact->SetParameter("aggregation: number of spatial dimensions", @@ -548,6 +552,7 @@ namespace MueLuTests { dropFact->SetFactory("UnAmalgamationInfo", amalgFact); RCP StructuredAggFact = rcp(new StructuredAggregationFactory()); StructuredAggFact->SetFactory("Graph", dropFact); + StructuredAggFact->SetFactory("DofsPerNode", dropFact); StructuredAggFact->SetParameter("aggregation: mesh layout", Teuchos::ParameterEntry(meshLayout)); StructuredAggFact->SetParameter("aggregation: number of spatial dimensions", @@ -655,6 +660,7 @@ namespace MueLuTests { dropFact->SetFactory("UnAmalgamationInfo", amalgFact); RCP StructuredAggFact = rcp(new StructuredAggregationFactory()); StructuredAggFact->SetFactory("Graph", dropFact); + StructuredAggFact->SetFactory("DofsPerNode", dropFact); StructuredAggFact->SetParameter("aggregation: mesh layout", Teuchos::ParameterEntry(meshLayout)); StructuredAggFact->SetParameter("aggregation: number of spatial dimensions", @@ -767,6 +773,7 @@ namespace MueLuTests { dropFact->SetFactory("UnAmalgamationInfo", amalgFact); RCP StructuredAggFact = rcp(new StructuredAggregationFactory()); StructuredAggFact->SetFactory("Graph", dropFact); + StructuredAggFact->SetFactory("DofsPerNode", dropFact); StructuredAggFact->SetParameter("aggregation: mesh layout", Teuchos::ParameterEntry(meshLayout)); StructuredAggFact->SetParameter("aggregation: coupling", @@ -881,6 +888,7 @@ namespace MueLuTests { dropFact->SetFactory("UnAmalgamationInfo", amalgFact); RCP StructuredAggFact = rcp(new StructuredAggregationFactory()); StructuredAggFact->SetFactory("Graph", dropFact); + StructuredAggFact->SetFactory("DofsPerNode", dropFact); StructuredAggFact->SetParameter("aggregation: mesh layout", Teuchos::ParameterEntry(meshLayout)); StructuredAggFact->SetParameter("aggregation: coupling", @@ -996,6 +1004,7 @@ namespace MueLuTests { dropFact->SetFactory("UnAmalgamationInfo", amalgFact); RCP StructuredAggFact = rcp(new StructuredAggregationFactory()); StructuredAggFact->SetFactory("Graph", dropFact); + StructuredAggFact->SetFactory("DofsPerNode", dropFact); StructuredAggFact->SetParameter("aggregation: mesh layout", Teuchos::ParameterEntry(meshLayout)); StructuredAggFact->SetParameter("aggregation: coupling", @@ -1126,6 +1135,7 @@ namespace MueLuTests { // Set interfactory dependencies CDropfact->SetFactory("UnAmalgamationInfo", AmalgFact); Aggfact->SetFactory("Graph", CDropfact); + Aggfact->SetFactory("DofsPerNode", CDropfact); coarseMapFact->SetFactory("Aggregates", Aggfact); Pfact->SetFactory("Aggregates", Aggfact); Pfact->SetFactory("CoarseMap", coarseMapFact);