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);