MueLu: fix issue #7411

packages/muelu/doc/UsersGuide/masterList.xml

@@ -88,6 +88,14 @@
       <visible>false</visible>
     </parameter>
 
+    <parameter>
+      <name>number of vectors</name>
+      <type>int</type>
+      <default>1</default>
+      <description>Number of columns in multivectors that are cached for Hierarchy apply phase.</description>
+      <visible>true</visible>
+    </parameter>
+
     <parameter>
       <name>problem: symmetric</name>
       <type>bool</type>

packages/muelu/doc/UsersGuide/options_general.tex

@@ -11,6 +11,8 @@
 
 \cbb{cycle type}{string}{"V"}{Multigrid cycle type. Possible values: "V", "W".}
 
+\cbb{number of vectors}{int}{1}{Number of columns in multivectors that are cached for Hierarchy apply phase.}
+
 \cbb{problem: symmetric}{bool}{true}{Symmetry of a problem. This setting affects the construction of a restrictor. If set to true, the restrictor is set to be the transpose of a prolongator. If set to false, underlying multigrid algorithm makes the decision.}
 
 \cbb{xml parameter file}{string}{""}{An XML file from which to read additional

packages/muelu/doc/UsersGuide/paramlist.tex

@@ -11,6 +11,8 @@
 
 \cbb{cycle type}{string}{"V"}{Multigrid cycle type. Possible values: "V", "W".}
 
+\cbb{number of vectors}{int}{1}{Number of columns in multivectors that are cached for Hierarchy apply phase.}
+
 \cbb{problem: symmetric}{bool}{true}{Symmetry of a problem. This setting affects the construction of a restrictor. If set to true, the restrictor is set to be the transpose of a prolongator. If set to false, underlying multigrid algorithm makes the decision.}
 
 \cbb{xml parameter file}{string}{""}{An XML file from which to read additional

packages/muelu/doc/UsersGuide/paramlist_hidden.tex

@@ -15,6 +15,8 @@
 
 \cbb{fuse prolongation and update}{bool}{false}{Fuse prolongation and update into one kernel call. Due to round-off error accumulation, this can in some cases result in slightly higher iteration counts. This only affects the solve phase.}
 
+\cbb{number of vectors}{int}{1}{Number of columns in multivectors that are cached for Hierarchy apply phase.}
+
 \cbb{problem: symmetric}{bool}{true}{Symmetry of a problem. This setting affects the construction of a restrictor. If set to true, the restrictor is set to be the transpose of a prolongator. If set to false, underlying multigrid algorithm makes the decision.}
 
 \cbb{xml parameter file}{string}{""}{An XML file from which to read additional

packages/muelu/src/Interface/MueLu_HierarchyManager.hpp

@@ -93,6 +93,7 @@ namespace MueLu {
         doPRrebalance_          (MasterList::getDefault<bool>("repartition: rebalance P and R")),
         implicitTranspose_      (MasterList::getDefault<bool>("transpose: use implicit")),
         fuseProlongationAndUpdate_ (MasterList::getDefault<bool>("fuse prolongation and update")),
+        sizeOfMultiVectors_     (MasterList::getDefault<int>("number of vectors")),
         graphOutputLevel_(-1) { }
 
     //!
@@ -249,7 +250,7 @@ namespace MueLu {
       if (!matvecParams_.is_null())
         H.SetMatvecParams(matvecParams_);
       // FIXME: Should allow specification of NumVectors on parameterlist
-      H.AllocateLevelMultiVectors(1);
+      H.AllocateLevelMultiVectors(sizeOfMultiVectors_);
       H.describe(H.GetOStream(Runtime0), verbosity_);
 
       // When we reuse hierarchy, it is necessary that we don't
@@ -316,6 +317,7 @@ namespace MueLu {
     bool                  doPRrebalance_;
     bool                  implicitTranspose_;
     bool                  fuseProlongationAndUpdate_;
+    int                   sizeOfMultiVectors_;
     int                   graphOutputLevel_;
     Teuchos::Array<int>   matricesToPrint_;
     Teuchos::Array<int>   prolongatorsToPrint_;

packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp

@@ -1942,6 +1942,11 @@ namespace MueLu {
         hieraList.remove("fuse prolongation and update");
       }
 
+      if (hieraList.isParameter("number of vectors")) {
+        this->numDesiredLevel_ = hieraList.get<int>("number of vectors");
+        hieraList.remove("number of vectors");
+      }
+
       if (hieraList.isSublist("matvec params"))
         this->matvecParams_ = Teuchos::parameterList(hieraList.sublist("matvec params"));
 

packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp

@@ -552,9 +552,9 @@ namespace MueLu {
     Levels_       .resize(levelID);
     levelManagers_.resize(levelID);
 
-    // NOTE: All reuse cases leave all of the maps the same, meaning that we do not
-    // need to reallocated the cached multivectors for Iterate().  If this were to change,
-    // we'd want to do a DeleteLevelMultiVectors() and AllocateLevelMultiVectors() here.
+    int sizeOfVecs = sizeOfAllocatedLevelMultiVectors_;
+    DeleteLevelMultiVectors();
+    AllocateLevelMultiVectors(sizeOfVecs);
 
     // since the # of levels, etc. may have changed, force re-determination of description during next call to description()
     ResetDescription();
@@ -909,7 +909,7 @@ namespace MueLu {
     // If the number of vectors is unchanged, this is a noop.
     // NOTE: We need to check against B because the tests in AllocateLevelMultiVectors
     // will fail on Stokhos Scalar types (due to the so-called 'hidden dimension')
-    const BlockedMultiVector * Bblocked = dynamic_cast<const BlockedMultiVector*>(&B);  
+    const BlockedMultiVector * Bblocked = dynamic_cast<const BlockedMultiVector*>(&B);
     if(residual_.size() > startLevel &&
        ( ( Bblocked && !Bblocked->isSameSize(*residual_[startLevel])) ||
          (!Bblocked && !residual_[startLevel]->isSameSize(B))))
@@ -1419,9 +1419,9 @@ namespace MueLu {
 
       for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
         std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
-	//	printf("[%d] CMS:   Hierarchy, adding edge (%d->%d) %d\n",rank,(int)eit->first.first,(int)eit->first.second,(int)boost_edge.second);
-	// Because xdot.py views 'Graph' as a keyword
-	if(eit->second==std::string("Graph")) boost::put("label", dp, boost_edge.first, std::string("Graph_"));
+        // printf("[%d] CMS:   Hierarchy, adding edge (%d->%d) %d\n",rank,(int)eit->first.first,(int)eit->first.second,(int)boost_edge.second);
+        // Because xdot.py views 'Graph' as a keyword
+        if(eit->second==std::string("Graph")) boost::put("label", dp, boost_edge.first, std::string("Graph_"));
         else boost::put("label", dp, boost_edge.first, eit->second);
         if (i == dumpLevel_)
           boost::put("color", dp, boost_edge.first, std::string("red"));
@@ -1525,7 +1525,7 @@ namespace MueLu {
 
   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
   void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AllocateLevelMultiVectors(int numvecs) {
-    int N = Levels_.size();  
+    int N = Levels_.size();
     if( (sizeOfAllocatedLevelMultiVectors_ == numvecs && residual_.size() == N) || numvecs<=0 ) return;
 
     // If, somehow, we changed the number of levels, delete everything first
@@ -1547,11 +1547,11 @@ namespace MueLu {
         RCP<const BlockedCrsMatrix> A_as_blocked = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(A);
         RCP<const Map> Arm = A->getRangeMap();
         RCP<const Map> Adm = A->getDomainMap();
-        if(!A_as_blocked.is_null()) { 
+        if(!A_as_blocked.is_null()) {
           Adm = A_as_blocked->getFullDomainMap();
         }
 
-        // This is zero'd by default since it is filled via an operator apply        
+        // This is zero'd by default since it is filled via an operator apply
         residual_[i] = MultiVectorFactory::Build(Arm, numvecs, true);
         correction_[i] = MultiVectorFactory::Build(Adm, numvecs, false);
       }
@@ -1570,7 +1570,7 @@ namespace MueLu {
         RCP<const Import> importer;
         if(Levels_[i+1]->IsAvailable("Importer"))
           importer = Levels_[i+1]->template Get< RCP<const Import> >("Importer");
-        if (doPRrebalance_ || importer.is_null()) 
+        if (doPRrebalance_ || importer.is_null())
           coarseX_[i] = MultiVectorFactory::Build(coarseRhs_[i]->getMap(),numvecs,true);
         else {
           coarseImport_[i] = MultiVectorFactory::Build(importer->getTargetMap(), numvecs,false);

packages/muelu/src/MueCentral/MueLu_MasterList.cpp

@@ -124,6 +124,7 @@ namespace MueLu {
     if (name == "max levels") { ss << "<Parameter name=\"max levels\" type=\"int\" value=" << value << "/>"; return ss.str(); }      
     if (name == "coarse grid correction scaling factor") { ss << "<Parameter name=\"coarse grid correction scaling factor\" type=\"double\" value=" << value << "/>"; return ss.str(); }      
     if (name == "fuse prolongation and update") { ss << "<Parameter name=\"fuse prolongation and update\" type=\"bool\" value=" << value << "/>"; return ss.str(); }      
+    if (name == "number of vectors") { ss << "<Parameter name=\"number of vectors\" type=\"int\" value=" << value << "/>"; return ss.str(); }      
     if (name == "problem: symmetric") { ss << "<Parameter name=\"problem: symmetric\" type=\"bool\" value=" << value << "/>"; return ss.str(); }      
     if (name == "hierarchy label") { ss << "<Parameter name=\"hierarchy label\" type=\"string\" value=" << value << "/>"; return ss.str(); }      
     if (name == "aggregation: drop tol") { ss << "<Parameter name=\"aggregation: drop tol\" type=\"double\" value=" << value << "/>"; return ss.str(); }      
@@ -168,6 +169,7 @@ namespace MueLu {
   "<Parameter name=\"cycle type\" type=\"string\" value=\"V\"/>"
   "<Parameter name=\"coarse grid correction scaling factor\" type=\"double\" value=\"1.0\"/>"
   "<Parameter name=\"fuse prolongation and update\" type=\"bool\" value=\"false\"/>"
+  "<Parameter name=\"number of vectors\" type=\"int\" value=\"1\"/>"
   "<Parameter name=\"problem: symmetric\" type=\"bool\" value=\"true\"/>"
   "<Parameter name=\"xml parameter file\" type=\"string\" value=\"\"/>"
   "<Parameter name=\"parameterlist: syntax\" type=\"string\" value=\"muelu\"/>"
@@ -522,6 +524,8 @@ namespace MueLu {
 
          ("fuse prolongation and update","fuse prolongation and update")
 
+         ("number of vectors","number of vectors")
+
          ("problem: symmetric","problem: symmetric")
 
          ("xml parameter file","xml parameter file")

packages/muelu/test/unit_tests/Adapters/CreatePreconditioner.cpp

@@ -599,11 +599,122 @@ namespace MueLuTests {
     }
   }
 
+  TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(PetraOperator, ReusePreconditioner2, Scalar, LocalOrdinal, GlobalOrdinal, Node)
+  {
+#   include "MueLu_UseShortNames.hpp"
+    MUELU_TESTING_SET_OSTREAM;
+    MUELU_TESTING_LIMIT_SCOPE(Scalar,GlobalOrdinal,Node);
+
+    out << "version: " << MueLu::Version() << std::endl;
+
+    using Teuchos::RCP;
+    using Teuchos::null;
+    typedef MueLu::Utilities<SC,LO,GO,NO> Utils;
+    typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType real_type;
+    typedef Xpetra::MultiVector<real_type,LO,GO,NO> dMultiVector;
+
+    Xpetra::UnderlyingLib          lib  = TestHelpers::Parameters::getLib();
+    RCP<const Teuchos::Comm<int> > comm = TestHelpers::Parameters::getDefaultComm();
+
+    Teuchos::ParameterList params;
+    params.set("aggregation: type","uncoupled");
+    params.set("aggregation: drop tol", 0.02);
+    params.set("coarse: max size", Teuchos::as<int>(500));
+
+    if (lib == Xpetra::UseTpetra) {
+      typedef Tpetra::Operator<SC,LO,GO,NO> tpetra_operator_type;
+
+      // Matrix
+      std::string matrixFile("TestMatrices/fuego0.mm");
+      RCP<const Map> rowmap        = MapFactory::Build(lib, Teuchos::as<size_t>(1500), Teuchos::as<int>(0), comm);
+      RCP<Matrix>     Op  = Xpetra::IO<SC,LO,GO,Node>::Read(matrixFile, rowmap, null, null, null);
+      RCP<const Map > map = Op->getRowMap();
+
+      // Normalized RHS
+      RCP<MultiVector> RHS1 = MultiVectorFactory::Build(Op->getRowMap(), 1);
+      RHS1->setSeed(846930886);
+      RHS1->randomize();
+      Teuchos::Array<typename Teuchos::ScalarTraits<SC>::magnitudeType> norms(1);
+      RHS1->norm2(norms);
+      RHS1->scale(1/norms[0]);
+
+      // Zero initial guess
+      RCP<MultiVector> X1   = MultiVectorFactory::Build(Op->getRowMap(), 1);
+      X1->putScalar(Teuchos::ScalarTraits<SC>::zero());
+
+      RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpA = MueLu::Utilities<SC,LO,GO,NO>::Op2NonConstTpetraCrs(Op);
+      RCP<MueLu::TpetraOperator<SC,LO,GO,NO> > tH = MueLu::CreateTpetraPreconditioner<SC,LO,GO,NO>(RCP<tpetra_operator_type>(tpA), params);
+      tH->apply(*(Utils::MV2TpetraMV(RHS1)), *(Utils::MV2NonConstTpetraMV(X1)));
+      out << "after apply, ||b-A*x||_2 = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) <<
+          Utils::ResidualNorm(*Op, *X1, *RHS1) << std::endl;
+
+      // Reuse preconditioner
+
+      matrixFile = "TestMatrices/fuego1.mm";
+      RCP<Matrix>     Op2  = Xpetra::IO<SC,LO,GO,Node>::Read(matrixFile, rowmap, null, null, null);
+      RCP<Tpetra::CrsMatrix<SC,LO,GO,NO> > tpA2 = MueLu::Utilities<SC,LO,GO,NO>::Op2NonConstTpetraCrs(Op2);
+
+      MueLu::ReuseTpetraPreconditioner(tpA2, *tH);
+
+      X1->putScalar(Teuchos::ScalarTraits<SC>::zero());
+      tH->apply(*(Utils::MV2TpetraMV(RHS1)), *(Utils::MV2NonConstTpetraMV(X1)));
+      out << "after apply, ||b-A*x||_2 = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) <<
+          Utils::ResidualNorm(*Op, *X1, *RHS1) << std::endl;
+
+    } else if (lib == Xpetra::UseEpetra) {
+#ifdef HAVE_MUELU_EPETRA
+      // Matrix
+      std::string matrixFile("TestMatrices/fuego0.mm");
+      RCP<const Map> rowmap        = MapFactory::Build(lib, Teuchos::as<size_t>(1500), Teuchos::as<int>(0), comm);
+      RCP<Matrix>     Op  = Xpetra::IO<SC,LO,GO,Node>::Read(matrixFile, rowmap, null, null, null);
+      RCP<const Map > map = Op->getRowMap();
+
+      // Normalized RHS
+      RCP<MultiVector> RHS1 = MultiVectorFactory::Build(Op->getRowMap(), 1);
+      RHS1->setSeed(846930886);
+      RHS1->randomize();
+      Teuchos::Array<typename Teuchos::ScalarTraits<SC>::magnitudeType> norms(1);
+      RHS1->norm2(norms);
+      RHS1->scale(1/norms[0]);
+
+      // Zero initial guess
+      RCP<MultiVector> X1   = MultiVectorFactory::Build(Op->getRowMap(), 1);
+      X1->putScalar(Teuchos::ScalarTraits<SC>::zero());
+
+      params.set("use kokkos refactor", false);
+      RCP<Epetra_CrsMatrix> epA = MueLu::Utilities<SC,LO,GO,NO>::Op2NonConstEpetraCrs(Op);
+      RCP<MueLu::EpetraOperator> eH = MueLu::CreateEpetraPreconditioner(epA, params);
+
+      eH->Apply(*(Utils::MV2EpetraMV(RHS1)), *(Utils::MV2NonConstEpetraMV(X1)));
+      out << "after apply, ||b-A*x||_2 = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) <<
+          Utils::ResidualNorm(*Op, *X1, *RHS1) << std::endl;
+
+      // Reuse preconditioner
+
+      matrixFile = "TestMatrices/fuego1.mm";
+      RCP<Matrix>     Op2  = Xpetra::IO<SC,LO,GO,Node>::Read(matrixFile, rowmap, null, null, null);
+      epA = MueLu::Utilities<SC,LO,GO,NO>::Op2NonConstEpetraCrs(Op);
+      MueLu::ReuseEpetraPreconditioner(epA, *eH);
+
+      X1->putScalar(Teuchos::ScalarTraits<SC>::zero());
+      eH->Apply(*(Utils::MV2EpetraMV(RHS1)), *(Utils::MV2NonConstEpetraMV(X1)));
+      out << "after apply, ||b-A*x||_2 = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) <<
+          Utils::ResidualNorm(*Op, *X1, *RHS1) << std::endl;
+#else
+      std::cout << "Skip PetraOperator::ReusePreconditioner: Epetra is not available" << std::endl;
+#endif
+
+    } else {
+      TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::InvalidArgument, "Unknown Xpetra lib");
+    }
+  }
+
 #  define MUELU_ETI_GROUP(Scalar, LO, GO, Node) \
   TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(PetraOperator, CreatePreconditioner, Scalar, LO, GO, Node) \
   TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(PetraOperator, CreatePreconditioner_XMLOnList, Scalar, LO, GO, Node) \
   TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(PetraOperator, CreatePreconditioner_PDESystem, Scalar, LO, GO, Node) \
   TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(PetraOperator, ReusePreconditioner, Scalar, LO, GO, Node) \
+  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(PetraOperator, ReusePreconditioner2, Scalar, LO, GO, Node) \
 
 #include <MueLu_ETI_4arg.hpp>
 

packages/muelu/test/unit_tests/CMakeLists.txt

@@ -91,6 +91,8 @@ TRIBITS_COPY_FILES_TO_BINARY_DIR(UnitTestsTestMatrices_cp
   SOURCE_FILES aniso2dx.mat
   SOURCE_FILES aniso2dy.mat
   SOURCE_FILES iso2d.mat
+  SOURCE_FILES fuego0.mm
+  SOURCE_FILES fuego1.mm
   DEST_DIR ${CMAKE_CURRENT_BINARY_DIR}/TestMatrices
 )