Merge 'trilinos/Trilinos:develop' (b76959a) into 'tcad-charon/Trilino…

…s:develop' (832f686).

* trilinos-develop:
  MueLu: tentative fix for issue trilinos#5962
  MueLu: fixing bug with kokkos refactor in FactoryManager, see issue trilinos#5961
  MueLu: adding unit-test for factoryManager
  mods to fix a bug in the find region code and to add more flexibility to the options that we can run and the output that we can dump.
  MueLu: Fix issue in nullspace fix
  Testing: Tpetra: enable ETI explicitly
  TrilinosCouplings: fix build errors w/o deprecated
  MueLu: remove regionMG coordinate debug output
  Sidafe's Branch Squashed in Place

cmake/ctest/drivers/rocketman/mpi_release_tpetra_deprecated_code_off_downstream_enabled.cmake

@@ -86,6 +86,7 @@ SET(Trilinos_CTEST_DO_ALL_AT_ONCE FALSE)
 # below.
 
 SET(EXTRA_CONFIGURE_OPTIONS
+  "-DTrilinos_ENABLE_EXPLICIT_INSTANTIATION=ON"
   "-DTpetra_ENABLE_DEPRECATED_CODE=OFF"
   "-DKOKKOS_ENABLE_DEPRECATED_CODE=OFF"
   "-DTPL_ENABLE_Matio=OFF"

packages/muelu/adapters/xpetra/MueLu_RefMaxwell_def.hpp

@@ -807,8 +807,19 @@ namespace MueLu {
         if (!reuse) {
           ParameterList& userParamList = precList22_.sublist("user data");
           userParamList.set<RCP<RealValuedMultiVector> >("Coordinates", Coords_);
-          // If we detected no boundary conditions, the (2,2) problem is singular
+          // If we detected no boundary conditions, the (2,2) problem is singular.
+          // Therefore, if we want to use a direct coarse solver, we need to fix up the nullspace.
+          std::string coarseType = "";
+          if (precList22_.isParameter("coarse: type")) {
+            coarseType = precList22_.get<std::string>("coarse: type");
+            // Transform string to "Abcde" notation
+            std::transform(coarseType.begin(),   coarseType.end(),   coarseType.begin(), ::tolower);
+            std::transform(coarseType.begin(), ++coarseType.begin(), coarseType.begin(), ::toupper);
+          }
           if (BCrowcount_ == 0 &&
+              (coarseType == "" ||
+               coarseType == "Klu" ||
+               coarseType == "Klu2") &&
               (!precList22_.isSublist("coarse: params") ||
                !precList22_.sublist("coarse: params").isParameter("fix nullspace")))
             precList22_.sublist("coarse: params").set("fix nullspace",true);

packages/muelu/research/mmayr/composite_to_regions/src/SetupRegionHierarchy_def.hpp

@@ -313,7 +313,6 @@ void MakeCoarseLevelMaps2(const int maxRegPerGID,
 
   using MT = typename Teuchos::ScalarTraits<SC>::magnitudeType;
 
-  RCP<const Teuchos::Comm<int> > comm = regProlong[1][0]->getRowMap()->getComm();
   const GO GO_INV = Teuchos::OrdinalTraits<GO>::invalid();
   const int numLevels = regProlong.size();
 
@@ -1048,87 +1047,95 @@ void vCycle(const int l, ///< ID of current level
 
     RCP<TimeMonitor> tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("vCycle: * - coarsest grid solve")));
 
-    // 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);
+    const bool useCoarseSmoother = false;
+    if (useCoarseSmoother) {
+      smootherApply(smootherParams[l], maxRegPerProc, fineRegX, fineRegB, regMatrices[l],
+                  regInterfaceScalings[l], compRowMaps[l],
+                  quasiRegRowMaps[l], regRowMaps[l], regRowImporters[l]);
     }
+    else {
+      // 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);
+        }
 
-    const bool useDirectSolver = coarseSolverData->get<bool>("use direct solver");
-    if (useDirectSolver)
-    {
-#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_AMESOS2)
+        regionalToComposite(fineRegB, compRhs, maxRegPerProc, quasiRegRowMaps[l],
+                            regRowImporters[l], Xpetra::ADD);
+      }
 
-      using DirectCoarseSolver = Amesos2::Solver<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >;
-      RCP<DirectCoarseSolver> coarseSolver = coarseSolverData->get<RCP<DirectCoarseSolver> >("direct solver object");
-
-      TEUCHOS_TEST_FOR_EXCEPT_MSG(coarseCompMat->getRowMap()->lib()!=Xpetra::UseTpetra,
-          "Coarse solver requires Tpetra/Amesos2 stack.");
-      TEUCHOS_ASSERT(!coarseSolver.is_null());
-
-      // using Utilities = MueLu::Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
-
-      // 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>;
-
-    //    *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);
-
-      /* Solve!
-       *
-       * Calling solve() on the coarseSolver should just do a triangular solve, since symbolic
-       * and numeric factorization are supposed to have happened during hierarchy setup.
-       * Here, we just check if they're done and print message if not.
-       *
-       * 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-continuous maps.
-       */
-      if (not coarseSolver->getStatus().symbolicFactorizationDone())
-        *fos << "Symbolic factorization should have been done during hierarchy setup, "
-            "but actually is missing. Anyway ... just do it right now." << std::endl;
-      if (not coarseSolver->getStatus().numericFactorizationDone())
-        *fos << "Numeric factorization should have been done during hierarchy setup, "
-            "but actually is missing. Anyway ... just do it right now." << std::endl;
-      coarseSolver->solve(tX.ptr(), tB.ptr());
+      const bool useDirectSolver = coarseSolverData->get<bool>("use direct solver");
+      if (useDirectSolver)
+      {
+#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_AMESOS2)
+
+        using DirectCoarseSolver = Amesos2::Solver<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >;
+        RCP<DirectCoarseSolver> coarseSolver = coarseSolverData->get<RCP<DirectCoarseSolver> >("direct solver object");
+
+        TEUCHOS_TEST_FOR_EXCEPT_MSG(coarseCompMat->getRowMap()->lib()!=Xpetra::UseTpetra,
+            "Coarse solver requires Tpetra/Amesos2 stack.");
+        TEUCHOS_ASSERT(!coarseSolver.is_null());
+
+        // using Utilities = MueLu::Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
+
+        // 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>;
+
+      //    *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);
+
+        /* Solve!
+         *
+         * Calling solve() on the coarseSolver should just do a triangular solve, since symbolic
+         * and numeric factorization are supposed to have happened during hierarchy setup.
+         * Here, we just check if they're done and print message if not.
+         *
+         * 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-continuous maps.
+         */
+        if (not coarseSolver->getStatus().symbolicFactorizationDone())
+          *fos << "Symbolic factorization should have been done during hierarchy setup, "
+              "but actually is missing. Anyway ... just do it right now." << std::endl;
+        if (not coarseSolver->getStatus().numericFactorizationDone())
+          *fos << "Numeric factorization should have been done during hierarchy setup, "
+              "but actually is missing. Anyway ... just do it right now." << std::endl;
+        coarseSolver->solve(tX.ptr(), tB.ptr());
 #else
-      *fos << "+++++++++++++++++++++++++++ WARNING +++++++++++++++++++++++++\n"
-           << "+ Coarse level direct solver requires Tpetra and Amesos2.   +\n"
-           << "+ Skipping the coarse level solve.                          +\n"
-           << "+++++++++++++++++++++++++++ WARNING +++++++++++++++++++++++++"
-           << std::endl;
+        *fos << "+++++++++++++++++++++++++++ WARNING +++++++++++++++++++++++++\n"
+             << "+ Coarse level direct solver requires Tpetra and Amesos2.   +\n"
+             << "+ Skipping the coarse level solve.                          +\n"
+             << "+++++++++++++++++++++++++++ WARNING +++++++++++++++++++++++++"
+             << std::endl;
 #endif
+      }
+      else // use AMG as coarse level solver
+      {
+
+        // Extract the hierarchy from the coarseSolverData
+        RCP<Hierarchy> amgHierarchy = coarseSolverData->get<RCP<Hierarchy>>("amg hierarchy object");
+
+        // Run a single V-cycle
+        amgHierarchy->Iterate(*compRhs, *compX, 1);
+      }
+
+      // Transform back to region format
+      Array<RCP<Vector> > quasiRegX(maxRegPerProc);
+      compositeToRegional(compX, quasiRegX, fineRegX,
+                          maxRegPerProc,
+                          quasiRegRowMaps[l],
+                          regRowMaps[l],
+                          regRowImporters[l]);
+
+      tm = Teuchos::null;
     }
-    else // use AMG as coarse level solver
-    {
-
-      // Extract the hierarchy from the coarseSolverData
-      RCP<Hierarchy> amgHierarchy = coarseSolverData->get<RCP<Hierarchy>>("amg hierarchy object");
-
-      // Run a single V-cycle
-      amgHierarchy->Iterate(*compRhs, *compX, 1);
-    }
-
-    // Transform back to region format
-    Array<RCP<Vector> > quasiRegX(maxRegPerProc);
-    compositeToRegional(compX, quasiRegX, fineRegX,
-                        maxRegPerProc,
-                        quasiRegRowMaps[l],
-                        regRowMaps[l],
-                        regRowImporters[l]);
-
-    tm = Teuchos::null;
-  }
+  } 
 
   return;
 } // vCycle

