Skip to content

GitLab

Unverified Commit 852d7e71 authored by Matthias Mayr's avatar Matthias Mayr
Browse files

Split setup & solve of regionMG direct solver 

- add new routine MakeCompositeDirectSolver() to create an Amesos2
  KLU2 solver based on a composite operator and pre-compute its
  factorization.
- pass the direct solver through the V-cycle
- just apple solver on the coarsest level

Related to #5072.
parent 67938857
Hide whitespace changes
Inline Side-by-side
packages/muelu/research/mmayr/composite_to_regions/src/SetupRegionHierarchy_def.hpp
View file @ 852d7e71
......@@ -1103,6 +1103,45 @@ void MakeCoarseCompositeOperator(const int maxRegPerProc, const int numLevels,
} // MakeCoarseCompositeOperator
/* Create a direct solver for a composite operator
*
* Create the solver object and compute symbolic and numeric factorization.
* Finally, the solver object will be ready to be applied during the V-cycle call.
*
* \note For now, we're limited to Tpetra/Amesos2. From Amesos2, we use KLU as direct solver.
*/
template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<Amesos2::Solver<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >
MakeCompositeDirectSolver(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& compOp)
{
using Tpetra_CrsMatrix = Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
using Tpetra_MultiVector = Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
using Utilities = MueLu::Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
// RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
// fos->setOutputToRootOnly(0);
// *fos << "Attempting to setup Amesos2 solver the composite coarse grid problem" << std::endl;
// convert matrix to Tpetra
RCP<Tpetra_CrsMatrix> tMat = Utilities::Op2NonConstTpetraCrs(compOp);
// Amesos2-specific key phrase that denote smoother type
std::string amesos2SolverName = "KLU2";
RCP<Amesos2::Solver<Tpetra_CrsMatrix, Tpetra_MultiVector> > coarseSolver;
TEUCHOS_ASSERT(Amesos2::query(amesos2SolverName));
coarseSolver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(amesos2SolverName, tMat);
Teuchos::ParameterList amesos2_params("Amesos2");
amesos2_params.sublist(amesos2SolverName).set("IsContiguous", false, "Are GIDs Contiguous");
coarseSolver->setParameters(Teuchos::rcpFromRef(amesos2_params));
coarseSolver->symbolicFactorization();
coarseSolver->numericFactorization();
return coarseSolver;
} // MakeCorseCompositeDirectSolver
// Make interface scaling factors recursively
template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MakeInterfaceScalingFactors(const int maxRegPerProc,
......@@ -1168,8 +1207,10 @@ void createRegionHierarchy(const int maxRegPerProc,
Array<std::vector<RCP<Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node> > > >& regRowImporters,
Array<std::vector<RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > > >& regInterfaceScalings,
RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& coarseCompOp,
const int maxRegPerGID = 0,
ArrayView<LocalOrdinal> compositeToRegionLIDs = {})
const int maxRegPerGID,
ArrayView<LocalOrdinal> compositeToRegionLIDs,
RCP<Amesos2::Solver<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >& coarseSolver
)
{
#include "Xpetra_UseShortNames.hpp"
......@@ -1314,6 +1355,10 @@ void createRegionHierarchy(const int maxRegPerProc,
regMatrices,
coarseCompOp);
std::cout << mapComp->getComm()->getRank() << " | MakeCoarseCompositeDirectSolver ..." << std::endl;
coarseSolver = MakeCompositeDirectSolver(coarseCompOp);
std::cout << mapComp->getComm()->getRank() << " | MakeInterfaceScalingFactors ..." << std::endl;
MakeInterfaceScalingFactors(maxRegPerProc,
......@@ -1326,6 +1371,48 @@ void createRegionHierarchy(const int maxRegPerProc,
} // createRegionHierarchy
// Wrapper to be used from the Matlab-based driver
template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void createRegionHierarchy(const int maxRegPerProc,
const int numDimensions,
const Array<Array<int> > lNodesPerDim,
const Array<std::string> aggregationRegionType,
const std::string xmlFileName,
Array<RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >& nullspace,
Array<RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > >& coordinates,
std::vector<RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >& regionGrpMats,
const RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > mapComp,
const std::vector<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > > rowMapPerGrp,
const std::vector<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > > colMapPerGrp,
const std::vector<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > > revisedRowMapPerGrp,
const std::vector<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > > revisedColMapPerGrp,
const std::vector<RCP<Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node> > > rowImportPerGrp,
Array<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > >& compRowMaps,
Array<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > >& compColMaps,
Array<std::vector<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > > >& regRowMaps,
Array<std::vector<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > > >& regColMaps,
Array<std::vector<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > > >& quasiRegRowMaps,
Array<std::vector<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > > >& quasiRegColMaps,
Array<std::vector<RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > > >& regMatrices,
Array<std::vector<RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > > >& regProlong,
Array<std::vector<RCP<Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node> > > >& regRowImporters,
Array<std::vector<RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > > >& regInterfaceScalings,
RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& coarseCompOp
)
{
// Define dummy values
const int maxRegPerGID = 0;
ArrayView<LocalOrdinal> compositeToRegionLIDs = {};
RCP<Amesos2::Solver<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>,
Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > > coarseSolver = Teuchos::null;
// Call the actual routine
createRegionHierarchy(maxRegPerProc, numDimensions, lNodesPerDim,
aggregationRegionType, xmlFileName, nullspace, coordinates,
regionGrpMats, mapComp, rowMapPerGrp, colMapPerGrp, revisedRowMapPerGrp, revisedColMapPerGrp, rowImportPerGrp,
compRowMaps, compColMaps, regRowMaps, regColMaps, quasiRegRowMaps, quasiRegColMaps, regMatrices, regProlong,
regRowImporters, regInterfaceScalings, coarseCompOp, maxRegPerGID, compositeToRegionLIDs, coarseSolver);
}
......@@ -1505,7 +1592,8 @@ void vCycle(const int l, ///< ID of current level
Array<std::vector<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > > > regRowMaps, ///< regional row maps
Array<std::vector<RCP<Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node> > > > regRowImporters, ///< regional row importers
Array<std::vector<RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > > > regInterfaceScalings, ///< regional interface scaling factors
RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > coarseCompMat ///< Coarsest level composite operator
RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > coarseCompMat, ///< Coarsest level composite operator
RCP<Amesos2::Solver<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > > coarseSolver = Teuchos::null ///< Coarsest level composite direct solver
)
{
#include "Xpetra_UseShortNames.hpp"
......@@ -1545,7 +1633,7 @@ void vCycle(const int l, ///< ID of current level
// Call V-cycle recursively
vCycle(l+1, numLevels, maxFineIter, maxCoarseIter, omega, maxRegPerProc,
coarseRegX, coarseRegB, regMatrices, regProlong, compRowMaps,
quasiRegRowMaps, regRowMaps, regRowImporters, regInterfaceScalings, coarseCompMat);
quasiRegRowMaps, regRowMaps, regRowImporters, regInterfaceScalings, coarseCompMat, coarseSolver);
// Transfer coarse level correction to fine level
std::vector<RCP<Vector> > regCorrection(maxRegPerProc);
......@@ -1577,153 +1665,56 @@ void vCycle(const int l, ///< ID of current level
#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_AMESOS2)
if(coarseCompMat->getRowMap()->lib() == Xpetra::UseTpetra) {
using Utilities = MueLu::Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
TEUCHOS_TEST_FOR_EXCEPT_MSG(coarseCompMat->getRowMap()->lib()!=Xpetra::UseTpetra,
"Coarse solver requires Tpetra/Amesos2 stack.");
TEUCHOS_ASSERT(!coarseSolver.is_null());
// First get the Xpetra vectors from region to composite format
// (the coarseCompMat should already exist)
RCP<Vector> compX = VectorFactory::Build(coarseCompMat->getRowMap(), true);
RCP<Vector> compRhs = VectorFactory::Build(coarseCompMat->getRowMap(), true);
{
for (int j = 0; j < maxRegPerProc; j++) {
RCP<Vector> inverseInterfaceScaling = VectorFactory::Build(regInterfaceScalings[l][j]->getMap());
inverseInterfaceScaling->reciprocal(*regInterfaceScalings[l][j]);
fineRegB[j]->elementWiseMultiply(SC_ONE, *fineRegB[j], *inverseInterfaceScaling, SC_ZERO);
}
regionalToComposite(fineRegB, compRhs, maxRegPerProc, quasiRegRowMaps[l],
regRowImporters[l], Xpetra::ADD);
}
// From here on we switch to Tpetra for simplicity
// we could also implement a similar Epetra branch
using Tpetra_CrsMatrix = Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
using Tpetra_MultiVector = Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
*fos << "Attempting to use Amesos2 to solve the coarse grid problem" << std::endl;
RCP<Tpetra_CrsMatrix> tMat = Utilities::Op2NonConstTpetraCrs(coarseCompMat);
RCP<Tpetra_MultiVector> tX = Utilities::MV2NonConstTpetraMV2(*compX);
RCP<Tpetra_MultiVector> tB = Utilities::MV2NonConstTpetraMV2(*compRhs);
// Now we can modify the row map of tMat, tX and tB to make it Amesos friendly
// With all the pieces in place we can call Amesos to solve the coarse grid problem
//! amesos2-specific key phrase that denote smoother type
std::string amesos2SolverName = "KLU2";
RCP<Amesos2::Solver<Tpetra_CrsMatrix, Tpetra_MultiVector> > coarseSolver;
TEUCHOS_ASSERT(Amesos2::query(amesos2SolverName));
coarseSolver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(amesos2SolverName, tMat);
Teuchos::ParameterList amesos2_params("Amesos2");
amesos2_params.sublist(amesos2SolverName).set("IsContiguous", false, "Are GIDs Contiguous");
coarseSolver->setParameters( Teuchos::rcpFromRef(amesos2_params) );
using Utilities = MueLu::Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
// First get the Xpetra vectors from region to composite format
// (the coarseCompMat should already exist)
RCP<Vector> compX = VectorFactory::Build(coarseCompMat->getRowMap(), true);
RCP<Vector> compRhs = VectorFactory::Build(coarseCompMat->getRowMap(), true);
{
for (int j = 0; j < maxRegPerProc; j++) {
RCP<Vector> inverseInterfaceScaling = VectorFactory::Build(regInterfaceScalings[l][j]->getMap());
inverseInterfaceScaling->reciprocal(*regInterfaceScalings[l][j]);
fineRegB[j]->elementWiseMultiply(SC_ONE, *fineRegB[j], *inverseInterfaceScaling, SC_ZERO);
}
coarseSolver->setX(tX);
coarseSolver->setB(tB);
regionalToComposite(fineRegB, compRhs, maxRegPerProc, quasiRegRowMaps[l],
regRowImporters[l], Xpetra::ADD);
}
coarseSolver->solve();
// From here on we switch to Tpetra for simplicity
// we could also implement a similar Epetra branch
using Tpetra_MultiVector = Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
coarseSolver->setX(Teuchos::null);
coarseSolver->setB(Teuchos::null);
// *fos << "Attempting to use Amesos2 to solve the coarse grid problem" << std::endl;
RCP<Tpetra_MultiVector> tX = Utilities::MV2NonConstTpetraMV2(*compX);
RCP<const Tpetra_MultiVector> tB = Utilities::MV2TpetraMV(compRhs);
std::vector<RCP<Vector> > quasiRegX(maxRegPerProc);
compositeToRegional(compX, quasiRegX, fineRegX,
maxRegPerProc,
quasiRegRowMaps[l],
regRowMaps[l],
regRowImporters[l]);
}
/* Solve!
*
* We don't have to change the map of tX and tB since we have configured the Amesos2 solver
* during its construction to work with non-continous maps.
*/
coarseSolver->solve(tX.ptr(), tB.ptr());
// Transform back to region format
std::vector<RCP<Vector> > quasiRegX(maxRegPerProc);
compositeToRegional(compX, quasiRegX, fineRegX,
maxRegPerProc,
quasiRegRowMaps[l],
regRowMaps[l],
regRowImporters[l]);
#else
*fos << "+++++++++++++++++++++++++++ WARNING +++++++++++++++++++++++++\n"
<< "+ Coarse level solver has not been migrated to Xpetra, yet. +\n"
<< "+ We skip it for now. +\n"
<< "+ Coarse level solver requires Tpetra and Amesos2. +\n"
<< "+ Skipping the coarse level solve. +\n"
<< "+++++++++++++++++++++++++++ WARNING +++++++++++++++++++++++++"
<< std::endl;
<< std::endl;
#endif
// TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Coarse level solver not migrated to Xpetra, yet.");
// {
// Level level;
// RCP<FactoryManager> factoryHandler = rcp(new FactoryManager());
// factoryHandler->SetKokkosRefactor(false);
// level.SetFactoryManager(factoryHandler);
// level.SetLevelID(0);
// level.Set("A", coarseCompMat);
// level.setlib(compX->getMap()->UnderlyingLib());
// }
{
/*
1. Create DirectSolver by calling its constructor
2. Call DirectSolver::Copy() to obtain Amesos/Amesos2 object wrapped into a SmootherPrototype
3. Call Setup() and Apply() on the SmootherPrototype
*/
}
// {
// using DirectSolver = MueLu::DirectSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
// using FactoryManager = MueLu::FactoryManager<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
// using Hierarchy = MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
// using SmootherPrototype = MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
// using SmootherFactory = MueLu::SmootherFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
//
// RCP<Hierarchy> H = rcp(new Hierarchy(coarseCompMat));
//
// RCP<SmootherPrototype> coarseProto = rcp(new DirectSolver("Klu"));
// RCP<SmootherFactory> coarseSolveFact = rcp(new SmootherFactory(coarseProto));
//
// FactoryManager M;
// M.SetFactory("CoarseSolver", coarseSolveFact);
//
// H->Setup(M, 0, 1);
// H->Iterate(*compRhs, *compX, 1);
// }
// Replace non-continuos maps by continuous maps
// compX->replaceMap(contigRowMap);
// compRhs->replaceMap(contigRowMap);
// coarseCompMat->replaceRowMap(*contigRowMap);
// TEUCHOS_ASSERT(err==0);
// err = coarseCompMat->ReplaceColMap(*contigColMap);
// TEUCHOS_ASSERT(err==0);
// err = coarseCompMat->ExpertStaticFillComplete(*contigRowMap, *contigRowMap);
// TEUCHOS_ASSERT(err==0);
//
// // create a linear problem object
// Epetra_LinearProblem problem(coarseCompMat.get(), &(*compX), &(*compRhs));
//
// // Direct solver
// {
// Teuchos::ParameterList pList;
// pList.set("PrintTiming",true);
// pList.set("PrintStatus",true);
// pList.set("MaxProcs", coarseCompMat->Comm().NumProc());
//
// Amesos Factory;
// Amesos_BaseSolver* solver = Factory.Create("Amesos_Umfpack", problem);
// TEUCHOS_ASSERT(solver!=NULL);
//
// solver->SetParameters(pList);
// solver->SetUseTranspose(false);
//
// solver->SymbolicFactorization();
// solver->NumericFactorization();
// solver->Solve();
// }
//
// // Replace maps with the original non-continuous maps
// compX->ReplaceMap(*noncontigColMap);
// compRhs->ReplaceMap(*noncontigRowMap);
}
return;
......
packages/muelu/research/regionMG/examples/structured/Driver_Structured_Regions.cpp
View file @ 852d7e71
......@@ -116,6 +116,11 @@
#include <MueLu_CreateXpetraPreconditioner.hpp>
#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_AMESOS2)
#include <Amesos2_config.h>
#include <Amesos2.hpp>
#endif
// Region MG headers
#include "SetupRegionHierarchy_def.hpp"
......@@ -1169,6 +1174,8 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
// create region coordinates vector
regionCoordinates[0] = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(revisedRowMapPerGrp[0],
numDimensions);
using Tpetra_CrsMatrix = Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
using Tpetra_MultiVector = Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
/* Stuff for multi-level algorithm
*
......@@ -1190,6 +1197,7 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
Array<std::vector<RCP<Import> > > regRowImporters; // regional row importers on each level
Array<std::vector<RCP<Vector> > > regInterfaceScalings; // regional interface scaling factors on each level
Teuchos::RCP<Matrix> coarseCompOp = Teuchos::null;
RCP<Amesos2::Solver<Tpetra_CrsMatrix, Tpetra_MultiVector> > coarseCompositeDirectSolver = Teuchos::null;
// Create multigrid hierarchy
......@@ -1219,7 +1227,8 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
regInterfaceScalings,
coarseCompOp,
maxRegPerGID,
compositeToRegionLIDs());
compositeToRegionLIDs(),
coarseCompositeDirectSolver);
comm->barrier();
......@@ -1326,7 +1335,7 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
vCycle(0, numLevels, maxFineIter, maxCoarseIter, omega, maxRegPerProc, regX, regB,
regMatrices,
regProlong, compRowMaps, quasiRegRowMaps, regRowMaps, regRowImporters,
regInterfaceScalings, coarseCompOp);
regInterfaceScalings, coarseCompOp, coarseCompositeDirectSolver);
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment