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 #if KOKKOS_VERSION >= 40799
21 #include "KokkosKernels_ArithTraits.hpp"
22 #else
23 #include "Kokkos_ArithTraits.hpp"
24 #endif
25 #include "Teuchos_RCP.hpp"
26 #include "Xpetra_Matrix.hpp"
27 #include "Xpetra_MultiVector.hpp"
28 #include "Xpetra_MultiVectorFactory.hpp"
29 
30 namespace MueLu::DistanceLaplacian {
31 
36 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
38  private:
40  using local_matrix_type = typename matrix_type::local_matrix_type;
41  using scalar_type = typename local_matrix_type::value_type;
43 #if KOKKOS_VERSION >= 40799
44  using ATS = KokkosKernels::ArithTraits<scalar_type>;
45 #else
46  using ATS = Kokkos::ArithTraits<scalar_type>;
47 #endif
48  using impl_scalar_type = typename ATS::val_type;
49 #if KOKKOS_VERSION >= 40799
50  using implATS = KokkosKernels::ArithTraits<impl_scalar_type>;
51 #else
52  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
53 #endif
54  using magnitudeType = typename implATS::magnitudeType;
55 #if KOKKOS_VERSION >= 40799
56  using magATS = KokkosKernels::ArithTraits<magnitudeType>;
57 #else
58  using magATS = Kokkos::ArithTraits<magnitudeType>;
59 #endif
61  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
62 
65 
68 
69  public:
71  coordsMV = coords_;
72  auto importer = A.getCrsGraph()->getImporter();
73  if (!importer.is_null()) {
75  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
76  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
77  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
78  } else {
79  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
81  }
82  }
83 
84  KOKKOS_FORCEINLINE_FUNCTION
86  magnitudeType d = magATS::zero();
87  magnitudeType s;
88  for (size_t j = 0; j < coords.extent(1); ++j) {
89  s = coords(row, j) - ghostedCoords(col, j);
90  d += s * s;
91  }
92  return d;
93  }
94 };
95 
100 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class weight_type>
102  private:
104  using local_matrix_type = typename matrix_type::local_matrix_type;
105  using scalar_type = typename local_matrix_type::value_type;
107 #if KOKKOS_VERSION >= 40799
108  using ATS = KokkosKernels::ArithTraits<scalar_type>;
109 #else
110  using ATS = Kokkos::ArithTraits<scalar_type>;
111 #endif
112  using impl_scalar_type = typename ATS::val_type;
113 #if KOKKOS_VERSION >= 40799
114  using implATS = KokkosKernels::ArithTraits<impl_scalar_type>;
115 #else
116  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
117 #endif
118  using magnitudeType = typename implATS::magnitudeType;
119 #if KOKKOS_VERSION >= 40799
120  using magATS = KokkosKernels::ArithTraits<magnitudeType>;
121 #else
122  using magATS = Kokkos::ArithTraits<magnitudeType>;
123 #endif
125  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
126 
129 
132 
133  weight_type weight;
134 
135  public:
136  WeightedDistanceFunctor(matrix_type& A, Teuchos::RCP<coords_type>& coords_, weight_type weight_) {
137  coordsMV = coords_;
138  auto importer = A.getCrsGraph()->getImporter();
139  if (!importer.is_null()) {
141  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
142  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
143  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
144  } else {
145  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
147  }
148  weight = weight_;
149  TEUCHOS_ASSERT(weight.extent(0) == coordsMV->getNumVectors());
150  }
151 
152  KOKKOS_FORCEINLINE_FUNCTION
154  magnitudeType d = magATS::zero();
155  magnitudeType w;
156  magnitudeType s;
157  for (size_t j = 0; j < coords.extent(1); ++j) {
158  s = coords(row, j) - ghostedCoords(col, j);
159  w = weight(j);
160  d += w * s * s;
161  }
162  return d;
163  }
164 };
165 
170 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class weight_type>
172  private:
174  using local_matrix_type = typename matrix_type::local_matrix_type;
175  using scalar_type = typename local_matrix_type::value_type;
177 #if KOKKOS_VERSION >= 40799
178  using ATS = KokkosKernels::ArithTraits<scalar_type>;
179 #else
180  using ATS = Kokkos::ArithTraits<scalar_type>;
181 #endif
182  using impl_scalar_type = typename ATS::val_type;
183 #if KOKKOS_VERSION >= 40799
184  using implATS = KokkosKernels::ArithTraits<impl_scalar_type>;
185 #else
186  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
187 #endif
188  using magnitudeType = typename implATS::magnitudeType;
189 #if KOKKOS_VERSION >= 40799
190  using magATS = KokkosKernels::ArithTraits<magnitudeType>;
191 #else
192  using magATS = Kokkos::ArithTraits<magnitudeType>;
193 #endif
195  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
196 
199 
202 
203  weight_type weight;
205 
206  public:
207  BlockWeightedDistanceFunctor(matrix_type& A, Teuchos::RCP<coords_type>& coords_, weight_type weight_, local_ordinal_type interleaved_blocksize_) {
208  coordsMV = coords_;
209  auto importer = A.getCrsGraph()->getImporter();
210  if (!importer.is_null()) {
212  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
213  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
214  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
215  } else {
216  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
218  }
219  weight = weight_;
220  interleaved_blocksize = interleaved_blocksize_;
221  TEUCHOS_ASSERT(weight.extent(0) == coordsMV->getNumVectors());
222  }
223 
224  KOKKOS_FORCEINLINE_FUNCTION
226  magnitudeType d = magATS::zero();
227  magnitudeType w;
228  magnitudeType s;
229  local_ordinal_type block_id = row % interleaved_blocksize;
230  local_ordinal_type block_start = block_id * interleaved_blocksize;
231  for (size_t j = 0; j < coords.extent(1); ++j) {
232  s = coords(row, j) - ghostedCoords(col, j);
233  w = weight(block_start + j);
234  d += w * s * s;
235  }
236  return d;
237  }
238 };
239 
240 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
242  private:
244  using local_matrix_type = typename matrix_type::local_matrix_type;
245  using scalar_type = typename local_matrix_type::value_type;
247 #if KOKKOS_VERSION >= 40799
248  using ATS = KokkosKernels::ArithTraits<scalar_type>;
249 #else
250  using ATS = Kokkos::ArithTraits<scalar_type>;
251 #endif
252  using impl_scalar_type = typename ATS::val_type;
253 #if KOKKOS_VERSION >= 40799
254  using implATS = KokkosKernels::ArithTraits<impl_scalar_type>;
255 #else
256  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
257 #endif
258  using magnitudeType = typename implATS::magnitudeType;
259 #if KOKKOS_VERSION >= 40799
260  using magATS = KokkosKernels::ArithTraits<magnitudeType>;
261 #else
262  using magATS = Kokkos::ArithTraits<magnitudeType>;
263 #endif
265  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
267  using local_material_type = typename material_type::dual_view_type_const::t_dev;
268 
271 
274 
277 
280 
281  public:
283  coordsMV = coords_;
284  materialMV = material_;
285  auto importer = A.getCrsGraph()->getImporter();
286  if (!importer.is_null()) {
288  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
289  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
290  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
291 
293  ghostedMaterialMV->doImport(*materialMV, *importer, Xpetra::INSERT);
294  material = materialMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
295  ghostedMaterial = ghostedMaterialMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
296  } else {
297  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
299 
300  material = materialMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
302  }
303  }
304 
305  KOKKOS_INLINE_FUNCTION
307  // || x_row - x_col ||_S^2
308  // where
309  // S = 1/material(row) * Identity
310  magnitudeType d = magATS::zero();
311  magnitudeType d_row = magATS::zero();
312  magnitudeType d_col = magATS::zero();
313  magnitudeType s;
314 
315  for (size_t j = 0; j < coords.extent(1); ++j) {
316  s = coords(row, j) - ghostedCoords(col, j);
317  d += s * s;
318  }
319 
320  d_row = d / implATS::magnitude(ghostedMaterial(row, 0));
321  d_col = d / implATS::magnitude(ghostedMaterial(col, 0));
322 
323  return Kokkos::max(d_row, d_col);
324  }
325 };
326 
327 template <class local_ordinal_type, class material_vector_type, class material_matrix_type>
329  private:
330  material_vector_type material_vector;
331  material_matrix_type material_matrix;
332 
333  public:
334  TensorInversion(material_vector_type& material_vector_, material_matrix_type& material_matrix_)
335  : material_vector(material_vector_)
336  , material_matrix(material_matrix_) {}
337 
338  KOKKOS_FORCEINLINE_FUNCTION
339  void operator()(local_ordinal_type i) const {
340  for (size_t j = 0; j < material_matrix.extent(1); ++j) {
341  for (size_t k = 0; k < material_matrix.extent(2); ++k) {
342  material_matrix(i, j, k) = material_vector(i, j * material_matrix.extent(1) + k);
343  }
344  }
345  auto matrix_material = Kokkos::subview(material_matrix, i, Kokkos::ALL(), Kokkos::ALL());
346  KokkosBatched::SerialLU<KokkosBatched::Algo::LU::Unblocked>::invoke(matrix_material);
347  }
348 };
349 
350 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
352  private:
354  using local_matrix_type = typename matrix_type::local_matrix_type;
355  using scalar_type = typename local_matrix_type::value_type;
357 #if KOKKOS_VERSION >= 40799
358  using ATS = KokkosKernels::ArithTraits<scalar_type>;
359 #else
360  using ATS = Kokkos::ArithTraits<scalar_type>;
361 #endif
362  using impl_scalar_type = typename ATS::val_type;
363 #if KOKKOS_VERSION >= 40799
364  using implATS = KokkosKernels::ArithTraits<impl_scalar_type>;
365 #else
366  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
367 #endif
368  using magnitudeType = typename implATS::magnitudeType;
369 #if KOKKOS_VERSION >= 40799
370  using magATS = KokkosKernels::ArithTraits<magnitudeType>;
371 #else
372  using magATS = Kokkos::ArithTraits<magnitudeType>;
373 #endif
375  using local_coords_type = typename coords_type::dual_view_type_const::t_dev;
377  using memory_space = typename local_matrix_type::memory_space;
378 
379  using local_material_type = Kokkos::View<impl_scalar_type***, memory_space>;
380  using local_dist_type = Kokkos::View<impl_scalar_type**, memory_space>;
381 
384 
387 
389 
391 
392  const scalar_type one = ATS::one();
393 
394  public:
396  coordsMV = coords_;
397 
398  auto importer = A.getCrsGraph()->getImporter();
399  if (!importer.is_null()) {
401  ghostedCoordsMV->doImport(*coordsMV, *importer, Xpetra::INSERT);
402  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
403  ghostedCoords = ghostedCoordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
404  } else {
405  coords = coordsMV->getLocalViewDevice(Xpetra::Access::ReadOnly);
407  }
408 
409  {
410  Teuchos::RCP<material_type> ghostedMaterial;
411  if (!importer.is_null()) {
412  ghostedMaterial = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(importer->getTargetMap(), material_->getNumVectors());
413  ghostedMaterial->doImport(*material_, *importer, Xpetra::INSERT);
414  } else {
415  ghostedMaterial = material_;
416  }
417 
418  using execution_space = typename Node::execution_space;
419  using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
420 
421  local_ordinal_type dim = std::sqrt(material_->getNumVectors());
422  auto lclMaterial = ghostedMaterial->getLocalViewDevice(Xpetra::Access::ReadOnly);
423  material = local_material_type("material", lclMaterial.extent(0), dim, dim);
424  lcl_dist = local_dist_type("material", lclMaterial.extent(0), dim);
426  Kokkos::parallel_for("MueLu:TensorMaterialDistanceFunctor::inversion", range_type(0, lclMaterial.extent(0)), functor);
427  }
428  }
429 
430  KOKKOS_INLINE_FUNCTION
432  // || x_row - x_col ||_S^2
433  // where
434  // S = inv(material(col))
435 
436  // row material
437  impl_scalar_type d_row = implATS::zero();
438  {
439  auto matrix_row_material = Kokkos::subview(material, row, Kokkos::ALL(), Kokkos::ALL());
440  auto dist = Kokkos::subview(lcl_dist, row, Kokkos::ALL());
441 
442  for (size_t j = 0; j < coords.extent(1); ++j) {
443  dist(j) = coords(row, j) - ghostedCoords(col, j);
444  }
445 
446  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Lower, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::Unit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_row_material, dist);
447  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Upper, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::NonUnit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_row_material, dist);
448 
449  for (size_t j = 0; j < coords.extent(1); ++j) {
450  d_row += dist(j) * (coords(row, j) - ghostedCoords(col, j));
451  }
452  }
453 
454  // column material
455  impl_scalar_type d_col = implATS::zero();
456  {
457  auto matrix_col_material = Kokkos::subview(material, col, Kokkos::ALL(), Kokkos::ALL());
458  auto dist = Kokkos::subview(lcl_dist, row, Kokkos::ALL());
459 
460  for (size_t j = 0; j < coords.extent(1); ++j) {
461  dist(j) = coords(row, j) - ghostedCoords(col, j);
462  }
463 
464  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Lower, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::Unit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_col_material, dist);
465  KokkosBatched::SerialTrsv<KokkosBatched::Uplo::Upper, KokkosBatched::Trans::NoTranspose, KokkosBatched::Diag::NonUnit, KokkosBatched::Algo::Trsv::Unblocked>::invoke(one, matrix_col_material, dist);
466 
467  for (size_t j = 0; j < coords.extent(1); ++j) {
468  d_col += dist(j) * (coords(row, j) - ghostedCoords(col, j));
469  }
470  }
471 
472  return Kokkos::max(implATS::magnitude(d_row), implATS::magnitude(d_col));
473  }
474 };
475 
479 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
482  DistanceFunctorType& distFunctor) {
483  using scalar_type = Scalar;
484  using local_ordinal_type = LocalOrdinal;
485 #if KOKKOS_VERSION >= 40799
486  using ATS = KokkosKernels::ArithTraits<scalar_type>;
487 #else
488  using ATS = Kokkos::ArithTraits<scalar_type>;
489 #endif
490  using impl_scalar_type = typename ATS::val_type;
491 #if KOKKOS_VERSION >= 40799
492  using implATS = KokkosKernels::ArithTraits<impl_scalar_type>;
493 #else
494  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
495 #endif
496  using magnitudeType = typename implATS::magnitudeType;
497  using execution_space = typename Node::execution_space;
498  using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
499 
501  {
502  auto lclA = A.getLocalMatrixDevice();
503  auto lclDiag = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
504 
505  Kokkos::parallel_for(
506  "MueLu:CoalesceDropF:Build:scalar_filter:laplacian_diag",
507  range_type(0, lclA.numRows()),
508  KOKKOS_LAMBDA(const local_ordinal_type& row) {
509  auto rowView = lclA.rowConst(row);
510  auto length = rowView.length;
511 
512  magnitudeType d;
513  impl_scalar_type d2 = implATS::zero();
514  bool haveAddedToDiag = false;
515  for (local_ordinal_type colID = 0; colID < length; colID++) {
516  auto col = rowView.colidx(colID);
517  if (row != col) {
518  d = distFunctor.distance2(row, col);
519  d2 += implATS::one() / d;
520  haveAddedToDiag = true;
521  }
522  }
523 
524  // Deal with the situation where boundary conditions have only been enforced on rows, but not on columns.
525  // We enforce dropping of these entries by assigning a very large number to the diagonal entries corresponding to BCs.
526  lclDiag(row, 0) = !haveAddedToDiag ? implATS::squareroot(implATS::rmax()) : d2;
527  });
528  }
529  auto importer = A.getCrsGraph()->getImporter();
530  if (!importer.is_null()) {
532  ghostedDiag->doImport(*diag, *importer, Xpetra::INSERT);
533  return ghostedDiag;
534  } else {
535  return diag;
536  }
537 }
538 
539 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
542  DistanceFunctorType& distFunctor) {
543  using scalar_type = Scalar;
544  using local_ordinal_type = LocalOrdinal;
545 #if KOKKOS_VERSION >= 40799
546  using ATS = KokkosKernels::ArithTraits<scalar_type>;
547 #else
548  using ATS = Kokkos::ArithTraits<scalar_type>;
549 #endif
550  using impl_scalar_type = typename ATS::val_type;
551 #if KOKKOS_VERSION >= 40799
552  using implATS = KokkosKernels::ArithTraits<impl_scalar_type>;
553 #else
554  using implATS = Kokkos::ArithTraits<impl_scalar_type>;
555 #endif
556  using magnitudeType = typename implATS::magnitudeType;
557  using execution_space = typename Node::execution_space;
558  using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
559 
561  {
562  auto lclA = A.getLocalMatrixDevice();
563  auto lclDiag = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
564 
565  Kokkos::parallel_for(
566  "MueLu:CoalesceDropF:Build:scalar_filter:laplacian_diag",
567  range_type(0, lclA.numRows()),
568  KOKKOS_LAMBDA(const local_ordinal_type& row) {
569  auto rowView = lclA.rowConst(row);
570  auto length = rowView.length;
571 
572  impl_scalar_type mymax = implATS::zero();
573  magnitudeType d;
574  impl_scalar_type d2;
575  for (local_ordinal_type colID = 0; colID < length; colID++) {
576  auto col = rowView.colidx(colID);
577  if (row != col) {
578  d = distFunctor.distance2(row, col);
579  d2 = implATS::one() / d;
580  if (implATS::magnitude(mymax) < implATS::magnitude(d2))
581  mymax = implATS::magnitude(d2);
582  }
583  }
584  lclDiag(row, 0) = mymax;
585  });
586  }
587  auto importer = A.getCrsGraph()->getImporter();
588  if (!importer.is_null()) {
590  ghostedDiag->doImport(*diag, *importer, Xpetra::INSERT);
591  return ghostedDiag;
592  } else {
593  return diag;
594  }
595 }
596 
607 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType, Misc::StrengthMeasure measure>
608 class DropFunctor {
609  public:
611  using local_matrix_type = typename matrix_type::local_matrix_type;
612  using scalar_type = typename local_matrix_type::value_type;
613  using local_ordinal_type = typename local_matrix_type::ordinal_type;
614  using memory_space = typename local_matrix_type::memory_space;
616  using diag_view_type = typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
617 
618  using results_view = Kokkos::View<DecisionType*, memory_space>;
619 
620 #if KOKKOS_VERSION >= 40799
621  using ATS = KokkosKernels::ArithTraits<scalar_type>;
622 #else
623  using ATS = Kokkos::ArithTraits<scalar_type>;
624 #endif
625  using magnitudeType = typename ATS::magnitudeType;
626  using boundary_nodes_view = Kokkos::View<const bool*, memory_space>;
627 #if KOKKOS_VERSION >= 40799
628  using mATS = KokkosKernels::ArithTraits<magnitudeType>;
629 #else
630  using mATS = Kokkos::ArithTraits<magnitudeType>;
631 #endif
632 
633  private:
637  diag_view_type diag; // corresponds to overlapped diagonal
638  DistanceFunctorType dist2;
640  const scalar_type one = ATS::one();
641 
642  public:
643  DropFunctor(matrix_type& A_, magnitudeType threshold, DistanceFunctorType& dist2_, results_view& results_)
644  : A(A_.getLocalMatrixDevice())
645  , eps(threshold)
646  , dist2(dist2_)
647  , results(results_) {
648  if constexpr ((measure == Misc::SmoothedAggregationMeasure) || (measure == Misc::SignedSmoothedAggregationMeasure)) {
649  diagVec = getDiagonal(A_, dist2);
650  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
651  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
652  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
654  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
655  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
656  }
657  }
658 
659  KOKKOS_FORCEINLINE_FUNCTION
660  void operator()(local_ordinal_type rlid) const {
661  auto row = A.rowConst(rlid);
662  const size_t offset = A.graph.row_map(rlid);
663 
664 #ifdef MUELU_COALESCE_DROP_DEBUG
665  {
666  Kokkos::printf("SoC: ");
667  for (local_ordinal_type k = 0; k < row.length; ++k) {
668  auto clid = row.colidx(k);
669 
670  scalar_type val;
671  if (rlid != clid) {
672  val = -one / dist2.distance2(rlid, clid);
673  } else {
674  val = diag(rlid);
675  }
676 
677  if constexpr (measure == Misc::SmoothedAggregationMeasure) {
678  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
679  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
680 
681  Kokkos::printf("%5f ", ATS::sqrt(aij2 / aiiajj));
682  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
683  auto neg_aij = -ATS::real(val);
684  auto max_neg_aik = ATS::real(diag(rlid));
685  Kokkos::printf("%5f ", neg_aij / max_neg_aik);
686  } else if constexpr (measure == Misc::SignedSmoothedAggregationMeasure) {
687  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
688  const bool is_nonpositive = ATS::real(val) <= mATS::zero();
689  magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
690  // + |a_ij|^2, if a_ij < 0, - |a_ij|^2 if a_ij >=0
691  if (!is_nonpositive)
692  aij2 = -aij2;
693  Kokkos::printf("%5f ", ATS::sqrt(aij2 / aiiajj));
694  }
695  }
696  Kokkos::printf("\n");
697  }
698 #endif
699 
700  for (local_ordinal_type k = 0; k < row.length; ++k) {
701  auto clid = row.colidx(k);
702 
703  scalar_type val;
704  if (rlid != clid) {
705  val = -one / dist2.distance2(rlid, clid);
706  } else {
707  val = diag(rlid);
708  }
709 
710  if constexpr (measure == Misc::SmoothedAggregationMeasure) {
711  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
712  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
713 
714  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
715  results(offset + k));
716  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
717  auto neg_aij = -ATS::real(val);
718  auto max_neg_aik = eps * ATS::real(diag(rlid));
719  results(offset + k) = Kokkos::max((neg_aij < max_neg_aik) ? DROP : KEEP,
720  results(offset + k));
721  } else if constexpr (measure == Misc::SignedSmoothedAggregationMeasure) {
722  auto aiiajj = ATS::magnitude(diag(rlid)) * ATS::magnitude(diag(clid)); // |a_ii|*|a_jj|
723  const bool is_nonpositive = ATS::real(val) <= mATS::zero();
724  magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
725  // + |a_ij|^2, if a_ij < 0, - |a_ij|^2 if a_ij >=0
726  if (!is_nonpositive)
727  aij2 = -aij2;
728  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
729  results(offset + k));
730  }
731  }
732  }
733 };
734 
735 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType, Misc::StrengthMeasure measure>
737  public:
739  using local_matrix_type = typename matrix_type::local_matrix_type;
740  using scalar_type = typename local_matrix_type::value_type;
741  using local_ordinal_type = typename local_matrix_type::ordinal_type;
742  using memory_space = typename local_matrix_type::memory_space;
743  using block_indices_view_type = Kokkos::View<local_ordinal_type*, memory_space>;
745  using diag_view_type = typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
746 
747  using results_view = Kokkos::View<DecisionType*, memory_space>;
748 
749 #if KOKKOS_VERSION >= 40799
750  using ATS = KokkosKernels::ArithTraits<scalar_type>;
751 #else
752  using ATS = Kokkos::ArithTraits<scalar_type>;
753 #endif
754  using magnitudeType = typename ATS::magnitudeType;
755  using boundary_nodes_view = Kokkos::View<const bool*, memory_space>;
756 #if KOKKOS_VERSION >= 40799
757  using mATS = KokkosKernels::ArithTraits<magnitudeType>;
758 #else
759  using mATS = Kokkos::ArithTraits<magnitudeType>;
760 #endif
761 
762  private:
766  diag_view_type diag; // corresponds to overlapped diagonal
767  DistanceFunctorType dist2;
771  const scalar_type one = ATS::one();
772 
773  public:
774  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_)
775  : A(A_.getLocalMatrixDevice())
776  , eps(threshold)
777  , dist2(dist2_)
778  , results(results_)
779  , point_to_block(point_to_block_)
780  , ghosted_point_to_block(ghosted_point_to_block_) {
781  if constexpr ((measure == Misc::SmoothedAggregationMeasure) || (measure == Misc::SignedSmoothedAggregationMeasure)) {
782  diagVec = getDiagonal(mergedA_, dist2);
783  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
784  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
785  } else if (measure == Misc::SignedRugeStuebenMeasure) {
787  auto lclDiag2d = diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
788  diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
789  }
790  }
791 
792  KOKKOS_FORCEINLINE_FUNCTION
793  void operator()(local_ordinal_type rlid) const {
794  auto brlid = point_to_block(rlid);
795  auto row = A.rowConst(rlid);
796  const size_t offset = A.graph.row_map(rlid);
797 
798 #ifdef MUELU_COALESCE_DROP_DEBUG
799  {
800  Kokkos::printf("SoC: ");
801  for (local_ordinal_type k = 0; k < row.length; ++k) {
802  auto clid = row.colidx(k);
803  auto bclid = ghosted_point_to_block(clid);
804 
805  scalar_type val;
806  if (brlid != bclid) {
807  val = -one / dist2.distance2(brlid, bclid);
808  } else {
809  val = diag(brlid);
810  }
811 
812  if constexpr (measure == Misc::SmoothedAggregationMeasure) {
813  auto aiiajj = ATS::magnitude(diag(brlid)) * ATS::magnitude(diag(bclid)); // |a_ii|*|a_jj|
814  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
815 
816  Kokkos::printf("%5f ", ATS::sqrt(aij2 / aiiajj));
817  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
818  auto neg_aij = -ATS::real(val);
819  auto max_neg_aik = eps * ATS::real(diag(brlid));
820  Kokkos::printf("%5f ", neg_aij / max_neg_aik);
821  } else if constexpr (measure == Misc::SignedSmoothedAggregationMeasure) {
822  auto aiiajj = ATS::magnitude(diag(brlid)) * ATS::magnitude(diag(bclid)); // |a_ii|*|a_jj|
823  const bool is_nonpositive = ATS::real(val) <= mATS::zero();
824  magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
825  // + |a_ij|^2, if a_ij < 0, - |a_ij|^2 if a_ij >=0
826  if (!is_nonpositive)
827  aij2 = -aij2;
828  Kokkos::printf("%5f ", ATS::sqrt(aij2 / aiiajj));
829  }
830  }
831  Kokkos::printf("\n");
832  }
833 #endif
834 
835  for (local_ordinal_type k = 0; k < row.length; ++k) {
836  auto clid = row.colidx(k);
837  auto bclid = ghosted_point_to_block(clid);
838 
839  scalar_type val;
840  if (brlid != bclid) {
841  val = -one / dist2.distance2(brlid, bclid);
842  } else {
843  val = diag(brlid);
844  }
845 
846  if constexpr (measure == Misc::SmoothedAggregationMeasure) {
847  auto aiiajj = ATS::magnitude(diag(brlid)) * ATS::magnitude(diag(bclid)); // |a_ii|*|a_jj|
848  auto aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
849 
850  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
851  results(offset + k));
852  } else if constexpr (measure == Misc::SignedRugeStuebenMeasure) {
853  auto neg_aij = -ATS::real(val);
854  auto max_neg_aik = eps * ATS::real(diag(brlid));
855  results(offset + k) = Kokkos::max((neg_aij < max_neg_aik) ? DROP : KEEP,
856  results(offset + k));
857  } else if constexpr (measure == Misc::SignedSmoothedAggregationMeasure) {
858  auto aiiajj = ATS::magnitude(diag(brlid)) * ATS::magnitude(diag(bclid)); // |a_ii|*|a_jj|
859  const bool is_nonpositive = ATS::real(val) <= mATS::zero();
860  magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val); // |a_ij|^2
861  // + |a_ij|^2, if a_ij < 0, - |a_ij|^2 if a_ij >=0
862  if (!is_nonpositive)
863  aij2 = -aij2;
864  results(offset + k) = Kokkos::max((aij2 <= eps * eps * aiiajj) ? DROP : KEEP,
865  results(offset + k));
866  }
867  }
868  }
869 };
870 
871 // helper function to allow partial template deduction
872 template <Misc::StrengthMeasure measure, class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
875  DistanceFunctorType& dist2_,
877  auto functor = DropFunctor<Scalar, LocalOrdinal, GlobalOrdinal, Node, DistanceFunctorType, measure>(A_, threshold, dist2_, results_);
878  return functor;
879 }
880 
881 template <Misc::StrengthMeasure measure, class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class DistanceFunctorType>
885  DistanceFunctorType& dist2_,
889  auto functor = VectorDropFunctor<Scalar, LocalOrdinal, GlobalOrdinal, Node, DistanceFunctorType, measure>(A_, mergedA_, threshold, dist2_, results_, point_to_block_, ghosted_point_to_block_);
890  return functor;
891 }
892 
893 } // namespace MueLu::DistanceLaplacian
894 
895 #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_)
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
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