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 
92 #ifdef MUELU_COALESCE_DROP_DEBUG
93  {
94  Kokkos::printf("SoC: ");
95  for (local_ordinal_type k = 0; k < row.length; ++k) {
96  auto clid = row.colidx(k);
97 
98  auto val = row.value(k);
99 
100  if constexpr (measure == Misc::SmoothedAggregationMeasure) {
101  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
102  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
103 
104  Kokkos::printf("%5f ", ATS::sqrt(aij2 / aiiajj));
105  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
106  auto neg_aij = -ATS::real(val);
107  auto max_neg_aik = eps * ATS::real(diag(rlid));
108  results(offset + k) = Kokkos::max((neg_aij < max_neg_aik) ? DROP : KEEP,
109  results(offset + k));
110  Kokkos::printf("%5f ", neg_aij / max_neg_aik);
111  } else if constexpr (measure == Misc::SignedSmoothedAggregationMeasure) {
112  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
113  const bool is_nonpositive = ATS::real(val) <= mATS::zero();
114  magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
115  // + |a_ij|^2, if a_ij < 0, - |a_ij|^2 if a_ij >=0
116  if (!is_nonpositive)
117  aij2 = -aij2;
118  Kokkos::printf("%5f ", ATS::sqrt(aij2 / aiiajj));
119  }
120  }
121  Kokkos::printf("\n");
122  }
123 #endif
124 
125  for (local_ordinal_type k = 0; k < row.length; ++k) {
126  auto clid = row.colidx(k);
127 
128  auto val = row.value(k);
129 
130  if constexpr (measure == Misc::SmoothedAggregationMeasure) {
131  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
132  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
133 
134  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
135  results(offset + k));
136  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
137  auto neg_aij = -ATS::real(val);
138  auto max_neg_aik = eps * ATS::real(diag(rlid));
139  results(offset + k) = Kokkos::max((neg_aij < max_neg_aik) ? DROP : KEEP,
140  results(offset + k));
141  } else if constexpr (measure == Misc::SignedSmoothedAggregationMeasure) {
142  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
143  const bool is_nonpositive = ATS::real(val) <= mATS::zero();
144  magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
145  // + |a_ij|^2, if a_ij < 0, - |a_ij|^2 if a_ij >=0
146  if (!is_nonpositive)
147  aij2 = -aij2;
148  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
149  results(offset + k));
150  }
151  }
152  }
153 };
154 
155 // helper function to allow partial template deduction
156 template <Misc::StrengthMeasure measure, class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
160  auto functor = DropFunctor<Scalar, LocalOrdinal, GlobalOrdinal, Node, measure>(A_, threshold, results_);
161  return functor;
162 }
163 
164 } // namespace MueLu::ClassicalDropping
165 
166 #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