diff --git a/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_kokkos_def.hpp b/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_kokkos_def.hpp index b28d2e832765..301e0dbf1c7b 100644 --- a/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_kokkos_def.hpp +++ b/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_kokkos_def.hpp @@ -217,6 +217,8 @@ std::tuple("filtered matrix: use root stencil"); const bool useSpreadLumping = pL.get("filtered matrix: use spread lumping"); + + const Scalar filteringDirichletThreshold = as(pL.get("filtered matrix: Dirichlet threshold")); TEUCHOS_ASSERT(!useRootStencil); TEUCHOS_ASSERT(!useSpreadLumping); @@ -692,18 +694,18 @@ std::tuple(lclA, results, lclFilteredA, lclGraph); + auto fillFunctor = MatrixConstruction::PointwiseFillReuseFunctor(lclA, results, lclFilteredA, lclGraph, filteringDirichletThreshold); Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_reuse", range, fillFunctor); } else { - auto fillFunctor = MatrixConstruction::PointwiseFillNoReuseFunctor(lclA, results, lclFilteredA); + auto fillFunctor = MatrixConstruction::PointwiseFillNoReuseFunctor(lclA, results, lclFilteredA, filteringDirichletThreshold); Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_noreuse", range, fillFunctor); } } else { if (reuseGraph) { - auto fillFunctor = MatrixConstruction::PointwiseFillReuseFunctor(lclA, results, lclFilteredA, lclGraph); + auto fillFunctor = MatrixConstruction::PointwiseFillReuseFunctor(lclA, results, lclFilteredA, lclGraph, filteringDirichletThreshold); Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_reuse", range, fillFunctor); } else { - auto fillFunctor = MatrixConstruction::PointwiseFillNoReuseFunctor(lclA, results, lclFilteredA); + auto fillFunctor = MatrixConstruction::PointwiseFillNoReuseFunctor(lclA, results, lclFilteredA, filteringDirichletThreshold); Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_noreuse", range, fillFunctor); } } @@ -854,6 +856,9 @@ std::tuple("filtered matrix: use root stencil"); const bool useSpreadLumping = pL.get("filtered matrix: use spread lumping"); + + const Scalar filteringDirichletThreshold = as(pL.get("filtered matrix: Dirichlet threshold")); + TEUCHOS_ASSERT(!useRootStencil); TEUCHOS_ASSERT(!useSpreadLumping); @@ -1095,18 +1100,18 @@ std::tuple(lclA, blkPartSize, colTranslation, results, lclFilteredA, lclGraph); + auto fillFunctor = MatrixConstruction::VectorFillFunctor(lclA, blkPartSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold); Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_reuse", range, fillFunctor); } else { - auto fillFunctor = MatrixConstruction::VectorFillFunctor(lclA, blkPartSize, colTranslation, results, lclFilteredA, lclGraph); + auto fillFunctor = MatrixConstruction::VectorFillFunctor(lclA, blkPartSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold); Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_lumped_noreuse", range, fillFunctor); } } else { if (reuseGraph) { - auto fillFunctor = MatrixConstruction::VectorFillFunctor(lclA, blkSize, colTranslation, results, lclFilteredA, lclGraph); + auto fillFunctor = MatrixConstruction::VectorFillFunctor(lclA, blkSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold); Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_reuse", range, fillFunctor); } else { - auto fillFunctor = MatrixConstruction::VectorFillFunctor(lclA, blkSize, colTranslation, results, lclFilteredA, lclGraph); + auto fillFunctor = MatrixConstruction::VectorFillFunctor(lclA, blkSize, colTranslation, results, lclFilteredA, lclGraph, filteringDirichletThreshold); Kokkos::parallel_for("MueLu::CoalesceDrop::Fill_unlumped_noreuse", range, fillFunctor); } } diff --git a/packages/muelu/src/Graph/MatrixTransformation/MueLu_MatrixConstruction.hpp b/packages/muelu/src/Graph/MatrixTransformation/MueLu_MatrixConstruction.hpp index 1a5f2729c72e..640e1a6973c4 100644 --- a/packages/muelu/src/Graph/MatrixTransformation/MueLu_MatrixConstruction.hpp +++ b/packages/muelu/src/Graph/MatrixTransformation/MueLu_MatrixConstruction.hpp @@ -255,14 +255,17 @@ class PointwiseFillReuseFunctor { results_view results; local_matrix_type filteredA; local_graph_type graph; + scalar_type dirichletThreshold; const scalar_type zero = ATS::zero(); + const scalar_type one = ATS::one(); public: - PointwiseFillReuseFunctor(local_matrix_type& A_, results_view& results_, local_matrix_type& filteredA_, local_graph_type& graph_) + PointwiseFillReuseFunctor(local_matrix_type& A_, results_view& results_, local_matrix_type& filteredA_, local_graph_type& graph_, scalar_type dirichletThreshold_) : A(A_) , results(results_) , filteredA(filteredA_) - , graph(graph_) {} + , graph(graph_) + , dirichletThreshold(dirichletThreshold_) {} KOKKOS_INLINE_FUNCTION void operator()(const local_ordinal_type rlid) const { @@ -300,6 +303,8 @@ class PointwiseFillReuseFunctor { } if constexpr (lumping) { rowFilteredA.value(diagOffset) += diagCorrection; + if ((dirichletThreshold >= 0.0) && (ATS::real(rowFilteredA.value(diagOffset)) <= dirichletThreshold)) + rowFilteredA.value(diagOffset) = one; } } }; @@ -323,13 +328,16 @@ class PointwiseFillNoReuseFunctor { local_matrix_type A; results_view results; local_matrix_type filteredA; + scalar_type dirichletThreshold; const scalar_type zero = ATS::zero(); + const scalar_type one = ATS::one(); public: - PointwiseFillNoReuseFunctor(local_matrix_type& A_, results_view& results_, local_matrix_type& filteredA_) + PointwiseFillNoReuseFunctor(local_matrix_type& A_, results_view& results_, local_matrix_type& filteredA_, scalar_type dirichletThreshold_) : A(A_) , results(results_) - , filteredA(filteredA_) {} + , filteredA(filteredA_) + , dirichletThreshold(dirichletThreshold_) {} KOKKOS_INLINE_FUNCTION void operator()(const local_ordinal_type rlid) const { @@ -356,6 +364,8 @@ class PointwiseFillNoReuseFunctor { } if constexpr (lumping) { rowFilteredA.value(diagOffset) += diagCorrection; + if ((dirichletThreshold >= 0.0) && (ATS::real(rowFilteredA.value(diagOffset)) <= dirichletThreshold)) + rowFilteredA.value(diagOffset) = one; } } }; @@ -727,16 +737,19 @@ class VectorFillFunctor { results_view results; local_matrix_type filteredA; local_graph_type graph; + scalar_type dirichletThreshold; const scalar_type zero = ATS::zero(); + const scalar_type one = ATS::one(); public: - VectorFillFunctor(local_matrix_type& A_, local_ordinal_type blockSize_, block_indices_view_type ghosted_point_to_block_, results_view& results_, local_matrix_type& filteredA_, local_graph_type& graph_) + VectorFillFunctor(local_matrix_type& A_, local_ordinal_type blockSize_, block_indices_view_type ghosted_point_to_block_, results_view& results_, local_matrix_type& filteredA_, local_graph_type& graph_, scalar_type dirichletThreshold_) : A(A_) , blockSize(blockSize_) , ghosted_point_to_block(ghosted_point_to_block_) , results(results_) , filteredA(filteredA_) - , graph(graph_) {} + , graph(graph_) + , dirichletThreshold(dirichletThreshold_) {} KOKKOS_INLINE_FUNCTION void operator()(const local_ordinal_type brlid) const { @@ -773,6 +786,8 @@ class VectorFillFunctor { } if constexpr (lumping) { rowFilteredA.value(diagOffset) += diagCorrection; + if ((dirichletThreshold >= 0.0) && (ATS::real(rowFilteredA.value(diagOffset)) <= dirichletThreshold)) + rowFilteredA.value(diagOffset) = one; } } diff --git a/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp b/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp index 36146faa491a..643b92518838 100644 --- a/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp +++ b/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp @@ -1059,6 +1059,11 @@ void ParameterListInterpreter:: MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use lumping", bool, dropParams); MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse graph", bool, dropParams); MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse eigenvalue", bool, dropParams); + MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use root stencil", bool, dropParams); + MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: Dirichlet threshold", double, dropParams); + MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use spread lumping", bool, dropParams); + MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom growth factor", double, dropParams); + MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom cap", double, dropParams); } if (!amalgFact.is_null()) @@ -1313,7 +1318,6 @@ void ParameterListInterpreter:: // Matrix analysis if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "matrix: compute analysis", bool, true)) { - RCP matrixAnalysisFact = rcp(new MatrixAnalysisFactory()); if (!RAP.is_null())