MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_ClassicalDropping.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_CLASSICALDROPPING_HPP
11 #define MUELU_CLASSICALDROPPING_HPP
12 
13 #include "MueLu_DroppingCommon.hpp"
14 #include "Kokkos_Core.hpp"
15 #include "Kokkos_ArithTraits.hpp"
16 #include "Xpetra_Matrix.hpp"
17 #include "MueLu_Utilities.hpp"
18 
19 namespace MueLu::ClassicalDropping {
20 
45 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, Misc::StrengthMeasure measure>
46 class DropFunctor {
47  public:
50  using local_matrix_type = typename matrix_type::local_matrix_type;
51  using scalar_type = typename local_matrix_type::value_type;
52  using local_ordinal_type = typename local_matrix_type::ordinal_type;
53  using memory_space = typename local_matrix_type::memory_space;
54  using diag_view_type = typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
55 
56  using results_view = Kokkos::View<DecisionType*, memory_space>;
57 
58  using ATS = Kokkos::ArithTraits<scalar_type>;
59  using magnitudeType = typename ATS::magnitudeType;
60  using mATS = Kokkos::ArithTraits<magnitudeType>;
61  using boundary_nodes_view = Kokkos::View<const bool*, memory_space>;
62 
63  private:
66  diag_view_type diag; // corresponds to overlapped diagonal
69 
70  public:
72  : A(A_.getLocalMatrixDevice())
73  , eps(threshold)
74  , results(results_) {
75  // Construct ghosted matrix diagonal
76  if constexpr ((measure == Misc::SmoothedAggregationMeasure) || (measure == Misc::SmoothedAggregationMeasure)) {
78  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
79  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
80  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
82  auto lcl2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
83  diag = Kokkos::subview(lcl2d, Kokkos::ALL(), 0);
84  }
85  }
86 
87  KOKKOS_FORCEINLINE_FUNCTION
88  void operator()(const local_ordinal_type rlid) const {
89  auto row = A.rowConst(rlid);
90  size_t offset = A.graph.row_map(rlid);
91  for (local_ordinal_type k = 0; k < row.length; ++k) {
92  auto clid = row.colidx(k);
93 
94  auto val = row.value(k);
95 
96  if constexpr (measure == Misc::SmoothedAggregationMeasure) {
97  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
98  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
99 
100  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
101  results(offset + k));
102  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
103  auto neg_aij = -ATS::real(val);
104  auto max_neg_aik = eps * ATS::real(diag(rlid));
105  results(offset + k) = Kokkos::max((neg_aij <= max_neg_aik) ? DROP : KEEP,
106  results(offset + k));
107  } else if constexpr (measure == Misc::SignedSmoothedAggregationMeasure) {
108  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
109  const bool is_nonpositive = ATS::real(val) <= mATS::zero();
110  magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
111  // + |a_ij|^2, if a_ij < 0, - |a_ij|^2 if a_ij >=0
112  if (is_nonpositive)
113  aij2 = -aij2;
114  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
115  results(offset + k));
116  }
117  }
118  }
119 };
120 
121 // helper function to allow partial template deduction
122 template <Misc::StrengthMeasure measure, class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
126  auto functor = DropFunctor<Scalar, LocalOrdinal, GlobalOrdinal, Node, measure>(A_, threshold, results_);
127  return functor;
128 }
129 
130 } // namespace MueLu::ClassicalDropping
131 
132 #endif
static Teuchos::RCP< Vector > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Return vector containing: max_{i=k}(-a_ik), for each for i in the matrix.
Kokkos::View< const bool *, memory_space > boundary_nodes_view
Kokkos::ArithTraits< magnitudeType > mATS
typename local_matrix_type::value_type scalar_type
typename local_matrix_type::memory_space memory_space
DropFunctor(matrix_type &A_, magnitudeType threshold, results_view &results_)
typename matrix_type::local_matrix_type local_matrix_type
Kokkos::View< DecisionType *, memory_space > results_view
typename Kokkos::DualView< const scalar_type *, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged >::t_dev diag_view_type
KOKKOS_FORCEINLINE_FUNCTION void operator()(const local_ordinal_type rlid) const
auto make_drop_functor(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A_, typename DropFunctor< Scalar, LocalOrdinal, GlobalOrdinal, Node, measure >::magnitudeType threshold, typename DropFunctor< Scalar, LocalOrdinal, GlobalOrdinal, Node, measure >::results_view &results_)
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
Kokkos::ArithTraits< scalar_type > ATS
typename local_matrix_type::ordinal_type local_ordinal_type