Skip to content

GitLab

Unverified Commit 4f15e6fb authored by William's avatar William Committed by GitHub
Browse files

Merge pull request #4215 from trilinos/master_merge_20190118_020615 

Trilinos Master Merge PR Generator: Auto PR created to promote from master_merge_20190118_020615 branch to master
parents 59b802f2 28ff0dfc
Expand all Hide whitespace changes
Inline Side-by-side
packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_decl.hpp
View file @ 4f15e6fb
......@@ -106,7 +106,7 @@ namespace MueLu {
/*! @brief Local aggregation. */
void BuildGraph(const GraphBase& graph, RCP<IndexManager>& geoData,
void BuildGraph(const GraphBase& graph, RCP<IndexManager>& geoData, const LO dofsPerNode,
RCP<CrsGraph>& myGraph, RCP<const Map>& coarseCoordinatesFineMap,
RCP<const Map>& coarseCoordinatesMap) const;
//@}
......@@ -116,12 +116,14 @@ namespace MueLu {
private:
void ComputeGraphDataConstant(const GraphBase& graph, RCP<IndexManager>& geoData,
const int numInterpolationPoints, ArrayRCP<size_t>& nnzOnRow,
Array<size_t>& rowPtr, Array<LO>& colIndex) const;
const LO dofsPerNode, const int numInterpolationPoints,
ArrayRCP<size_t>& nnzOnRow, Array<size_t>& rowPtr,
Array<LO>& colIndex) const;
void ComputeGraphDataLinear(const GraphBase& graph, RCP<IndexManager>& geoData,
const int numInterpolationPoints, ArrayRCP<size_t>& nnzOnRow,
Array<size_t>& rowPtr, Array<LO>& colIndex) const;
const LO dofsPerNode, const int numInterpolationPoints,
ArrayRCP<size_t>& nnzOnRow, Array<size_t>& rowPtr,
Array<LO>& colIndex) const;
};
......
packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_def.hpp
View file @ 4f15e6fb
......@@ -129,8 +129,9 @@ namespace MueLu {
template <class LocalOrdinal, class GlobalOrdinal, class Node>
void AggregationStructuredAlgorithm<LocalOrdinal, GlobalOrdinal, Node>::
BuildGraph(const GraphBase& graph, RCP<IndexManager>& geoData, RCP<CrsGraph>& myGraph,
RCP<const Map>& coarseCoordinatesFineMap, RCP<const Map>& coarseCoordinatesMap) const {
BuildGraph(const GraphBase& graph, RCP<IndexManager>& geoData, const LO dofsPerNode,
RCP<CrsGraph>& myGraph, RCP<const Map>& coarseCoordinatesFineMap,
RCP<const Map>& coarseCoordinatesMap) const {
Monitor m(*this, "BuildGraphP");
RCP<Teuchos::FancyOStream> out;
......@@ -153,20 +154,23 @@ namespace MueLu {
}
*out << "numInterpolationPoints=" << numInterpolationPoints << std::endl;
Array<LO> colIndex( geoData->getNumLocalCoarseNodes() + numInterpolationPoints*
(geoData->getNumLocalFineNodes() - geoData->getNumLocalCoarseNodes()) );
Array<size_t> rowPtr(geoData->getNumLocalFineNodes()+1);
Array<LO> colIndex((geoData->getNumLocalCoarseNodes() + numInterpolationPoints*
(geoData->getNumLocalFineNodes() - geoData->getNumLocalCoarseNodes()))*dofsPerNode);
Array<size_t> rowPtr(geoData->getNumLocalFineNodes()*dofsPerNode + 1);
rowPtr[0] = 0;
ArrayRCP<size_t> nnzOnRow(geoData->getNumLocalFineNodes());
ArrayRCP<size_t> 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<Map> rowMap = MapFactory::Build(graph.GetDomainMap(), dofsPerNode);
RCP<Map> 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<GO>::invalid()
// and it will assign GIDs to the local coarse nodes.
colMap = MapFactory::Build(graph.GetDomainMap()->lib(),
geoData->getNumGlobalCoarseNodes(),
geoData->getNumLocalCoarseNodes(),
Teuchos::OrdinalTraits<GO>::invalid(),
geoData->getNumLocalCoarseNodes()*dofsPerNode,
graph.GetDomainMap()->getIndexBase(),
graph.GetDomainMap()->getComm(),
graph.GetDomainMap()->getNode());
......@@ -229,13 +232,13 @@ namespace MueLu {
Array<GO> coarseNodeFineGIDs(geoData->getNumLocalCoarseNodes());
geoData->getCoarseNodesData(graph.GetDomainMap(), coarseNodeCoarseGIDs, coarseNodeFineGIDs);
coarseCoordinatesMap = MapFactory::Build(graph.GetDomainMap()->lib(),
geoData->getNumGlobalCoarseNodes(),
Teuchos::OrdinalTraits<GO>::invalid(),
geoData->getNumLocalCoarseNodes(),
graph.GetDomainMap()->getIndexBase(),
graph.GetDomainMap()->getComm(),
graph.GetDomainMap()->getNode());
coarseCoordinatesFineMap = MapFactory::Build(graph.GetDomainMap()->lib(),
geoData->getNumGlobalCoarseNodes(),
Teuchos::OrdinalTraits<GO>::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 <class LocalOrdinal, class GlobalOrdinal, class Node>
void AggregationStructuredAlgorithm<LocalOrdinal, GlobalOrdinal, Node>::
ComputeGraphDataConstant(const GraphBase& graph, RCP<IndexManager>& geoData,
const int numInterpolationPoints, ArrayRCP<size_t>& nnzOnRow,
Array<size_t>& rowPtr, Array<LO>& colIndex) const {
const LO dofsPerNode, const int numInterpolationPoints,
ArrayRCP<size_t>& nnzOnRow, Array<size_t>& rowPtr,
Array<LO>& colIndex) const {
RCP<Teuchos::FancyOStream> out;
if(const char* dbg = std::getenv("MUELU_STRUCTUREDALGORITHM_DEBUG")) {
......@@ -283,9 +291,6 @@ namespace MueLu {
LO ghostedCoarseNodeCoarseLID, rem, rate;
Array<LO> 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<size_t>(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 <class LocalOrdinal, class GlobalOrdinal, class Node>
void AggregationStructuredAlgorithm<LocalOrdinal, GlobalOrdinal, Node>::
ComputeGraphDataLinear(const GraphBase& graph, RCP<IndexManager>& geoData,
const int numInterpolationPoints, ArrayRCP<size_t>& nnzOnRow,
Array<size_t>& rowPtr, Array<LO>& colIndex) const {
const LO dofsPerNode, const int numInterpolationPoints,
ArrayRCP<size_t>& nnzOnRow, Array<size_t>& rowPtr,
Array<LO>& colIndex) const {
RCP<Teuchos::FancyOStream> out;
if(const char* dbg = std::getenv("MUELU_STRUCTUREDALGORITHM_DEBUG")) {
......@@ -332,6 +344,8 @@ namespace MueLu {
Array<LO> coarseIdx(3,0);
Array<LO> 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<size_t>(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<size_t>( numInterpolationPoints );
rowPtr[nodeIdx + 1] = rowPtr[nodeIdx] + Teuchos::as<LO>( 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<size_t>( numInterpolationPoints );
rowPtr[rowIdx + 1] = rowPtr[rowIdx] + Teuchos::as<LO>(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()
......
packages/muelu/src/Graph/StructuredAggregation/MueLu_StructuredAggregationFactory_def.hpp
View file @ 4f15e6fb
......@@ -97,7 +97,6 @@ namespace MueLu {
validParamList->set<RCP<const FactoryBase> >("aggregation: mesh data", Teuchos::null,
"Mesh ordering associated data");
validParamList->set<RCP<const FactoryBase> >("Graph", Teuchos::null,
"Graph of the matrix after amalgamation but without dropping.");
validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
......@@ -106,6 +105,8 @@ namespace MueLu {
"Global number of nodes per spatial dimension provided by CoordinatesTransferFactory.");
validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null,
"Local number of nodes per spatial dimension provided by CoordinatesTransferFactory.");
validParamList->set<RCP<const FactoryBase> >("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<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
DeclareInput(Level& currentLevel) const {
Input(currentLevel, "Graph");
Input(currentLevel, "DofsPerNode");
ParameterList pL = GetParameterList();
std::string coupling = pL.get<std::string>("aggregation: coupling");
......@@ -175,10 +177,11 @@ namespace MueLu {
// General problem informations are gathered from data stored in the problem matix.
RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel, "Graph");
RCP<const Map> fineMap = graph->GetDomainMap();
const int myRank = fineMap->getComm()->getRank();
const int numRanks = fineMap->getComm()->getSize();
const GO minGlobalIndex = fineMap->getMinGlobalIndex();
RCP<const Map> fineMap = graph->GetDomainMap();
const int myRank = fineMap->getComm()->getRank();
const int numRanks = fineMap->getComm()->getSize();
const GO minGlobalIndex = fineMap->getMinGlobalIndex();
const LO dofsPerNode = Get<LO>(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<CrsGraph> myGraph;
myStructuredAlgorithm->BuildGraph(*graph, geoData, myGraph, coarseCoordinatesFineMap,
coarseCoordinatesMap);
myStructuredAlgorithm->BuildGraph(*graph, geoData, dofsPerNode, myGraph,
coarseCoordinatesFineMap, coarseCoordinatesMap);
Set(currentLevel, "prolongatorGraph", myGraph);
}
......
packages/muelu/src/Transfers/GeneralGeometric/MueLu_GeometricInterpolationPFactory_decl.hpp
View file @ 4f15e6fb
......@@ -98,8 +98,8 @@ namespace MueLu{
//@}
private:
void BuildConstantP(RCP<Matrix>& P, RCP<CrsGraph>& prolongatorGraph, RCP<Matrix>& A) const;
void BuildLinearP(RCP<Matrix>& A, RCP<CrsGraph>& prolongatorGraph,
void BuildConstantP(RCP<Matrix>& P, RCP<const CrsGraph>& prolongatorGraph, RCP<Matrix>& A) const;
void BuildLinearP(RCP<Matrix>& A, RCP<const CrsGraph>& prolongatorGraph,
RCP<realvaluedmultivector_type>& fineCoordinates,
RCP<realvaluedmultivector_type>& ghostCoordinates,
const int numDimensions, RCP<Matrix>& P) const;
......
packages/muelu/src/Transfers/GeneralGeometric/MueLu_GeometricInterpolationPFactory_def.hpp
View file @ 4f15e6fb
......@@ -138,7 +138,7 @@ namespace MueLu {
// Declared main input/outputs to be retrieved and placed on the fine resp. coarse level
RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
RCP<CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel, "prolongatorGraph");
RCP<const CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel, "prolongatorGraph");
RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
RCP<Matrix> P;
const int interpolationOrder = pL.get<int>("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<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel, "Nullspace");
......@@ -201,7 +201,7 @@ namespace MueLu {
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
BuildConstantP(RCP<Matrix>& P, RCP<CrsGraph>& prolongatorGraph, RCP<Matrix>& A) const {
BuildConstantP(RCP<Matrix>& P, RCP<const CrsGraph>& prolongatorGraph, RCP<Matrix>& A) const {
// Set debug outputs based on environment variable
RCP<Teuchos::FancyOStream> 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<const GO> initialDomainMapLIDs =
prolongatorGraph->getDomainMap()->getNodeElementList();
Array<GO> domainMapLIDs(initialDomainMapLIDs.size()*dofsPerNode);
for(LO elementIdx = 0; elementIdx < as<LO>(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<Map> domainMap = MapFactory::Build(prolongatorGraph->getRowMap()->lib(),
prolongatorGraph->getGlobalNumCols()*dofsPerNode,
domainMapLIDs(),
prolongatorGraph->getIndexBase(),
prolongatorGraph->getComm(),
prolongatorGraph->getRowMap()->getNode());
ArrayView<const GO> initialColMapLIDs =
prolongatorGraph->getColMap()->getNodeElementList();
Array<GO> colMapLIDs(initialColMapLIDs.size()*dofsPerNode);
for(LO elementIdx = 0; elementIdx < as<LO>(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<Map> colMap = MapFactory::Build(prolongatorGraph->getColMap()->lib(),
prolongatorGraph->getGlobalNumCols()*dofsPerNode,
colMapLIDs(),
prolongatorGraph->getIndexBase(),
prolongatorGraph->getComm(),
prolongatorGraph->getColMap()->getNode());
std::vector<size_t> strideInfo(1);
strideInfo[0] = dofsPerNode;
RCP<const StridedMap> stridedDomainMap = StridedMapFactory::Build(domainMap, strideInfo);
strideInfo[0] = A->GetFixedBlockSize();
RCP<const StridedMap> 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<ParameterList> dummyList = rcp(new ParameterList());
P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
ArrayRCP<size_t> iaP;
ArrayRCP<LO> jaP;
ArrayRCP<SC> valP;
PCrs->allocateAllValues(A->getDomainMap()->getNodeNumElements(), iaP, jaP, valP);
ArrayView<size_t> ia = iaP();
ArrayView<LO> ja = jaP();
ArrayView<SC> 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<const LO> colIdx;
for(LO rowIdx = 0; rowIdx < static_cast<LO>(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 <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
BuildLinearP(RCP<Matrix>& A, RCP<CrsGraph>& prolongatorGraph,
BuildLinearP(RCP<Matrix>& A, RCP<const CrsGraph>& prolongatorGraph,
RCP<realvaluedmultivector_type>& fineCoordinates,
RCP<realvaluedmultivector_type>& ghostCoordinates,
const int numDimensions, RCP<Matrix>& 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<ArrayRCP<const real_type> > fineCoords(3);
Array<ArrayRCP<const real_type> > ghostCoords(3);
const real_type realZero = Teuchos::as<real_type>(0.0);
ArrayRCP<real_type> fineZero(fineCoordinates->getLocalLength(), realZero);
ArrayRCP<real_type> ghostZero(ghostCoordinates->getLocalLength(), realZero);
ArrayRCP<real_type> fineZero(numFineNodes, realZero);
ArrayRCP<real_type> 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<const GO> initialDomainMapLIDs =
prolongatorGraph->getDomainMap()->getNodeElementList();
Array<GO> domainMapLIDs(initialDomainMapLIDs.size()*dofsPerNode);
for(LO elementIdx = 0; elementIdx < as<LO>(initialDomainMapLIDs.size()); ++elementIdx) {
for(int dof = 0; dof < dofsPerNode; ++dof) {
domainMapLIDs[elementIdx*dofsPerNode + dof] =
initialDomainMapLIDs[elementIdx]*dofsPerNode + dof;
}
}
RCP<Map> domainMap = MapFactory::Build(prolongatorGraph->getRowMap()->lib(),
prolongatorGraph->getGlobalNumCols()*dofsPerNode,
domainMapLIDs(),
prolongatorGraph->getIndexBase(),
prolongatorGraph->getComm(),
prolongatorGraph->getRowMap()->getNode());
ArrayView<const GO> initialColMapLIDs =
prolongatorGraph->getColMap()->getNodeElementList();
Array<GO> colMapLIDs(initialColMapLIDs.size()*dofsPerNode);
for(LO elementIdx = 0; elementIdx < as<LO>(initialColMapLIDs.size()); ++elementIdx) {
for(int dof = 0; dof < dofsPerNode; ++dof) {
colMapLIDs[elementIdx*dofsPerNode + dof] =
initialColMapLIDs[elementIdx]*dofsPerNode + dof;
}
}
RCP<Map> colMap = MapFactory::Build(prolongatorGraph->getColMap()->lib(),
prolongatorGraph->getGlobalNumCols()*dofsPerNode,
colMapLIDs(),
prolongatorGraph->getIndexBase(),
prolongatorGraph->getComm(),
prolongatorGraph->getColMap()->getNode());
std::vector<size_t> strideInfo(1);
strideInfo[0] = dofsPerNode;
RCP<const StridedMap> stridedDomainMap = StridedMapFactory::Build(domainMap, strideInfo);
strideInfo[0] = dofsPerNode;
RCP<const StridedMap> 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<ParameterList> dummyList = rcp(new ParameterList());
P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
PCrs->resumeFill(); // The Epetra matrix is considered filled at this point.
ArrayRCP<size_t> iaP;
ArrayRCP<LO> jaP;
ArrayRCP<SC> 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<size_t> ia = iaP();
ArrayView<LO> ja = jaP();
ArrayView<SC> val = valP();
ia[0] = 0;
LO rowDofIdx, colDofIdx;
LO interpolationNodeIdx = 0, rowIdx = 0;
ArrayView<const LO> colIndices;
Array<SC> values;
Array<Array<real_type> > coords(numInterpolationPoints + 1);
Array<real_type> stencil(numInterpolationPoints);
for(LO rowIdx = 0; rowIdx < static_cast<LO>(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