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 
84 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class weight_type>
86  private:
88  using local_matrix_type = typename matrix_type::local_matrix_type;
89  using scalar_type = typename local_matrix_type::value_type;
91  using ATS = Kokkos::ArithTraits<scalar_type>;
92  using impl_scalar_type = typename ATS::val_type;
93  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
94  using magnitudeType = typename implATS::magnitudeType;
95  using magATS = Kokkos::ArithTraits<magnitudeType>;
97  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
98 
101 
104 
105  weight_type weight;
106 
107  public:
108  WeightedDistanceFunctor(matrix_type& A, Teuchos::RCP<coords_type>& coords_, weight_type weight_) {
109  coordsMV = coords_;
110  auto importer = A.getCrsGraph()->getImporter();
111  if (!importer.is_null()) {
113  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
114  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
115  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
116  } else {
117  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
119  }
120  weight = weight_;
121  TEUCHOS_ASSERT(weight.extent(0) == coordsMV->getNumVectors());
122  }
123 
124  KOKKOS_FORCEINLINE_FUNCTION
126  magnitudeType d = magATS::zero();
127  magnitudeType w;
128  magnitudeType s;
129  for (size_t j = 0; j < coords.extent(1); ++j) {
130  s = coords(row, j) - ghostedCoords(col, j);
131  w = weight(j);
132  d += w * s * s;
133  }
134  return d;
135  }
136 };
137 
142 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class weight_type>
144  private:
146  using local_matrix_type = typename matrix_type::local_matrix_type;
147  using scalar_type = typename local_matrix_type::value_type;
149  using ATS = Kokkos::ArithTraits<scalar_type>;
150  using impl_scalar_type = typename ATS::val_type;
151  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
152  using magnitudeType = typename implATS::magnitudeType;
153  using magATS = Kokkos::ArithTraits<magnitudeType>;
155  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
156 
159 
162 
163  weight_type weight;
165 
166  public:
167  BlockWeightedDistanceFunctor(matrix_type& A, Teuchos::RCP<coords_type>& coords_, weight_type weight_, local_ordinal_type interleaved_blocksize_) {
168  coordsMV = coords_;
169  auto importer = A.getCrsGraph()->getImporter();
170  if (!importer.is_null()) {
172  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
173  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
174  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
175  } else {
176  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
178  }
179  weight = weight_;
180  interleaved_blocksize = interleaved_blocksize_;
181  TEUCHOS_ASSERT(weight.extent(0) == coordsMV->getNumVectors());
182  }
183 
184  KOKKOS_FORCEINLINE_FUNCTION
186  magnitudeType d = magATS::zero();
187  magnitudeType w;
188  magnitudeType s;
189  local_ordinal_type block_id = row % interleaved_blocksize;
190  local_ordinal_type block_start = block_id * interleaved_blocksize;
191  for (size_t j = 0; j < coords.extent(1); ++j) {
192  s = coords(row, j) - ghostedCoords(col, j);
193  w = weight(block_start + j);
194  d += w * s * s;
195  }
196  return d;
197  }
198 };
199 
200 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
202  private:
204  using local_matrix_type = typename matrix_type::local_matrix_type;
205  using scalar_type = typename local_matrix_type::value_type;
207  using ATS = Kokkos::ArithTraits<scalar_type>;
208  using impl_scalar_type = typename ATS::val_type;
209  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
210  using magnitudeType = typename implATS::magnitudeType;
211  using magATS = Kokkos::ArithTraits<magnitudeType>;
213  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
215  using local_material_type = typename material_type::dual_view_type_const::t_dev;
216 
219 
222 
225 
228 
229  public:
231  coordsMV = coords_;
232  materialMV = material_;
233  auto importer = A.getCrsGraph()->getImporter();
234  if (!importer.is_null()) {
236  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
237  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
238  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
239 
241  ghostedMaterialMV->doImport(*materialMV, *importer, Xpetra::INSERT);
242  material = materialMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
243  ghostedMaterial = ghostedMaterialMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
244  } else {
245  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
247 
248  material = materialMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
250  }
251  }
252 
253  KOKKOS_INLINE_FUNCTION
255  // || x_row - x_col ||_S^2
256  // where
257  // S = 1/material(row) * Identity
258  magnitudeType d = magATS::zero();
259  magnitudeType d_row = magATS::zero();
260  magnitudeType d_col = magATS::zero();
261  magnitudeType s;
262 
263  for (size_t j = 0; j < coords.extent(1); ++j) {
264  s = coords(row, j) - ghostedCoords(col, j);
265  d += s * s;
266  }
267 
268  d_row = d / implATS::magnitude(ghostedMaterial(row, 0));
269  d_col = d / implATS::magnitude(ghostedMaterial(col, 0));
270 
271  return Kokkos::max(d_row, d_col);
272  }
273 };
274 
275 template <class local_ordinal_type, class material_vector_type, class material_matrix_type>
277  private:
278  material_vector_type material_vector;
279  material_matrix_type material_matrix;
280 
281  public:
282  TensorInversion(material_vector_type& material_vector_, material_matrix_type& material_matrix_)
283  : material_vector(material_vector_)
284  , material_matrix(material_matrix_) {}
285 
286  KOKKOS_FORCEINLINE_FUNCTION
287  void operator()(local_ordinal_type i) const {
288  for (size_t j = 0; j < material_matrix.extent(1); ++j) {
289  for (size_t k = 0; k < material_matrix.extent(2); ++k) {
290  material_matrix(i, j, k) = material_vector(i, j * material_matrix.extent(1) + k);
291  }
292  }
293  auto matrix_material = Kokkos::subview(material_matrix, i, Kokkos::ALL(), Kokkos::ALL());
294  KokkosBatched::SerialLU<KokkosBatched::Algo::LU::Unblocked>::invoke(matrix_material);
295  }
296 };
297 
298 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
300  private:
302  using local_matrix_type = typename matrix_type::local_matrix_type;
303  using scalar_type = typename local_matrix_type::value_type;
305  using ATS = Kokkos::ArithTraits<scalar_type>;
306  using impl_scalar_type = typename ATS::val_type;
307  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
308  using magnitudeType = typename implATS::magnitudeType;
309  using magATS = Kokkos::ArithTraits<magnitudeType>;
311  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
313  using memory_space = typename local_matrix_type::memory_space;
314 
315  using local_material_type = Kokkos::View<impl_scalar_type***, memory_space>;
316  using local_dist_type = Kokkos::View<impl_scalar_type**, memory_space>;
317 
320 
323 
325 
327 
328  const scalar_type one = ATS::one();
329 
330  public:
332  coordsMV = coords_;
333 
334  auto importer = A.getCrsGraph()->getImporter();
335  if (!importer.is_null()) {
337  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
338  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
339  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
340  } else {
341  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
343  }
344 
345  {
346  Teuchos::RCP<material_type> ghostedMaterial;
347  if (!importer.is_null()) {
348  ghostedMaterial = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(importer->getTargetMap(), material_->getNumVectors());
349  ghostedMaterial->doImport(*material_, *importer, Xpetra::INSERT);
350  } else {
351  ghostedMaterial = material_;
352  }
353 
354  using execution_space = typename Node::execution_space;
355  using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
356 
357  local_ordinal_type dim = std::sqrt(material_->getNumVectors());
358  auto lclMaterial = ghostedMaterial->getLocalViewDevice(Xpetra::Access::ReadOnly);
359  material = local_material_type("material", lclMaterial.extent(0), dim, dim);
360  lcl_dist = local_dist_type("material", lclMaterial.extent(0), dim);
362  Kokkos::parallel_for("MueLu:TensorMaterialDistanceFunctor::inversion", range_type(0, lclMaterial.extent(0)), functor);
363  }
364  }
365 
366  KOKKOS_INLINE_FUNCTION
368  // || x_row - x_col ||_S^2
369  // where
370  // S = inv(material(col))
371 
372  // row material
373  impl_scalar_type d_row = implATS::zero();
374  {
375  auto matrix_row_material = Kokkos::subview(material, row, Kokkos::ALL(), Kokkos::ALL());
376  auto dist = Kokkos::subview(lcl_dist, row, Kokkos::ALL());
377 
378  for (size_t j = 0; j < coords.extent(1); ++j) {
379  dist(j) = coords(row, j) - ghostedCoords(col, j);
380  }
381 
382  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Lower, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::Unit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_row_material, dist);
383  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Upper, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::NonUnit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_row_material, dist);
384 
385  for (size_t j = 0; j < coords.extent(1); ++j) {
386  d_row += dist(j) * (coords(row, j) - ghostedCoords(col, j));
387  }
388  }
389 
390  // column material
391  impl_scalar_type d_col = implATS::zero();
392  {
393  auto matrix_col_material = Kokkos::subview(material, col, Kokkos::ALL(), Kokkos::ALL());
394  auto dist = Kokkos::subview(lcl_dist, row, Kokkos::ALL());
395 
396  for (size_t j = 0; j < coords.extent(1); ++j) {
397  dist(j) = coords(row, j) - ghostedCoords(col, j);
398  }
399 
400  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Lower, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::Unit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_col_material, dist);
401  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Upper, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::NonUnit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_col_material, dist);
402 
403  for (size_t j = 0; j < coords.extent(1); ++j) {
404  d_col += dist(j) * (coords(row, j) - ghostedCoords(col, j));
405  }
406  }
407 
408  return Kokkos::max(implATS::magnitude(d_row), implATS::magnitude(d_col));
409  }
410 };
411 
415 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
418  DistanceFunctorType& distFunctor) {
419  using scalar_type = Scalar;
420  using local_ordinal_type = LocalOrdinal;
421  using global_ordinal_type = GlobalOrdinal;
422  using node_type = Node;
423  using ATS = Kokkos::ArithTraits<scalar_type>;
424  using impl_scalar_type = typename ATS::val_type;
425  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
426  using magnitudeType = typename implATS::magnitudeType;
427  using execution_space = typename Node::execution_space;
428  using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
429 
431  {
432  auto lclA = A.getLocalMatrixDevice();
433  auto lclDiag = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
434 
435  Kokkos::parallel_for(
436  "MueLu:CoalesceDropF:Build:scalar_filter:laplacian_diag",
437  range_type(0, lclA.numRows()),
438  KOKKOS_LAMBDA(const local_ordinal_type& row) {
439  auto rowView = lclA.rowConst(row);
440  auto length = rowView.length;
441 
442  magnitudeType d;
443  impl_scalar_type d2 = implATS::zero();
444  bool haveAddedToDiag = false;
445  for (local_ordinal_type colID = 0; colID < length; colID++) {
446  auto col = rowView.colidx(colID);
447  if (row != col) {
448  d = distFunctor.distance2(row, col);
449  d2 += implATS::one() / d;
450  haveAddedToDiag = true;
451  }
452  }
453 
454  // Deal with the situation where boundary conditions have only been enforced on rows, but not on columns.
455  // We enforce dropping of these entries by assigning a very large number to the diagonal entries corresponding to BCs.
456  lclDiag(row, 0) = !haveAddedToDiag ? implATS::squareroot(implATS::rmax()) : d2;
457  });
458  }
459  auto importer = A.getCrsGraph()->getImporter();
460  if (!importer.is_null()) {
462  ghostedDiag->doImport(*diag, *importer, Xpetra::INSERT);
463  return ghostedDiag;
464  } else {
465  return diag;
466  }
467 }
468 
469 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
472  DistanceFunctorType& distFunctor) {
473  using scalar_type = Scalar;
474  using local_ordinal_type = LocalOrdinal;
475  using global_ordinal_type = GlobalOrdinal;
476  using node_type = Node;
477  using ATS = Kokkos::ArithTraits<scalar_type>;
478  using impl_scalar_type = typename ATS::val_type;
479  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
480  using magnitudeType = typename implATS::magnitudeType;
481  using execution_space = typename Node::execution_space;
482  using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
483  using mATS = Kokkos::ArithTraits<magnitudeType>;
484 
486  {
487  auto lclA = A.getLocalMatrixDevice();
488  auto lclDiag = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
489 
490  Kokkos::parallel_for(
491  "MueLu:CoalesceDropF:Build:scalar_filter:laplacian_diag",
492  range_type(0, lclA.numRows()),
493  KOKKOS_LAMBDA(const local_ordinal_type& row) {
494  auto rowView = lclA.rowConst(row);
495  auto length = rowView.length;
496 
497  impl_scalar_type mymax = implATS::zero();
498  magnitudeType d;
499  impl_scalar_type d2;
500  for (local_ordinal_type colID = 0; colID < length; colID++) {
501  auto col = rowView.colidx(colID);
502  if (row != col) {
503  d = distFunctor.distance2(row, col);
504  d2 = implATS::one() / d;
505  if (implATS::magnitude(mymax) < -implATS::magnitude(d2))
506  mymax = -implATS::magnitude(d2);
507  }
508  }
509  lclDiag(row, 0) = mymax;
510  });
511  }
512  auto importer = A.getCrsGraph()->getImporter();
513  if (!importer.is_null()) {
515  ghostedDiag->doImport(*diag, *importer, Xpetra::INSERT);
516  return ghostedDiag;
517  } else {
518  return diag;
519  }
520 }
521 
532 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType, Misc::StrengthMeasure measure>
533 class DropFunctor {
534  public:
536  using local_matrix_type = typename matrix_type::local_matrix_type;
537  using scalar_type = typename local_matrix_type::value_type;
538  using local_ordinal_type = typename local_matrix_type::ordinal_type;
539  using memory_space = typename local_matrix_type::memory_space;
541  using diag_view_type = typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
542 
543  using results_view = Kokkos::View<DecisionType*, memory_space>;
544 
545  using ATS = Kokkos::ArithTraits<scalar_type>;
546  using magnitudeType = typename ATS::magnitudeType;
547  using boundary_nodes_view = Kokkos::View<const bool*, memory_space>;
548  using mATS = Kokkos::ArithTraits<magnitudeType>;
549 
550  private:
554  diag_view_type diag; // corresponds to overlapped diagonal
555  DistanceFunctorType dist2;
557  const scalar_type one = ATS::one();
558 
559  public:
560  DropFunctor(matrix_type& A_, magnitudeType threshold, DistanceFunctorType& dist2_, results_view& results_)
561  : A(A_.getLocalMatrixDevice())
562  , eps(threshold)
563  , dist2(dist2_)
564  , results(results_) {
565  if constexpr ((measure == Misc::SmoothedAggregationMeasure) || (measure == Misc::SignedSmoothedAggregationMeasure)) {
566  diagVec = getDiagonal(A_, dist2);
567  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
568  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
569  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
571  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
572  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
573  }
574  }
575 
576  KOKKOS_FORCEINLINE_FUNCTION
577  void operator()(local_ordinal_type rlid) const {
578  auto row = A.rowConst(rlid);
579  const size_t offset = A.graph.row_map(rlid);
580  for (local_ordinal_type k = 0; k < row.length; ++k) {
581  auto clid = row.colidx(k);
582 
583  scalar_type val;
584  if (rlid != clid) {
585  val = one / dist2.distance2(rlid, clid);
586  } else {
587  val = diag(rlid);
588  }
589 
590  if constexpr (measure == Misc::SmoothedAggregationMeasure) {
591  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
592  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
593 
594  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
595  results(offset + k));
596  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
597  auto neg_aij = -ATS::real(val);
598  auto max_neg_aik = eps * ATS::real(diag(rlid));
599  results(offset + k) = Kokkos::max((neg_aij <= max_neg_aik) ? DROP : KEEP,
600  results(offset + k));
601  } else if constexpr (measure == Misc::SignedSmoothedAggregationMeasure) {
602  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
603  const bool is_nonpositive = ATS::real(val) <= mATS::zero();
604  magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
605  // + |a_ij|^2, if a_ij < 0, - |a_ij|^2 if a_ij >=0
606  if (is_nonpositive)
607  aij2 = -aij2;
608  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
609  results(offset + k));
610  }
611  }
612  }
613 };
614 
615 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType, Misc::StrengthMeasure measure>
617  public:
619  using local_matrix_type = typename matrix_type::local_matrix_type;
620  using scalar_type = typename local_matrix_type::value_type;
621  using local_ordinal_type = typename local_matrix_type::ordinal_type;
622  using memory_space = typename local_matrix_type::memory_space;
623  using block_indices_view_type = Kokkos::View<local_ordinal_type*, memory_space>;
625  using diag_view_type = typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
626 
627  using results_view = Kokkos::View<DecisionType*, memory_space>;
628 
629  using ATS = Kokkos::ArithTraits<scalar_type>;
630  using magnitudeType = typename ATS::magnitudeType;
631  using boundary_nodes_view = Kokkos::View<const bool*, memory_space>;
632  using mATS = Kokkos::ArithTraits<magnitudeType>;
633 
634  private:
638  diag_view_type diag; // corresponds to overlapped diagonal
639  DistanceFunctorType dist2;
643  const scalar_type one = ATS::one();
644 
645  public:
646  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_)
647  : A(A_.getLocalMatrixDevice())
648  , eps(threshold)
649  , dist2(dist2_)
650  , results(results_)
651  , point_to_block(point_to_block_)
652  , ghosted_point_to_block(ghosted_point_to_block_) {
653  if constexpr ((measure == Misc::SmoothedAggregationMeasure) || (measure == Misc::SignedSmoothedAggregationMeasure)) {
654  diagVec = getDiagonal(mergedA_, dist2);
655  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
656  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
657  } else if (measure == Misc::SignedRugeStuebenMeasure) {
659  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
660  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
661  }
662  }
663 
664  KOKKOS_FORCEINLINE_FUNCTION
665  void operator()(local_ordinal_type rlid) const {
666  auto brlid = point_to_block(rlid);
667  auto row = A.rowConst(rlid);
668  const size_t offset = A.graph.row_map(rlid);
669  for (local_ordinal_type k = 0; k < row.length; ++k) {
670  auto clid = row.colidx(k);
671  auto bclid = ghosted_point_to_block(clid);
672 
673  scalar_type val;
674  if (brlid != bclid) {
675  val = one / dist2.distance2(brlid, bclid);
676  } else {
677  val = diag(brlid);
678  }
679 
680  if constexpr (measure == Misc::SmoothedAggregationMeasure) {
681  auto aiiajj = ATS::magnitude(diag(brlid)) * ATS::magnitude(diag(bclid)); // |a_ii|*|a_jj|
682  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
683 
684  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
685  results(offset + k));
686  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
687  auto neg_aij = -ATS::real(val);
688  auto max_neg_aik = eps * ATS::real(diag(brlid));
689  results(offset + k) = Kokkos::max((neg_aij <= max_neg_aik) ? DROP : KEEP,
690  results(offset + k));
691  } else if constexpr (measure == Misc::SignedSmoothedAggregationMeasure) {
692  auto aiiajj = ATS::magnitude(diag(brlid)) * ATS::magnitude(diag(bclid)); // |a_ii|*|a_jj|
693  const bool is_nonpositive = ATS::real(val) <= mATS::zero();
694  magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
695  // + |a_ij|^2, if a_ij < 0, - |a_ij|^2 if a_ij >=0
696  if (is_nonpositive)
697  aij2 = -aij2;
698  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
699  results(offset + k));
700  }
701  }
702  }
703 };
704 
705 // helper function to allow partial template deduction
706 template <Misc::StrengthMeasure measure, class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
709  DistanceFunctorType& dist2_,
711  auto functor = DropFunctor<Scalar, LocalOrdinal, GlobalOrdinal, Node, DistanceFunctorType, measure>(A_, threshold, dist2_, results_);
712  return functor;
713 }
714 
715 template <Misc::StrengthMeasure measure, class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
719  DistanceFunctorType& dist2_,
723  auto functor = VectorDropFunctor<Scalar, LocalOrdinal, GlobalOrdinal, Node, DistanceFunctorType, measure>(A_, mergedA_, threshold, dist2_, results_, point_to_block_, ghosted_point_to_block_);
724  return functor;
725 }
726 
727 } // namespace MueLu::DistanceLaplacian
728 
729 #endif
Drops entries the unscaled distance Laplacian.
typename material_type::dual_view_type_const::t_dev local_material_type
WeightedDistanceFunctor(matrix_type &A, Teuchos::RCP< coords_type > &coords_, weight_type weight_)
MueLu::DefaultLocalOrdinal LocalOrdinal
UnweightedDistanceFunctor(matrix_type &A, Teuchos::RCP< coords_type > &coords_)
Kokkos::View< local_ordinal_type *, memory_space > block_indices_view_type
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
typename local_matrix_type::memory_space memory_space
typename Kokkos::DualView< const scalar_type *, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged >::t_dev diag_view_type
typename coords_type::dual_view_type_const::t_dev local_coords_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
Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getMaxMinusOffDiagonal(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, DistanceFunctorType &distFunctor)
Kokkos::View< DecisionType *, memory_space > results_view
BlockWeightedDistanceFunctor(matrix_type &A, Teuchos::RCP< coords_type > &coords_, weight_type weight_, local_ordinal_type interleaved_blocksize_)
MueLu::DefaultNode Node
typename local_matrix_type::ordinal_type local_ordinal_type
Kokkos::View< const bool *, memory_space > boundary_nodes_view
DropFunctor(matrix_type &A_, magnitudeType threshold, DistanceFunctorType &dist2_, results_view &results_)
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_)
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
Kokkos::View< DecisionType *, memory_space > results_view
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
TensorInversion(material_vector_type &material_vector_, material_matrix_type &material_matrix_)
typename local_matrix_type::value_type scalar_type
typename Kokkos::DualView< const scalar_type *, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged >::t_dev diag_view_type
typename matrix_type::local_matrix_type local_matrix_type
typename coords_type::dual_view_type_const::t_dev local_coords_type
auto make_drop_functor(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A_, typename DropFunctor< Scalar, LocalOrdinal, GlobalOrdinal, Node, DistanceFunctorType, measure >::magnitudeType threshold, DistanceFunctorType &dist2_, typename DropFunctor< Scalar, LocalOrdinal, GlobalOrdinal, Node, DistanceFunctorType, measure >::results_view &results_)
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_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
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
typename coords_type::dual_view_type_const::t_dev local_coords_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 local_matrix_type::ordinal_type local_ordinal_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
#define TEUCHOS_ASSERT(assertion_test)
virtual size_t getNumVectors() const =0
typename local_matrix_type::memory_space memory_space
auto make_vector_drop_functor(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A_, Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &mergedA_, typename VectorDropFunctor< Scalar, LocalOrdinal, GlobalOrdinal, Node, DistanceFunctorType, measure >::magnitudeType threshold, DistanceFunctorType &dist2_, typename VectorDropFunctor< Scalar, LocalOrdinal, GlobalOrdinal, Node, DistanceFunctorType, measure >::results_view &results_, typename VectorDropFunctor< Scalar, LocalOrdinal, GlobalOrdinal, Node, DistanceFunctorType, measure >::block_indices_view_type point_to_block_, typename VectorDropFunctor< Scalar, LocalOrdinal, GlobalOrdinal, Node, DistanceFunctorType, measure >::block_indices_view_type ghosted_point_to_block_)
KOKKOS_FORCEINLINE_FUNCTION magnitudeType distance2(const local_ordinal_type row, const local_ordinal_type col) const
typename matrix_type::local_matrix_type local_matrix_type