packages/muelu/research/mmayr/composite_to_regions/src/SetupRegionMatrix_def.hpp

@@ -87,7 +87,7 @@ Teuchos::Array<int> findCommonRegions(const GlobalOrdinal nodeA, ///< GID of fir
   Array<int> regionsA, regionsB;
   {
     LO nodeALID = nodesToRegionsMap->getLocalElement(nodeA);
-    LO nodeBLID = nodesToRegionsMap->getLocalElement(nodeA);
+    LO nodeBLID = nodesToRegionsMap->getLocalElement(nodeB);
     for (int i = 0; i < nodesToRegions.size(); ++i) {
       regionsA.push_back(nodesToRegions[i][nodeALID]);
       regionsB.push_back(nodesToRegions[i][nodeBLID]);

packages/muelu/research/regionMG/examples/structured/Driver_Structured_Regions.cpp

@@ -1286,7 +1286,7 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
   // std::cout << "p=" << myRank << " | compositeToRegionLIDs: " << compositeToRegionLIDs() << std::endl;
   // std::cout << "p=" << myRank << " | quasiRegionGIDs: " << quasiRegionGIDs << std::endl;
   // std::cout << "p=" << myRank << " | interfaceLIDs: " << interfaceLIDs() << std::endl;
-  std::cout << "p=" << myRank << " | quasiRegionCoordGIDs: " << quasiRegionCoordGIDs() << std::endl;
+  // std::cout << "p=" << myRank << " | quasiRegionCoordGIDs: " << quasiRegionCoordGIDs() << std::endl;
 
   // In our very particular case we know that a node is at most shared by 4 regions.
   // Other geometries will certainly have different constrains and a parallel reduction using MAX
@@ -1534,6 +1534,41 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
     Array<RCP<Vector> > regB(maxRegPerProc);
     compositeToRegional(B, quasiRegB, regB, maxRegPerProc, rowMapPerGrp,
                         revisedRowMapPerGrp, rowImportPerGrp);
+#ifdef DUMP_LOCALX_AND_A
+    FILE *fp;
+    char str[80];
+    sprintf(str,"theMatrix.%d",myRank);
+    fp = fopen(str,"w");
+    fprintf(fp, "%%%%MatrixMarket matrix coordinate real general\n");
+    LO numNzs = 0;
+    for (size_t kkk = 0; kkk < regionGrpMats[0]->getNodeNumRows(); kkk++) {
+      ArrayView<const LO> AAcols;
+      ArrayView<const SC> AAvals;
+      regionGrpMats[0]->getLocalRowView(kkk, AAcols, AAvals);
+      const int *Acols    = AAcols.getRawPtr();
+      const SC  *Avals = AAvals.getRawPtr();
+      numNzs += AAvals.size();
+    }
+    fprintf(fp, "%d %d %d\n",regionGrpMats[0]->getNodeNumRows(),regionGrpMats[0]->getNodeNumRows(),numNzs);
+
+    for (size_t kkk = 0; kkk < regionGrpMats[0]->getNodeNumRows(); kkk++) {
+      ArrayView<const LO> AAcols;
+      ArrayView<const SC> AAvals;
+      regionGrpMats[0]->getLocalRowView(kkk, AAcols, AAvals);
+      const int *Acols    = AAcols.getRawPtr();
+      const SC  *Avals = AAvals.getRawPtr();
+      LO RowLeng = AAvals.size();
+      for (LO kk = 0; kk < RowLeng; kk++) {
+          fprintf(fp, "%d %d %22.16e\n",kkk+1,Acols[kk]+1,Avals[kk]);
+      }
+    }
+    fclose(fp);
+    sprintf(str,"theX.%d",myRank);
+    fp = fopen(str,"w");
+    ArrayRCP<SC> lX= regX[0]->getDataNonConst(0);
+    for (size_t kkk = 0; kkk < regionGrpMats[0]->getNodeNumRows(); kkk++) fprintf(fp, "%22.16e\n",lX[kkk]);
+    fclose(fp);
+#endif
 
     //    printRegionalObject<Vector>("regB 0", regB, myRank, *fos);
 

packages/muelu/src/Graph/UncoupledAggregation/MueLu_AggregationPhase2bAlgorithm_kokkos_def.hpp

@@ -141,7 +141,7 @@ namespace MueLu {
                                     // (aggStat[neigh] == AGGREGATED)
                                     if (graph.isLocalNeighborVertex(neigh) &&
                                         aggStat(neigh) == AGGREGATED)
-                                      Kokkos::atomic_add(&aggWeight(vertex2AggId(neigh, 0), 0),
+                                      Kokkos::atomic_add(&aggWeight(vertex2AggId(neigh, 0)),
                                                          connectWeight(neigh));
                                   }
 
@@ -169,7 +169,7 @@ namespace MueLu {
                                     }
                                   }
                                   if (bestScore >= 0) {
-                                    aggStat(i, 0)      = AGGREGATED;
+                                    aggStat(i)         = AGGREGATED;
                                     vertex2AggId(i, 0) = bestAggId;
                                     procWinner(i, 0)   = myRank;
 
@@ -244,7 +244,7 @@ namespace MueLu {
               // (aggStat[neigh] == AGGREGATED)
               if (graph.isLocalNeighborVertex(neigh) &&
                   aggStat(neigh) == AGGREGATED)
-              Kokkos::atomic_add(&aggWeight(vertex2AggId(neigh, 0), 0),
+              Kokkos::atomic_add(&aggWeight(vertex2AggId(neigh, 0)),
                   connectWeight(neigh));
             }
           });

packages/muelu/src/Interface/MueLu_HierarchyManager.hpp

@@ -353,7 +353,7 @@ namespace MueLu {
             if (!M.is_null()) {
               Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(fileName,* M);
             }
-          }	  
+          }
 	  else if (L->IsAvailable(name)) {
 	    // Try nofactory
             RCP<T> M = L->template Get< RCP<T> >(name);

packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp

@@ -421,6 +421,7 @@ namespace MueLu {
     // FIXME: should it be here, or higher up
     RCP<FactoryManager> defaultManager = rcp(new FactoryManager());
     defaultManager->SetVerbLevel(this->verbosity_);
+    defaultManager->SetKokkosRefactor(useKokkos_);
 
     // We will ignore keeps0
     std::vector<keep_pair> keeps0;

packages/muelu/src/MueCentral/MueLu_FactoryManager_decl.hpp

@@ -176,6 +176,8 @@ namespace MueLu {
       useKokkos_ = useKokkos;
     }
 
+    bool GetKokkosRefactor() const { return useKokkos_; }
+
     //@}
 
     void Clean() const { defaultFactoryTable_.clear(); }

packages/muelu/src/Smoothers/MueLu_Amesos2Smoother_def.hpp

@@ -256,7 +256,7 @@ namespace MueLu {
       // add A
       for (size_t i = 0; i < N; i++) {
         for (size_t jj = rowPointers[i]; jj < rowPointers[i+1]; jj++) {
-          LO j = A->getColMap()->getGlobalElement(colIndices[jj]);
+          LO j = colMap->getLocalElement(A->getColMap()->getGlobalElement(colIndices[jj]));
           SC v = values[jj];
           newValues[i*M+j] += v;
         }

packages/muelu/test/interface/ParameterListInterpreter.cpp

@@ -78,8 +78,8 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
   using Teuchos::rcp;
   using Teuchos::TimeMonitor;
 
-  typedef typename Teuchos::ScalarTraits<SC>::magnitudeType real_type;
-  typedef Xpetra::MultiVector<real_type,LO,GO,NO> RealValuedMultiVector;
+  using real_type             = typename Teuchos::ScalarTraits<SC>::coordinateType;
+  using RealValuedMultiVector = Xpetra::MultiVector<real_type,LO,GO,NO>;
 
   // =========================================================================
   // MPI initialization using Teuchos
@@ -119,7 +119,7 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar
   matrixParameters.set("nx",         Teuchos::as<GO>(9999));
   matrixParameters.set("matrixType", "Laplace1D");
   RCP<Matrix>      A           = MueLuTests::TestHelpers::TestFactory<SC, LO, GO, NO>::Build1DPoisson(matrixParameters.get<GO>("nx"), lib);
-  RCP<RealValuedMultiVector> coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates<double,LO,GO,Map,RealValuedMultiVector>("1D", A->getRowMap(), matrixParameters);
+  RCP<RealValuedMultiVector> coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates<real_type,LO,GO,Map,RealValuedMultiVector>("1D", A->getRowMap(), matrixParameters);
 
   std::string prefix;
   if (useKokkos) {

packages/muelu/test/unit_tests_kokkos/CMakeLists.txt

@@ -26,6 +26,7 @@ APPEND_SET(SOURCES
   StructuredAggregation_kokkos.cpp
   #SaPFactory_kokkos.cpp
   TentativePFactory_kokkos.cpp
+  UseKokkos_kokkos.cpp
 )
 
 ### Tests that require other Trilinos packages