diff --git a/packages/ifpack2/src/Ifpack2_Relaxation_def.hpp b/packages/ifpack2/src/Ifpack2_Relaxation_def.hpp index fe996f1baad0..afb0253ee507 100644 --- a/packages/ifpack2/src/Ifpack2_Relaxation_def.hpp +++ b/packages/ifpack2/src/Ifpack2_Relaxation_def.hpp @@ -2658,6 +2658,8 @@ std::string Relaxation::description () const else { os << "INVALID"; } + if(hasBlockCrsMatrix_) + os<<", BlockCrs"; os << ", " << "sweeps: " << NumSweeps_ << ", " << "damping factor: " << DampingFactor_ << ", "; diff --git a/packages/muelu/src/Smoothers/MueLu_Ifpack2Smoother_def.hpp b/packages/muelu/src/Smoothers/MueLu_Ifpack2Smoother_def.hpp index b63515b9aa22..833a1d74e815 100644 --- a/packages/muelu/src/Smoothers/MueLu_Ifpack2Smoother_def.hpp +++ b/packages/muelu/src/Smoothers/MueLu_Ifpack2Smoother_def.hpp @@ -65,16 +65,21 @@ #include #include #include +#include #include #include #include +#include + #include "MueLu_Ifpack2Smoother_decl.hpp" #include "MueLu_Level.hpp" #include "MueLu_FactoryManagerBase.hpp" #include "MueLu_Utilities.hpp" #include "MueLu_Monitor.hpp" + + #ifdef HAVE_MUELU_INTREPID2 #include "MueLu_IntrepidPCoarsenFactory_decl.hpp" #include "MueLu_IntrepidPCoarsenFactory_def.hpp" @@ -192,6 +197,31 @@ namespace MueLu { void Ifpack2Smoother::Setup(Level& currentLevel) { FactoryMonitor m(*this, "Setup Smoother", currentLevel); A_ = Factory::Get< RCP >(currentLevel, "A"); + ParameterList& paramList = const_cast(this->GetParameterList()); + + // If the user asked us to convert the matrix into BlockCrsMatrix form, we do that now + if(paramList.isParameter("smoother: use blockcrsmatrix storage") && paramList.get("smoother: use blockcrsmatrix storage") && A_->GetFixedBlockSize()) { + // NOTE: Don't think you can move this out of the if block. You can't. The test MueLu_MeshTyingBlocked_SimpleSmoother_2dof_medium_MPI_1 will fail + int blocksize = A_->GetFixedBlockSize(); + using TpetraBlockCrsMatrix = Xpetra::TpetraBlockCrsMatrix; + RCP AcrsWrap = rcp_dynamic_cast(A_); + if(AcrsWrap.is_null()) + throw std::runtime_error("Ifpack2Smoother: Cannot convert matrix A to CrsMatrixWrap object."); + RCP Acrs = AcrsWrap->getCrsMatrix(); + if(Acrs.is_null()) + throw std::runtime_error("Ifpack2Smoother: Cannot extract CrsMatrix from matrix A."); + RCP At = rcp_dynamic_cast(Acrs); + if(At.is_null()) + throw std::runtime_error("Ifpack2Smoother: Cannot extract TpetraCrsMatrix from matrix A."); + + RCP > blockCrs = Tpetra::convertToBlockCrsMatrix(*At->getTpetra_CrsMatrix(),blocksize); + RCP blockCrs_as_crs = rcp(new TpetraBlockCrsMatrix(blockCrs)); + RCP blockWrap = rcp(new CrsMatrixWrap(blockCrs_as_crs)); + A_ = blockWrap; + this->GetOStream(Statistics0) << "Ifpack2Smoother: Using BlockCrsMatrix storage with blocksize "< + MUELU_TESTING_LIMIT_SCOPE(Scalar,GlobalOrdinal,Node); + MUELU_TESTING_SET_OSTREAM; + MUELU_TEST_ONLY_FOR(Xpetra::UseTpetra) { + +# if !defined(HAVE_MUELU_AMESOS) || !defined(HAVE_MUELU_IFPACK) + MUELU_TESTING_DO_NOT_TEST(Xpetra::UseEpetra, "Amesos, Ifpack"); +# endif +# if !defined(HAVE_MUELU_AMESOS2) || !defined(HAVE_MUELU_IFPACK2) + MUELU_TESTING_DO_NOT_TEST(Xpetra::UseTpetra, "Amesos2, Ifpack2"); +# endif + + RCP > comm = TestHelpers::Parameters::getDefaultComm(); + GO nx = 10*comm->getSize(); + Teuchos::ParameterList galeriList, ifpack2Params; + galeriList.set("nx", nx); + RCP A = TestHelpers::TpetraTestFactory::BuildBlockMatrixAsPoint(galeriList,Xpetra::UseTpetra); + + ifpack2Params.set("smoother: use blockcrsmatrix storage",true); + + // Multigrid Hierarchy + Hierarchy H(A); + H.setVerbLevel(Teuchos::VERB_HIGH); + H.SetMaxCoarseSize(2); + + // Multigrid setup phase (using default parameters) + FactoryManager M0; // how to build aggregates and smoother of the first level + M0.SetKokkosRefactor(false); + M0.SetFactory("Smoother",rcp(new SmootherFactory(rcp(new TrilinosSmoother("RELAXATION", ifpack2Params, 0))))); + + FactoryManager M1; // first coarse level (Plain aggregation) + M1.SetKokkosRefactor(false); + M1.SetFactory("A", rcp(new RAPFactory())); + RCP P = rcp(new TentativePFactory()); + M1.SetFactory("P", P); + M1.SetFactory("Ptent", P); //FIXME: can it be done automatically in FactoryManager? + M1.SetFactory("Smoother",rcp(new SmootherFactory(rcp(new TrilinosSmoother("RELAXATION", ifpack2Params, 0))))); + + FactoryManager M2; // last level (SA) + M2.SetKokkosRefactor(false); + M2.SetFactory("A", rcp(new RAPFactory())); + M2.SetFactory("P", rcp(new SaPFactory())); + + bool r; // cf. bug Teuchos Bug 5214 + r = H.Setup(0, Teuchos::null, rcpFromRef(M0), rcpFromRef(M1)); TEST_EQUALITY(r, false); + r = H.Setup(1, rcpFromRef(M0), rcpFromRef(M1), rcpFromRef(M2)); TEST_EQUALITY(r, false); + r = H.Setup(2, rcpFromRef(M1), rcpFromRef(M2), Teuchos::null ); TEST_EQUALITY(r, true); + + RCP l0 = H.GetLevel(0); + RCP l1 = H.GetLevel(1); + RCP l2 = H.GetLevel(2); + + TEST_EQUALITY(l0->IsAvailable("PreSmoother", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l1->IsAvailable("PreSmoother", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l2->IsAvailable("PreSmoother", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l0->IsAvailable("PostSmoother", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l1->IsAvailable("PostSmoother", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l2->IsAvailable("PostSmoother", MueLu::NoFactory::get()), false); // direct solve + TEST_EQUALITY(l1->IsAvailable("P", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l2->IsAvailable("P", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l1->IsAvailable("R", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l2->IsAvailable("R", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l0->IsAvailable("A", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l1->IsAvailable("A", MueLu::NoFactory::get()), true); + TEST_EQUALITY(l2->IsAvailable("A", MueLu::NoFactory::get()), true); + + TEST_EQUALITY(l0->GetKeepFlag("A", MueLu::NoFactory::get()), MueLu::UserData); + TEST_EQUALITY(l0->GetKeepFlag("PreSmoother", MueLu::NoFactory::get()), MueLu::Final); + TEST_EQUALITY(l0->GetKeepFlag("PostSmoother", MueLu::NoFactory::get()), MueLu::Final); + + TEST_EQUALITY(l1->GetKeepFlag("A", MueLu::NoFactory::get()), MueLu::Final); + TEST_EQUALITY(l1->GetKeepFlag("P", MueLu::NoFactory::get()), MueLu::Final); + TEST_EQUALITY(l1->GetKeepFlag("R", MueLu::NoFactory::get()), MueLu::Final); + TEST_EQUALITY(l1->GetKeepFlag("PreSmoother", MueLu::NoFactory::get()), MueLu::Final); + TEST_EQUALITY(l1->GetKeepFlag("PostSmoother", MueLu::NoFactory::get()), MueLu::Final); + + TEST_EQUALITY(l2->GetKeepFlag("A", MueLu::NoFactory::get()), MueLu::Final); + TEST_EQUALITY(l2->GetKeepFlag("P", MueLu::NoFactory::get()), MueLu::Final); + TEST_EQUALITY(l2->GetKeepFlag("R", MueLu::NoFactory::get()), MueLu::Final); + TEST_EQUALITY(l2->GetKeepFlag("PreSmoother", MueLu::NoFactory::get()), MueLu::Final); + // TEST_EQUALITY(l2->GetKeepFlag("PostSmoother", MueLu::NoFactory::get()), MueLu::Final); // direct solve + + RCP RHS = MultiVectorFactory::Build(A->getRowMap(), 1); + RCP X = MultiVectorFactory::Build(A->getRowMap(), 1); + RHS->setSeed(846930886); + RHS->randomize(); + + X->putScalar( (Scalar) 0.0); + + int iterations=10; + H.Iterate(*RHS, *X, iterations); + } + } + + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Hierarchy, SetupHierarchy3level_NoPreSmooth, Scalar, LocalOrdinal, GlobalOrdinal, Node) { # include @@ -1475,6 +1572,7 @@ namespace MueLuTests { TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Hierarchy, SetupHierarchy1levelv2, Scalar, LO, GO, Node) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Hierarchy, SetupHierarchy2level, Scalar, LO, GO, Node) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Hierarchy, SetupHierarchy3level, Scalar, LO, GO, Node) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Hierarchy, SetupHierarchy3level_BlockSmooth, Scalar, LO, GO, Node) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Hierarchy, SetupHierarchy3level_NoPreSmooth, Scalar, LO, GO, Node) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Hierarchy, SetupHierarchy3levelFacManagers, Scalar, LO, GO, Node) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Hierarchy, SetupHierarchyTestBreakCondition, Scalar, LO, GO, Node) \ diff --git a/packages/muelu/test/unit_tests/MueLu_TestHelpers.hpp b/packages/muelu/test/unit_tests/MueLu_TestHelpers.hpp index 3c9fdf105327..7b4f041f2a9f 100644 --- a/packages/muelu/test/unit_tests/MueLu_TestHelpers.hpp +++ b/packages/muelu/test/unit_tests/MueLu_TestHelpers.hpp @@ -735,58 +735,6 @@ namespace MueLuTests { return bop; } - // Create a matrix as specified by parameter list options - /*static RCP BuildBlockMatrix(Teuchos::ParameterList &matrixList, Xpetra::UnderlyingLib lib=Xpetra::NotSpecified) { - RCP > comm = TestHelpers::Parameters::getDefaultComm(); - RCP Op; - - if (lib == Xpetra::NotSpecified) - lib = TestHelpers::Parameters::getLib(); - - // This only works for Tpetra - if (lib!=Xpetra::UseTpetra) return Op; - -#if defined(HAVE_MUELU_TPETRA) -#ifdef HAVE_MUELU_BROKEN_TESTS - // Thanks for the code, Travis! - - // Make the graph - RCP FirstMatrix = BuildMatrix(matrixList,lib); - RCP > Graph = FirstMatrix->getCrsGraph(); - - int blocksize = 3; - RCP > TGraph = rcp_dynamic_cast >(Graph); - RCP > TTGraph = TGraph->getTpetra_CrsGraph(); - - RCP > bcrsmatrix = rcp(new Tpetra::BlockCrsMatrix (*TTGraph, blocksize)); - - const Tpetra::Map& meshRowMap = *bcrsmatrix->getRowMap(); - const Scalar zero = Teuchos::ScalarTraits::zero(); - const Scalar one = Teuchos::ScalarTraits::one(); - const Scalar two = one+one; - const Scalar three = two+one; - - Teuchos::Array basematrix(blocksize*blocksize, zero); - basematrix[0] = two; - basematrix[2] = three; - basematrix[3] = three; - basematrix[4] = two; - basematrix[7] = three; - basematrix[8] = two; - Teuchos::Array lclColInds(1); - for (LocalOrdinal lclRowInd = meshRowMap.getMinLocalIndex (); lclRowInd <= meshRowMap.getMaxLocalIndex(); ++lclRowInd) { - lclColInds[0] = lclRowInd; - bcrsmatrix->replaceLocalValues(lclRowInd, lclColInds.getRawPtr(), &basematrix[0], 1); - } - - RCP > temp = rcp(new Xpetra::TpetraBlockCrsMatrix(bcrsmatrix)); - Op = rcp(new Xpetra::CrsMatrixWrap(temp)); -#endif -#endif - return Op; - } // BuildMatrix()*/ - - // Needed to initialize correctly a level used for testing SingleLevel factory Build() methods. // This method initializes LevelID and linked list of level static void createSingleLevelHierarchy(Level& currentLevel) { @@ -827,7 +775,7 @@ namespace MueLuTests { // We put this into an extra helper class as we need partial specializations and // do not want to introduce partial specializations for the full TestFactory class // - // The BuildBlockMatrix is only available with Teptra. However, if both Epetra + // The BuildBlockMatrix is only available with Tpetra. However, if both Epetra // and Tpetra are enabled it may be that Tpetra is not instantiated on either // GO=int/long long and/or Node=Serial/OpenMP. We need partial specializations // with an empty BuildBlockMatrix routine for all instantiations Teptra is not @@ -886,6 +834,68 @@ namespace MueLuTests { return Op; } // BuildBlockMatrix() + // Create a matrix as specified by parameter list options + static RCP BuildBlockMatrixAsPoint(Teuchos::ParameterList &matrixList, Xpetra::UnderlyingLib lib) { + RCP > comm = TestHelpers::Parameters::getDefaultComm(); + GO GO_INVALID = Teuchos::OrdinalTraits::invalid(); + RCP Op; + + if (lib == Xpetra::NotSpecified) + lib = TestHelpers::Parameters::getLib(); + + // This only works for Tpetra + if (lib!=Xpetra::UseTpetra) return Op; + +#if defined(HAVE_MUELU_TPETRA) + // Make the base graph + RCP old_matrix = TestHelpers::TestFactory::BuildMatrix(matrixList,lib); + RCP old_graph = old_matrix->getCrsGraph(); + RCP old_rowmap = old_graph->getRowMap(); + RCP old_colmap = old_graph->getColMap(); + int blocksize = 3; + + // Block Map + LO orig_num_rows = (LO) old_graph->getRowMap()->getNodeNumElements(); + Teuchos::Array owned_rows(blocksize*orig_num_rows); + for(LO i=0; igetGlobalElement(i); + for(int j=0; j new_map = Xpetra::MapFactory::Build(lib,GO_INVALID,owned_rows(),0,comm); + if(new_map.is_null()) throw std::runtime_error("BuildBlockMatrixAsPoint: Map constructor failed"); + + + // Block Graph / Matrix + RCP new_matrix = Xpetra::CrsMatrixFactory::Build(new_map,blocksize*old_graph->getNodeMaxNumRowEntries()); + if(new_matrix.is_null()) throw std::runtime_error("BuildBlockMatrixAsPoint: Matrix constructor failed"); + for(LO i=0; i old_indices; + Teuchos::ArrayView old_values; + Teuchos::Array new_indices(1); + Teuchos::Array new_values(1); + old_matrix->getLocalRowView(i,old_indices,old_values); + for(int ii=0; iigetGlobalElement(i*blocksize+ii); + for(LO j=0; j<(LO)old_indices.size(); j++) { + for(int jj=0; jjgetGlobalElement(old_indices[j]) * blocksize + jj; + new_values[0] = old_values[j] * (SC)( (ii == jj && i == old_indices[j] ) ? blocksize*blocksize : 1 ); + new_matrix->insertGlobalValues(GRID,new_indices(),new_values); + } + } + } + } + new_matrix->fillComplete(); + Op = rcp(new CrsMatrixWrap(new_matrix)); + if(new_map.is_null()) throw std::runtime_error("BuildBlockMatrixAsPoint: CrsMatrixWrap constructor failed"); + Op->SetFixedBlockSize(blocksize); +#endif + return Op; + } // BuildBlockMatrixAsPoint() + + private: TpetraTestFactory() {} // static class @@ -902,6 +912,7 @@ namespace MueLuTests { #include "MueLu_UseShortNames.hpp" public: static RCP BuildBlockMatrix(Teuchos::ParameterList &matrixList, Xpetra::UnderlyingLib lib) { return Teuchos::null; } + static RCP BuildBlockMatrixAsPoint(Teuchos::ParameterList &matrixList, Xpetra::UnderlyingLib lib) { return Teuchos::null; } private: TpetraTestFactory() {} // static class }; // class TpetraTestFactory @@ -915,6 +926,7 @@ namespace MueLuTests { #include "MueLu_UseShortNames.hpp" public: static RCP BuildBlockMatrix(Teuchos::ParameterList &matrixList, Xpetra::UnderlyingLib lib) { return Teuchos::null; } + static RCP BuildBlockMatrixAsPoint(Teuchos::ParameterList &matrixList, Xpetra::UnderlyingLib lib) { return Teuchos::null; } private: TpetraTestFactory() {} // static class }; // class TpetraTestFactory @@ -930,6 +942,7 @@ namespace MueLuTests { #include "MueLu_UseShortNames.hpp" public: static RCP BuildBlockMatrix(Teuchos::ParameterList &matrixList, Xpetra::UnderlyingLib lib) { return Teuchos::null; } + static RCP BuildBlockMatrixAsPoint(Teuchos::ParameterList &matrixList, Xpetra::UnderlyingLib lib) { return Teuchos::null; } private: TpetraTestFactory() {} // static class }; // class TpetraTestFactory diff --git a/packages/muelu/test/unit_tests/Smoothers/Ifpack2Smoother.cpp b/packages/muelu/test/unit_tests/Smoothers/Ifpack2Smoother.cpp index 1112964178a0..658c273febc1 100644 --- a/packages/muelu/test/unit_tests/Smoothers/Ifpack2Smoother.cpp +++ b/packages/muelu/test/unit_tests/Smoothers/Ifpack2Smoother.cpp @@ -457,6 +457,31 @@ namespace MueLuTests { } // banded + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Ifpack2Smoother, BlockCrsMatrix_Relaxation, Scalar, LocalOrdinal, GlobalOrdinal, Node) + { +# include + MUELU_TESTING_SET_OSTREAM; + MUELU_TESTING_LIMIT_SCOPE(Scalar,GlobalOrdinal,Node); + + MUELU_TEST_ONLY_FOR(Xpetra::UseTpetra) { + Teuchos::ParameterList matrixParams, ifpack2Params; + + matrixParams.set("matrixType","Laplace1D"); + matrixParams.set("nx",(GlobalOrdinal)20);// needs to be even + + RCP A = TestHelpers::TpetraTestFactory::BuildBlockMatrixAsPoint(matrixParams,Xpetra::UseTpetra); + ifpack2Params.set("smoother: use blockcrsmatrix storage",true); + + Ifpack2Smoother smoother("RELAXATION",ifpack2Params); + + Level level; TestHelpers::TestFactory::createSingleLevelHierarchy(level); + level.Set("A", A); + smoother.Setup(level); + + TEST_EQUALITY(1,1); + } + } + #define MUELU_ETI_GROUP(SC,LO,GO,NO) \ @@ -468,7 +493,8 @@ namespace MueLuTests { TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Ifpack2Smoother,ILU_TwoSweeps,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Ifpack2Smoother,BandedRelaxation,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Ifpack2Smoother,TriDiRelaxation,SC,LO,GO,NO) \ - TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Ifpack2Smoother,BlockRelaxation_Autosize,SC,LO,GO,NO) + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Ifpack2Smoother,BlockRelaxation_Autosize,SC,LO,GO,NO) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Ifpack2Smoother,BlockCrsMatrix_Relaxation,SC,LO,GO,NO) #include diff --git a/packages/stokhos/src/CMakeLists.txt b/packages/stokhos/src/CMakeLists.txt index f8095978169d..f5a5667e6069 100644 --- a/packages/stokhos/src/CMakeLists.txt +++ b/packages/stokhos/src/CMakeLists.txt @@ -804,6 +804,7 @@ IF (Stokhos_ENABLE_Sacado) Details::packCrsMatrix Details::unpackCrsMatrixAndCombine replaceDiagonalCrsMatrix + BlockCrsMatrix_Helpers ) SET(TPETRAEXT_ETI_CLASSES MatrixMatrix