diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_TekoSmoother_decl.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_TekoSmoother_decl.hpp index b30515984cff..e83daca5e696 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_TekoSmoother_decl.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_TekoSmoother_decl.hpp @@ -174,7 +174,7 @@ class TekoSmoother : public SmootherPrototype< RCP GetValidParameterList() const { RCP validParamList = rcp(new ParameterList()); - validParamList->set >("A", null, "Generating factory of the matrix A"); + validParamList->set>("A", null, "Generating factory of the matrix A"); validParamList->set("Inverse Type", "", "Name of parameter list within 'Teko parameters' containing the Teko smoother parameters."); return validParamList; @@ -200,7 +200,7 @@ class TekoSmoother : public SmootherPrototype< this->GetOStream(Warnings0) << "MueLu::TekoSmoother::Setup(): Setup() has already been called"; // extract blocked operator A from current level - A_ = Factory::Get >(currentLevel, "A"); // A needed for extracting map extractors + A_ = Factory::Get>(currentLevel, "A"); // A needed for extracting map extractors bA_ = Teuchos::rcp_dynamic_cast(A_); TEUCHOS_TEST_FOR_EXCEPTION(bA_.is_null(), Exceptions::BadCast, "MueLu::TekoSmoother::Build: input matrix A is not of type BlockedCrsMatrix."); @@ -209,7 +209,7 @@ class TekoSmoother : public SmootherPrototype< TEUCHOS_TEST_FOR_EXCEPTION(bThyOp_.is_null(), Exceptions::BadCast, "MueLu::TekoSmoother::Build: Could not extract thyra operator from BlockedCrsMatrix."); - Teuchos::RCP > thyOp = Teuchos::rcp_dynamic_cast >(bThyOp_); + Teuchos::RCP> thyOp = Teuchos::rcp_dynamic_cast>(bThyOp_); TEUCHOS_TEST_FOR_EXCEPTION(thyOp.is_null(), Exceptions::BadCast, "MueLu::TekoSmoother::Build: Downcast of Thyra::BlockedLinearOpBase to Teko::LinearOp failed."); @@ -237,51 +237,64 @@ class TekoSmoother : public SmootherPrototype< @param X initial guess @param B right-hand side - @param InitialGuessIsZero This option has no effect. + @param InitialGuessIsZero whether the initial guess is zero */ - void Apply(MultiVector &X, const MultiVector &B, bool /* InitialGuessIsZero */ = false) const { - TEUCHOS_TEST_FOR_EXCEPTION(this->IsSetup() == false, Exceptions::RuntimeError, + void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero = false) const { + TEUCHOS_TEST_FOR_EXCEPTION(!this->IsSetup(), Exceptions::RuntimeError, "MueLu::TekoSmoother::Apply(): Setup() has not been called"); - Teuchos::RCP > comm = X.getMap()->getComm(); - - Teuchos::RCP rgMapExtractor = bA_->getRangeMapExtractor(); + auto comm = X.getMap()->getComm(); + auto rgMapExtractor = bA_->getRangeMapExtractor(); TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgMapExtractor)); - // copy initial solution vector X to Ptr YY - - // create a Thyra RHS vector - Teuchos::RCP > thyB = Thyra::createMembers(Teuchos::rcp_dynamic_cast >(bThyOp_->productRange()), Teuchos::as(B.getNumVectors())); - Teuchos::RCP > thyProdB = - Teuchos::rcp_dynamic_cast >(thyB); - TEUCHOS_TEST_FOR_EXCEPTION(thyProdB.is_null(), Exceptions::BadCast, - "MueLu::TekoSmoother::Apply: Failed to cast range space to product range space."); - - // copy RHS vector B to Thyra::MultiVectorBase thyProdB - Xpetra::ThyraUtils::updateThyra(Teuchos::rcpFromRef(B), rgMapExtractor, thyProdB); - - // create a Thyra SOL vector - Teuchos::RCP > thyX = Thyra::createMembers(Teuchos::rcp_dynamic_cast >(bThyOp_->productDomain()), Teuchos::as(X.getNumVectors())); - Teuchos::RCP > thyProdX = - Teuchos::rcp_dynamic_cast >(thyX); - TEUCHOS_TEST_FOR_EXCEPTION(thyProdX.is_null(), Exceptions::BadCast, - "MueLu::TekoSmoother::Apply: Failed to cast domain space to product domain space."); - - // copy RHS vector X to Thyra::MultiVectorBase thyProdX - Xpetra::ThyraUtils::updateThyra(Teuchos::rcpFromRef(X), rgMapExtractor, thyProdX); - - inverseOp_->apply( - Thyra::NOTRANS, - *thyB, // const MultiVectorBase &X, - thyX.ptr(), // const Ptr > &Y, - 1.0, - 0.0); - - // copy back content of Ptr thyX into X - Teuchos::RCP > XX = - Xpetra::ThyraUtils::toXpetra(thyX, comm); - - X.update(Teuchos::ScalarTraits::one(), *XX, Teuchos::ScalarTraits::zero()); + using STS = Teuchos::ScalarTraits; + + auto createProductMultiVector = [](const Teuchos::RCP> &space, + int numVecs, + const char *castErrorMsg) { + auto mv = Thyra::createMembers(space, numVecs); + auto prodMv = Teuchos::rcp_dynamic_cast>(mv); + TEUCHOS_TEST_FOR_EXCEPTION(prodMv.is_null(), Exceptions::BadCast, castErrorMsg); + return std::make_pair(mv, prodMv); + }; + + auto solveWithThyra = [&](const Teuchos::RCP> &rhs, + bool initializeFromX, + const typename STS::magnitudeType alphaX) { + auto [thyB, thyProdB] = createProductMultiVector( + Teuchos::rcp_dynamic_cast>(bThyOp_->productRange()), + Teuchos::as(rhs->getNumVectors()), + "MueLu::TekoSmoother::Apply: Failed to cast range space to product range space."); + + Xpetra::ThyraUtils::updateThyra(rhs, rgMapExtractor, thyProdB); + + auto [thyX, thyProdX] = createProductMultiVector( + Teuchos::rcp_dynamic_cast>(bThyOp_->productDomain()), + Teuchos::as(X.getNumVectors()), + "MueLu::TekoSmoother::Apply: Failed to cast domain space to product domain space."); + + if (initializeFromX) { + Xpetra::ThyraUtils::updateThyra(Teuchos::rcpFromRef(X), rgMapExtractor, thyProdX); + } + + inverseOp_->apply( + Thyra::NOTRANS, + *thyB, + thyX.ptr(), + STS::one(), + STS::zero()); + + auto XX = Xpetra::ThyraUtils::toXpetra(thyX, comm); + X.update(STS::one(), *XX, alphaX); + }; + + if (InitialGuessIsZero) { + solveWithThyra(Teuchos::rcpFromRef(B), /*initializeFromX=*/true, STS::zero()); + return; + } + + auto residual = Utilities::Residual(*A_, X, B); + solveWithThyra(residual, /*initializeFromX=*/false, STS::one()); } //@} @@ -325,7 +338,7 @@ class TekoSmoother : public SmootherPrototype< //! block operator RCP A_; // < ! internal blocked operator "A" generated by AFact_ RCP bA_; - RCP > bThyOp_; + RCP> bThyOp_; //! Teko parameters RCP tekoParams_; // < ! parameter list containing Teko parameters. These parameters are not administrated by the factory and not validated. diff --git a/packages/muelu/test/unit_tests/Smoothers/BlockedSmoother.cpp b/packages/muelu/test/unit_tests/Smoothers/BlockedSmoother.cpp index f74608c0aa6a..3d1943fbd995 100644 --- a/packages/muelu/test/unit_tests/Smoothers/BlockedSmoother.cpp +++ b/packages/muelu/test/unit_tests/Smoothers/BlockedSmoother.cpp @@ -30,6 +30,10 @@ #include #include +#ifdef HAVE_MUELU_TEKO +#include +#endif + namespace MueLuTests { // this namespace already has: #include "MueLu_UseShortNames.hpp" @@ -6774,6 +6778,176 @@ TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, SplitReorder, Scalar, LocalOr } // end UseTpetra } +#ifdef HAVE_MUELU_TEKO +TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, Teko_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) { +#include + MUELU_TESTING_SET_OSTREAM; + MUELU_TESTING_LIMIT_SCOPE(Scalar, GlobalOrdinal, Node); + + using Teuchos::ParameterList; + + MUELU_TEST_ONLY_FOR(Xpetra::UseTpetra) { + if (!std::is_same::value || !std::is_same::value) { + out << "Skipping Teko test: TekoSmoother is only available for Scalar=Teko::ST, LocalOrdinal=Teko::LO." << std::endl; + return; + } + + if (!std::is_same::value) { + out << "Skipping Teko test: TekoSmoother is only available for Node=Teko::NT." << std::endl; + return; + } + + RCP > comm = Parameters::getDefaultComm(); + + // Teko wants a blocked Thyra-capable operator + int noBlocks = 2; + Teuchos::RCP bop = CreateBlockDiagonalExampleMatrixThyra(noBlocks, *comm); + Teuchos::RCP Aconst = Teuchos::rcp_dynamic_cast(bop); + Teuchos::RCP A = Teuchos::rcp_const_cast(Aconst); + + Level level; + TestHelpers::TestFactory::createSingleLevelHierarchy(level); + level.Set("A", A); + + ////////////////////////////////////////////////////////////////////// + // Teko smoother + RCP smootherPrototype = + rcp(new TekoSmoother()); + + smootherPrototype->SetParameter("Inverse Type", Teuchos::ParameterEntry(std::string("TekoDiag"))); + + RCP tekoParams = rcp(new ParameterList()); + { + ParameterList& tekoDiag = tekoParams->sublist("TekoDiag"); + tekoDiag.set("Type", "Block Jacobi"); + } + + smootherPrototype->SetTekoParameters(tekoParams); + + RCP smootherFact = rcp(new SmootherFactory(smootherPrototype)); + + FactoryManager M; + M.SetFactory("Smoother", smootherFact); + + MueLu::SetFactoryManager SFM(Teuchos::rcpFromRef(level), Teuchos::rcpFromRef(M)); + + level.Request("Smoother", smootherFact.get()); + level.Request("PreSmoother", smootherFact.get()); + level.Request("PostSmoother", smootherFact.get()); + + smootherFact->Build(level); + + level.print(std::cout, Teuchos::VERB_EXTREME); + + RCP tekoSmoother = level.Get >("PreSmoother", smootherFact.get()); + + RCP X = MultiVectorFactory::Build(A->getDomainMap(), 1); + RCP RHS = MultiVectorFactory::Build(A->getRangeMap(), 1); + + // Random exact solution + X->setSeed(846930886); + X->randomize(); + + typedef typename Teuchos::ScalarTraits::magnitudeType magnitude_type; + + // Normalize X + Array norms(1); + X->norm2(norms); + X->scale(1 / norms[0]); + + // Build RHS = A * X_exact + A->apply(*X, *RHS, Teuchos::NO_TRANS, (SC)1.0, (SC)0.0); + + // Reset X to 0 + X->putScalar((SC)0.0); + + RHS->norm2(norms); + out << "||RHS|| = " << std::setiosflags(std::ios::fixed) + << std::setprecision(10) << norms[0] << std::endl; + + magnitude_type tol = 100. * Teuchos::ScalarTraits::eps(); + + // + // Case 1: zero initial guess + // + out << "solve with zero initial guess" << std::endl; + Teuchos::Array initialNorms(1); + X->norm2(initialNorms); + out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) + << std::setprecision(10) << initialNorms[0] << std::endl; + + tekoSmoother->Apply(*X, *RHS, true); + + Teuchos::Array finalNorms(1); + X->norm2(finalNorms); + Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*A, *X, *RHS); + out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) + << std::setprecision(20) << residualNorm1[0] << std::endl; + out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) + << std::setprecision(10) << finalNorms[0] << std::endl; + + TEUCHOS_TEST_COMPARE(residualNorm1[0], <, tol, out, success); + TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + + // + // Case 2: InitialGuessIsZero=true with nonzero/garbage X + // + out << "solve with zero initial guess, and unreliable nonzeroed vector X" << std::endl; + X->randomize(); + X->norm2(initialNorms); + out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) + << std::setprecision(10) << initialNorms[0] << std::endl; + + tekoSmoother->Apply(*X, *RHS, true); + + X->norm2(finalNorms); + Teuchos::Array residualNorm2 = Utilities::ResidualNorm(*A, *X, *RHS); + out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) + << std::setprecision(20) << residualNorm2[0] << std::endl; + out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) + << std::setprecision(10) << finalNorms[0] << std::endl; + + TEUCHOS_TEST_COMPARE(residualNorm2[0], <, tol, out, success); + TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + + // + // Case 3: random initial guess, InitialGuessIsZero=false + // + out << "solve with random initial guess" << std::endl; + X->randomize(); + X->norm2(initialNorms); + out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) + << std::setprecision(10) << initialNorms[0] << std::endl; + + tekoSmoother->Apply(*X, *RHS, false); + + X->norm2(finalNorms); + Teuchos::Array residualNorm3 = Utilities::ResidualNorm(*A, *X, *RHS); + out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) + << std::setprecision(20) << residualNorm3[0] << std::endl; + out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) + << std::setprecision(10) << finalNorms[0] << std::endl; + + TEUCHOS_TEST_COMPARE(residualNorm3[0], <, tol, out, success); + TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + + if (comm->getSize() == 1) { + TEST_EQUALITY(residualNorm1[0] == residualNorm2[0], true); + TEST_EQUALITY(residualNorm1[0] != residualNorm3[0], true); + } else { + out << "Pass/Fail is only checked in serial." << std::endl; + } + } // end UseTpetra +} +#endif // HAVE_MUELU_TEKO + +#ifdef HAVE_MUELU_TEKO +#define MUELU_ETI_GROUP_TEKO(SC, LO, GO, NO) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother, Teko_Setup_Apply, SC, LO, GO, NO) +#else +#define MUELU_ETI_GROUP_TEKO(SC, LO, GO, NO) +#endif + #define MUELU_ETI_GROUP(SC, LO, GO, NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother, Jacobi_Setup_Apply, SC, LO, GO, NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother, BGS_Setup_Apply, SC, LO, GO, NO) \ @@ -6812,7 +6986,8 @@ TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, SplitReorder, Scalar, LocalOr TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother, NestedI2I01II_Thyra_Indef_Setup_Apply, SC, LO, GO, NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother, NestedII01I2I_Thyra_Indef_Setup_Apply2, SC, LO, GO, NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother, NestedI0I21II_Thyra_Indef_Setup_Apply3, SC, LO, GO, NO) \ - TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother, SplitReorder, SC, LO, GO, NO) + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother, SplitReorder, SC, LO, GO, NO) \ + MUELU_ETI_GROUP_TEKO(SC, LO, GO, NO) #include