MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_DistanceLaplacianDropping.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_DISTANCELAPLACIANDROPPING_HPP
11 #define MUELU_DISTANCELAPLACIANDROPPING_HPP
12 
13 #include "Kokkos_Macros.hpp"
14 #include "KokkosBatched_LU_Decl.hpp"
15 #include "KokkosBatched_LU_Serial_Impl.hpp"
16 #include "KokkosBatched_Trsv_Decl.hpp"
17 #include "KokkosBatched_Trsv_Serial_Impl.hpp"
18 #include "MueLu_DroppingCommon.hpp"
19 #include "Kokkos_Core.hpp"
20 #include "Kokkos_ArithTraits.hpp"
21 #include "Teuchos_RCP.hpp"
22 #include "Xpetra_Matrix.hpp"
23 #include "Xpetra_MultiVector.hpp"
24 #include "Xpetra_MultiVectorFactory.hpp"
25 
26 namespace MueLu::DistanceLaplacian {
27 
32 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
34  private:
36  using local_matrix_type = typename matrix_type::local_matrix_type;
37  using scalar_type = typename local_matrix_type::value_type;
39  using ATS = Kokkos::ArithTraits<scalar_type>;
40  using impl_scalar_type = typename ATS::val_type;
41  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
42  using magnitudeType = typename implATS::magnitudeType;
43  using magATS = Kokkos::ArithTraits<magnitudeType>;
45  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
46 
49 
52 
53  public:
55  coordsMV = coords_;
56  auto importer = A.getCrsGraph()->getImporter();
57  if (!importer.is_null()) {
59  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
60  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
61  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
62  } else {
63  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
65  }
66  }
67 
68  KOKKOS_FORCEINLINE_FUNCTION
70  magnitudeType d = magATS::zero();
71  magnitudeType s;
72  for (size_t j = 0; j < coords.extent(1); ++j) {
73  s = coords(row, j) - ghostedCoords(col, j);
74  d += s * s;
75  }
76  return d;
77  }
78 };
79 
80 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
82  private:
84  using local_matrix_type = typename matrix_type::local_matrix_type;
85  using scalar_type = typename local_matrix_type::value_type;
87  using ATS = Kokkos::ArithTraits<scalar_type>;
88  using impl_scalar_type = typename ATS::val_type;
89  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
90  using magnitudeType = typename implATS::magnitudeType;
91  using magATS = Kokkos::ArithTraits<magnitudeType>;
93  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
95  using local_material_type = typename material_type::dual_view_type_const::t_dev;
96 
99 
102 
105 
108 
109  public:
111  coordsMV = coords_;
112  materialMV = material_;
113  auto importer = A.getCrsGraph()->getImporter();
114  if (!importer.is_null()) {
116  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
117  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
118  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
119 
121  ghostedMaterialMV->doImport(*materialMV, *importer, Xpetra::INSERT);
122  material = materialMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
123  ghostedMaterial = ghostedMaterialMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
124  } else {
125  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
127 
128  material = materialMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
130  }
131  }
132 
133  KOKKOS_INLINE_FUNCTION
135  // || x_row - x_col ||_S^2
136  // where
137  // S = 1/material(row) * Identity
138  magnitudeType d = magATS::zero();
139  magnitudeType d_row = magATS::zero();
140  magnitudeType d_col = magATS::zero();
141  magnitudeType s;
142 
143  for (size_t j = 0; j < coords.extent(1); ++j) {
144  s = coords(row, j) - ghostedCoords(col, j);
145  d += s * s;
146  }
147 
148  d_row = d / implATS::magnitude(ghostedMaterial(row, 0));
149  d_col = d / implATS::magnitude(ghostedMaterial(col, 0));
150 
151  return Kokkos::max(d_row, d_col);
152  }
153 };
154 
155 template <class local_ordinal_type, class material_vector_type, class material_matrix_type>
157  private:
158  material_vector_type material_vector;
159  material_matrix_type material_matrix;
160 
161  public:
162  TensorInversion(material_vector_type& material_vector_, material_matrix_type& material_matrix_)
163  : material_vector(material_vector_)
164  , material_matrix(material_matrix_) {}
165 
166  KOKKOS_FORCEINLINE_FUNCTION
167  void operator()(local_ordinal_type i) const {
168  for (size_t j = 0; j < material_matrix.extent(1); ++j) {
169  for (size_t k = 0; k < material_matrix.extent(2); ++k) {
170  material_matrix(i, j, k) = material_vector(i, j * material_matrix.extent(1) + k);
171  }
172  }
173  auto matrix_material = Kokkos::subview(material_matrix, i, Kokkos::ALL(), Kokkos::ALL());
174  KokkosBatched::SerialLU<KokkosBatched::Algo::LU::Unblocked>::invoke(matrix_material);
175  }
176 };
177 
178 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
180  private:
182  using local_matrix_type = typename matrix_type::local_matrix_type;
183  using scalar_type = typename local_matrix_type::value_type;
185  using ATS = Kokkos::ArithTraits<scalar_type>;
186  using impl_scalar_type = typename ATS::val_type;
187  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
188  using magnitudeType = typename implATS::magnitudeType;
189  using magATS = Kokkos::ArithTraits<magnitudeType>;
191  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
193  using memory_space = typename local_matrix_type::memory_space;
194 
195  using local_material_type = Kokkos::View<impl_scalar_type***, memory_space>;
196  using local_dist_type = Kokkos::View<impl_scalar_type**, memory_space>;
197 
200 
203 
205 
207 
208  const scalar_type one = ATS::one();
209 
210  public:
212  coordsMV = coords_;
213 
214  auto importer = A.getCrsGraph()->getImporter();
215  if (!importer.is_null()) {
217  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
218  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
219  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
220  } else {
221  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
223  }
224 
225  {
226  Teuchos::RCP<material_type> ghostedMaterial;
227  if (!importer.is_null()) {
228  ghostedMaterial = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(importer->getTargetMap(), material_->getNumVectors());
229  ghostedMaterial->doImport(*material_, *importer, Xpetra::INSERT);
230  } else {
231  ghostedMaterial = material_;
232  }
233 
234  using execution_space = typename Node::execution_space;
235  using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
236 
237  local_ordinal_type dim = std::sqrt(material_->getNumVectors());
238  auto lclMaterial = ghostedMaterial->getLocalViewDevice(Xpetra::Access::ReadOnly);
239  material = local_material_type("material", lclMaterial.extent(0), dim, dim);
240  lcl_dist = local_dist_type("material", lclMaterial.extent(0), dim);
242  Kokkos::parallel_for("MueLu:TensorMaterialDistanceFunctor::inversion", range_type(0, lclMaterial.extent(0)), functor);
243  }
244  }
245 
246  KOKKOS_INLINE_FUNCTION
248  // || x_row - x_col ||_S^2
249  // where
250  // S = inv(material(col))
251 
252  // row material
253  impl_scalar_type d_row = implATS::zero();
254  {
255  auto matrix_row_material = Kokkos::subview(material, row, Kokkos::ALL(), Kokkos::ALL());
256  auto dist = Kokkos::subview(lcl_dist, row, Kokkos::ALL());
257 
258  for (size_t j = 0; j < coords.extent(1); ++j) {
259  dist(j) = coords(row, j) - ghostedCoords(col, j);
260  }
261 
262  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Lower, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::Unit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_row_material, dist);
263  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Upper, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::NonUnit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_row_material, dist);
264 
265  for (size_t j = 0; j < coords.extent(1); ++j) {
266  d_row += dist(j) * (coords(row, j) - ghostedCoords(col, j));
267  }
268  }
269 
270  // column material
271  impl_scalar_type d_col = implATS::zero();
272  {
273  auto matrix_col_material = Kokkos::subview(material, col, Kokkos::ALL(), Kokkos::ALL());
274  auto dist = Kokkos::subview(lcl_dist, row, Kokkos::ALL());
275 
276  for (size_t j = 0; j < coords.extent(1); ++j) {
277  dist(j) = coords(row, j) - ghostedCoords(col, j);
278  }
279 
280  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Lower, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::Unit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_col_material, dist);
281  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Upper, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::NonUnit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_col_material, dist);
282 
283  for (size_t j = 0; j < coords.extent(1); ++j) {
284  d_col += dist(j) * (coords(row, j) - ghostedCoords(col, j));
285  }
286  }
287 
288  return Kokkos::max(implATS::magnitude(d_row), implATS::magnitude(d_col));
289  }
290 };
291 
295 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
298  DistanceFunctorType& distFunctor) {
299  using scalar_type = Scalar;
300  using local_ordinal_type = LocalOrdinal;
301  using global_ordinal_type = GlobalOrdinal;
302  using node_type = Node;
303  using ATS = Kokkos::ArithTraits<scalar_type>;
304  using impl_scalar_type = typename ATS::val_type;
305  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
306  using magnitudeType = typename implATS::magnitudeType;
307  using execution_space = typename Node::execution_space;
308  using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
309 
311  {
312  auto lclA = A.getLocalMatrixDevice();
313  auto lclDiag = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
314 
315  Kokkos::parallel_for(
316  "MueLu:CoalesceDropF:Build:scalar_filter:laplacian_diag",
317  range_type(0, lclA.numRows()),
318  KOKKOS_LAMBDA(const local_ordinal_type& row) {
319  auto rowView = lclA.rowConst(row);
320  auto length = rowView.length;
321 
322  magnitudeType d;
323  impl_scalar_type d2 = implATS::zero();
324  bool haveAddedToDiag = false;
325  for (local_ordinal_type colID = 0; colID < length; colID++) {
326  auto col = rowView.colidx(colID);
327  if (row != col) {
328  d = distFunctor.distance2(row, col);
329  d2 += implATS::one() / d;
330  haveAddedToDiag = true;
331  }
332  }
333 
334  // Deal with the situation where boundary conditions have only been enforced on rows, but not on columns.
335  // We enforce dropping of these entries by assigning a very large number to the diagonal entries corresponding to BCs.
336  lclDiag(row, 0) = !haveAddedToDiag ? implATS::squareroot(implATS::rmax()) : d2;
337  });
338  }
339  auto importer = A.getCrsGraph()->getImporter();
340  if (!importer.is_null()) {
342  ghostedDiag->doImport(*diag, *importer, Xpetra::INSERT);
343  return ghostedDiag;
344  } else {
345  return diag;
346  }
347 }
348 
359 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
360 class DropFunctor {
361  private:
363  using local_matrix_type = typename matrix_type::local_matrix_type;
364  using scalar_type = typename local_matrix_type::value_type;
365  using local_ordinal_type = typename local_matrix_type::ordinal_type;
366  using memory_space = typename local_matrix_type::memory_space;
368  using diag_view_type = typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
369 
370  using results_view = Kokkos::View<DecisionType*, memory_space>;
371 
372  using ATS = Kokkos::ArithTraits<scalar_type>;
373  using magnitudeType = typename ATS::magnitudeType;
374  using boundary_nodes_view = Kokkos::View<const bool*, memory_space>;
375 
379  diag_view_type diag; // corresponds to overlapped diagonal
380  DistanceFunctorType dist2;
382  const scalar_type one = ATS::one();
383 
384  public:
385  DropFunctor(matrix_type& A_, magnitudeType threshold, DistanceFunctorType& dist2_, results_view& results_)
386  : A(A_.getLocalMatrixDevice())
387  , eps(threshold)
388  , dist2(dist2_)
389  , results(results_) {
390  diagVec = getDiagonal(A_, dist2);
391  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
392  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
393  }
394 
395  KOKKOS_FORCEINLINE_FUNCTION
396  void operator()(local_ordinal_type rlid) const {
397  auto row = A.rowConst(rlid);
398  const size_t offset = A.graph.row_map(rlid);
399  for (local_ordinal_type k = 0; k < row.length; ++k) {
400  auto clid = row.colidx(k);
401 
402  scalar_type val;
403  if (rlid != clid) {
404  val = one / dist2.distance2(rlid, clid);
405  } else {
406  val = diag(rlid);
407  }
408  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
409  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
410 
411  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
412  results(offset + k));
413  }
414  }
415 };
416 
417 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
419  private:
421  using local_matrix_type = typename matrix_type::local_matrix_type;
422  using scalar_type = typename local_matrix_type::value_type;
423  using local_ordinal_type = typename local_matrix_type::ordinal_type;
424  using memory_space = typename local_matrix_type::memory_space;
425  using block_indices_view_type = Kokkos::View<local_ordinal_type*, memory_space>;
427  using diag_view_type = typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
428 
429  using results_view = Kokkos::View<DecisionType*, memory_space>;
430 
431  using ATS = Kokkos::ArithTraits<scalar_type>;
432  using magnitudeType = typename ATS::magnitudeType;
433  using boundary_nodes_view = Kokkos::View<const bool*, memory_space>;
434 
438  diag_view_type diag; // corresponds to overlapped diagonal
439  DistanceFunctorType dist2;
443  const scalar_type one = ATS::one();
444 
445  public:
446  VectorDropFunctor(matrix_type& A_, matrix_type& mergedA_, magnitudeType threshold, DistanceFunctorType& dist2_, results_view& results_, block_indices_view_type point_to_block_, block_indices_view_type ghosted_point_to_block_)
447  : A(A_.getLocalMatrixDevice())
448  , eps(threshold)
449  , dist2(dist2_)
450  , results(results_)
451  , point_to_block(point_to_block_)
452  , ghosted_point_to_block(ghosted_point_to_block_) {
453  diagVec = getDiagonal(mergedA_, dist2);
454  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
455  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
456  }
457 
458  KOKKOS_FORCEINLINE_FUNCTION
459  void operator()(local_ordinal_type rlid) const {
460  auto brlid = point_to_block(rlid);
461  auto row = A.rowConst(rlid);
462  const size_t offset = A.graph.row_map(rlid);
463  for (local_ordinal_type k = 0; k < row.length; ++k) {
464  auto clid = row.colidx(k);
465  auto bclid = ghosted_point_to_block(clid);
466 
467  scalar_type val;
468  if (brlid != bclid) {
469  val = one / dist2.distance2(brlid, bclid);
470  } else {
471  val = diag(brlid);
472  }
473  auto aiiajj = ATS::magnitude(diag(brlid)) * ATS::magnitude(diag(bclid)); // |a_ii|*|a_jj|
474  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
475 
476  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
477  results(offset + k));
478  }
479  }
480 };
481 
482 } // namespace MueLu::DistanceLaplacian
483 
484 #endif
Drops entries the unscaled distance Laplacian.
typename material_type::dual_view_type_const::t_dev local_material_type
MueLu::DefaultLocalOrdinal LocalOrdinal
typename matrix_type::local_matrix_type local_matrix_type
UnweightedDistanceFunctor(matrix_type &A, Teuchos::RCP< coords_type > &coords_)
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
typename coords_type::dual_view_type_const::t_dev local_coords_type
DropFunctor(matrix_type &A_, magnitudeType threshold, DistanceFunctorType &dist2_, results_view &results_)
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
VectorDropFunctor(matrix_type &A_, matrix_type &mergedA_, magnitudeType threshold, DistanceFunctorType &dist2_, results_view &results_, block_indices_view_type point_to_block_, block_indices_view_type ghosted_point_to_block_)
typename coords_type::dual_view_type_const::t_dev local_coords_type
typename matrix_type::local_matrix_type local_matrix_type
KOKKOS_INLINE_FUNCTION magnitudeType distance2(const local_ordinal_type row, const local_ordinal_type col) const
KOKKOS_INLINE_FUNCTION magnitudeType distance2(const local_ordinal_type row, const local_ordinal_type col) const
Kokkos::View< local_ordinal_type *, memory_space > block_indices_view_type
typename local_matrix_type::memory_space memory_space
MueLu::DefaultNode Node
virtual void doImport(const DistObject< Scalar, LocalOrdinal, GlobalOrdinal, Node > &source, const Import< LocalOrdinal, GlobalOrdinal, Node > &importer, CombineMode CM)=0
TensorMaterialDistanceFunctor(matrix_type &A, Teuchos::RCP< coords_type > &coords_, Teuchos::RCP< material_type > &material_)
typename local_matrix_type::ordinal_type local_ordinal_type
typename local_matrix_type::memory_space memory_space
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
TensorInversion(material_vector_type &material_vector_, material_matrix_type &material_matrix_)
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
Kokkos::View< DecisionType *, memory_space > results_view
typename coords_type::dual_view_type_const::t_dev local_coords_type
ScalarMaterialDistanceFunctor(matrix_type &A, Teuchos::RCP< coords_type > &coords_, Teuchos::RCP< material_type > &material_)
Kokkos::View< const bool *, memory_space > boundary_nodes_view
KOKKOS_FORCEINLINE_FUNCTION magnitudeType distance2(const local_ordinal_type row, const local_ordinal_type col) const
Kokkos::View< impl_scalar_type **, memory_space > local_dist_type
virtual RCP< const CrsGraph > getCrsGraph() const =0
typename local_matrix_type::ordinal_type local_ordinal_type
typename local_matrix_type::value_type scalar_type
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type i) const
Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getDiagonal(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, DistanceFunctorType &distFunctor)
Kokkos::View< impl_scalar_type ***, memory_space > local_material_type
Kokkos::View< const bool *, memory_space > boundary_nodes_view
typename Kokkos::DualView< const scalar_type *, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged >::t_dev diag_view_type
virtual size_t getNumVectors() const =0
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const