10 #ifndef MUELU_DISTANCELAPLACIANDROPPING_HPP
11 #define MUELU_DISTANCELAPLACIANDROPPING_HPP
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"
19 #include "Kokkos_Core.hpp"
20 #include "Kokkos_ArithTraits.hpp"
24 #include "Xpetra_MultiVectorFactory.hpp"
26 namespace MueLu::DistanceLaplacian {
32 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
39 using ATS = Kokkos::ArithTraits<scalar_type>;
41 using implATS = Kokkos::ArithTraits<impl_scalar_type>;
43 using magATS = Kokkos::ArithTraits<magnitudeType>;
57 if (!importer.is_null()) {
68 KOKKOS_FORCEINLINE_FUNCTION
72 for (
size_t j = 0; j <
coords.extent(1); ++j) {
84 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class weight_type>
91 using ATS = Kokkos::ArithTraits<scalar_type>;
93 using implATS = Kokkos::ArithTraits<impl_scalar_type>;
95 using magATS = Kokkos::ArithTraits<magnitudeType>;
111 if (!importer.is_null()) {
124 KOKKOS_FORCEINLINE_FUNCTION
129 for (
size_t j = 0; j <
coords.extent(1); ++j) {
142 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class weight_type>
149 using ATS = Kokkos::ArithTraits<scalar_type>;
151 using implATS = Kokkos::ArithTraits<impl_scalar_type>;
153 using magATS = Kokkos::ArithTraits<magnitudeType>;
170 if (!importer.is_null()) {
184 KOKKOS_FORCEINLINE_FUNCTION
191 for (
size_t j = 0; j <
coords.extent(1); ++j) {
193 w =
weight(block_start + j);
200 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
207 using ATS = Kokkos::ArithTraits<scalar_type>;
209 using implATS = Kokkos::ArithTraits<impl_scalar_type>;
211 using magATS = Kokkos::ArithTraits<magnitudeType>;
234 if (!importer.is_null()) {
253 KOKKOS_INLINE_FUNCTION
263 for (
size_t j = 0; j <
coords.extent(1); ++j) {
271 return Kokkos::max(d_row, d_col);
275 template <
class local_ordinal_type,
class material_vector_type,
class material_matrix_type>
282 TensorInversion(material_vector_type& material_vector_, material_matrix_type& material_matrix_)
286 KOKKOS_FORCEINLINE_FUNCTION
293 auto matrix_material = Kokkos::subview(
material_matrix, i, Kokkos::ALL(), Kokkos::ALL());
294 KokkosBatched::SerialLU<KokkosBatched::Algo::LU::Unblocked>::invoke(matrix_material);
298 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
305 using ATS = Kokkos::ArithTraits<scalar_type>;
307 using implATS = Kokkos::ArithTraits<impl_scalar_type>;
309 using magATS = Kokkos::ArithTraits<magnitudeType>;
335 if (!importer.is_null()) {
347 if (!importer.is_null()) {
351 ghostedMaterial = material_;
354 using execution_space =
typename Node::execution_space;
355 using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
358 auto lclMaterial = ghostedMaterial->getLocalViewDevice(Xpetra::Access::ReadOnly);
362 Kokkos::parallel_for(
"MueLu:TensorMaterialDistanceFunctor::inversion", range_type(0, lclMaterial.extent(0)), functor);
366 KOKKOS_INLINE_FUNCTION
375 auto matrix_row_material = Kokkos::subview(
material, row, Kokkos::ALL(), Kokkos::ALL());
376 auto dist = Kokkos::subview(
lcl_dist, row, Kokkos::ALL());
378 for (
size_t j = 0; j <
coords.extent(1); ++j) {
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);
385 for (
size_t j = 0; j <
coords.extent(1); ++j) {
393 auto matrix_col_material = Kokkos::subview(
material, col, Kokkos::ALL(), Kokkos::ALL());
394 auto dist = Kokkos::subview(
lcl_dist, row, Kokkos::ALL());
396 for (
size_t j = 0; j <
coords.extent(1); ++j) {
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);
403 for (
size_t j = 0; j <
coords.extent(1); ++j) {
408 return Kokkos::max(implATS::magnitude(d_row), implATS::magnitude(d_col));
415 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType>
418 DistanceFunctorType& distFunctor) {
419 using scalar_type =
Scalar;
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>;
432 auto lclA = A.getLocalMatrixDevice();
433 auto lclDiag = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
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;
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);
448 d = distFunctor.distance2(row, col);
449 d2 += implATS::one() / d;
450 haveAddedToDiag =
true;
456 lclDiag(row, 0) = !haveAddedToDiag ? implATS::squareroot(implATS::rmax()) : d2;
460 if (!importer.is_null()) {
469 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType>
472 DistanceFunctorType& distFunctor) {
473 using scalar_type =
Scalar;
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>;
487 auto lclA = A.getLocalMatrixDevice();
488 auto lclDiag = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
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;
497 impl_scalar_type mymax = implATS::zero();
500 for (local_ordinal_type colID = 0; colID < length; colID++) {
501 auto col = rowView.colidx(colID);
503 d = distFunctor.distance2(row, col);
504 d2 = implATS::one() / d;
505 if (implATS::magnitude(mymax) < -implATS::magnitude(d2))
506 mymax = -implATS::magnitude(d2);
509 lclDiag(row, 0) = mymax;
513 if (!importer.is_null()) {
532 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType, Misc::StrengthMeasure measure>
541 using diag_view_type =
typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
545 using ATS = Kokkos::ArithTraits<scalar_type>;
548 using mATS = Kokkos::ArithTraits<magnitudeType>;
561 :
A(A_.getLocalMatrixDevice())
567 auto lclDiag2d =
diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
568 diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
571 auto lclDiag2d =
diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
572 diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
576 KOKKOS_FORCEINLINE_FUNCTION
578 auto row =
A.rowConst(rlid);
579 const size_t offset =
A.graph.row_map(rlid);
581 auto clid = row.colidx(k);
585 val =
one /
dist2.distance2(rlid, clid);
591 auto aiiajj = ATS::magnitude(
diag(rlid)) * ATS::magnitude(
diag(clid));
592 auto aij2 = ATS::magnitude(val) * ATS::magnitude(val);
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,
602 auto aiiajj = ATS::magnitude(
diag(rlid)) * ATS::magnitude(
diag(clid));
603 const bool is_nonpositive = ATS::real(val) <= mATS::zero();
604 magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val);
615 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType, Misc::StrengthMeasure measure>
625 using diag_view_type =
typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
629 using ATS = Kokkos::ArithTraits<scalar_type>;
632 using mATS = Kokkos::ArithTraits<magnitudeType>;
647 :
A(A_.getLocalMatrixDevice())
655 auto lclDiag2d =
diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
656 diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
659 auto lclDiag2d =
diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
660 diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
664 KOKKOS_FORCEINLINE_FUNCTION
667 auto row =
A.rowConst(rlid);
668 const size_t offset =
A.graph.row_map(rlid);
670 auto clid = row.colidx(k);
674 if (brlid != bclid) {
675 val =
one /
dist2.distance2(brlid, bclid);
681 auto aiiajj = ATS::magnitude(
diag(brlid)) * ATS::magnitude(
diag(bclid));
682 auto aij2 = ATS::magnitude(val) * ATS::magnitude(val);
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,
692 auto aiiajj = ATS::magnitude(
diag(brlid)) * ATS::magnitude(
diag(bclid));
693 const bool is_nonpositive = ATS::real(val) <= mATS::zero();
694 magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val);
706 template <Misc::StrengthMeasure measure,
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType>
709 DistanceFunctorType& dist2_,
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_);
Drops entries the unscaled distance Laplacian.
typename matrix_type::local_matrix_type local_matrix_type
typename matrix_type::local_matrix_type local_matrix_type
local_coords_type ghostedCoords
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_)
Teuchos::RCP< coords_type > coordsMV
material_vector_type material_vector
Teuchos::RCP< material_type > materialMV
local_coords_type ghostedCoords
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
Kokkos::ArithTraits< magnitudeType > magATS
typename local_matrix_type::memory_space memory_space
Kokkos::ArithTraits< scalar_type > ATS
Kokkos::ArithTraits< scalar_type > ATS
typename Kokkos::DualView< const scalar_type *, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged >::t_dev diag_view_type
Kokkos::ArithTraits< scalar_type > ATS
typename local_matrix_type::memory_space memory_space
typename ATS::val_type impl_scalar_type
typename coords_type::dual_view_type_const::t_dev local_coords_type
typename ATS::val_type impl_scalar_type
KOKKOS_INLINE_FUNCTION magnitudeType distance2(const local_ordinal_type row, const local_ordinal_type col) const
typename local_matrix_type::value_type scalar_type
typename implATS::magnitudeType magnitudeType
KOKKOS_INLINE_FUNCTION magnitudeType distance2(const local_ordinal_type row, const local_ordinal_type col) const
Teuchos::RCP< material_type > ghostedMaterialMV
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
Teuchos::RCP< coords_type > ghostedCoordsMV
typename implATS::magnitudeType magnitudeType
local_material_type ghostedMaterial
BlockWeightedDistanceFunctor(matrix_type &A, Teuchos::RCP< coords_type > &coords_, weight_type weight_, local_ordinal_type interleaved_blocksize_)
Teuchos::RCP< coords_type > coordsMV
typename local_matrix_type::ordinal_type local_ordinal_type
typename matrix_type::local_matrix_type local_matrix_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
LocalOrdinal local_ordinal_type
Kokkos::ArithTraits< impl_scalar_type > implATS
Teuchos::RCP< coords_type > coordsMV
Teuchos::RCP< coords_type > coordsMV
Kokkos::ArithTraits< magnitudeType > magATS
Teuchos::RCP< coords_type > ghostedCoordsMV
TensorMaterialDistanceFunctor(matrix_type &A, Teuchos::RCP< coords_type > &coords_, Teuchos::RCP< material_type > &material_)
typename local_matrix_type::value_type scalar_type
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
local_coords_type ghostedCoords
local_material_type material
Kokkos::View< DecisionType *, memory_space > results_view
Kokkos::ArithTraits< magnitudeType > magATS
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Kokkos::ArithTraits< impl_scalar_type > implATS
LocalOrdinal local_ordinal_type
typename local_matrix_type::value_type scalar_type
TensorInversion(material_vector_type &material_vector_, material_matrix_type &material_matrix_)
typename ATS::val_type impl_scalar_type
local_ordinal_type interleaved_blocksize
typename local_matrix_type::value_type scalar_type
block_indices_view_type point_to_block
typename Kokkos::DualView< const scalar_type *, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged >::t_dev diag_view_type
block_indices_view_type ghosted_point_to_block
Kokkos::ArithTraits< magnitudeType > mATS
typename matrix_type::local_matrix_type local_matrix_type
typename ATS::magnitudeType magnitudeType
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_)
typename implATS::magnitudeType magnitudeType
typename local_matrix_type::value_type scalar_type
typename ATS::magnitudeType magnitudeType
ScalarMaterialDistanceFunctor(matrix_type &A, Teuchos::RCP< coords_type > &coords_, Teuchos::RCP< material_type > &material_)
LocalOrdinal local_ordinal_type
typename implATS::magnitudeType magnitudeType
Computes the weighted distance Laplacian.
Kokkos::View< const bool *, memory_space > boundary_nodes_view
typename ATS::val_type impl_scalar_type
KOKKOS_FORCEINLINE_FUNCTION magnitudeType distance2(const local_ordinal_type row, const local_ordinal_type col) const
Kokkos::ArithTraits< impl_scalar_type > implATS
KOKKOS_FORCEINLINE_FUNCTION magnitudeType distance2(const local_ordinal_type row, const local_ordinal_type col) const
Teuchos::RCP< coords_type > ghostedCoordsMV
typename implATS::magnitudeType magnitudeType
Kokkos::ArithTraits< scalar_type > ATS
Computes the weighted distance Laplacian.
Kokkos::ArithTraits< magnitudeType > mATS
Kokkos::View< impl_scalar_type **, memory_space > local_dist_type
virtual RCP< const CrsGraph > getCrsGraph() const =0
local_material_type material
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type rlid) const
DistanceFunctorType dist2
Kokkos::ArithTraits< scalar_type > ATS
Kokkos::ArithTraits< magnitudeType > magATS
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
typename ATS::val_type impl_scalar_type
Kokkos::ArithTraits< scalar_type > ATS
Computes the unscaled distance Laplacian.
LocalOrdinal local_ordinal_type
typename matrix_type::local_matrix_type local_matrix_type
local_coords_type ghostedCoords
KOKKOS_FORCEINLINE_FUNCTION void operator()(local_ordinal_type i) const
local_coords_type ghostedCoords
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
typename local_matrix_type::value_type scalar_type
Teuchos::RCP< coords_type > ghostedCoordsMV
LocalOrdinal local_ordinal_type
DistanceFunctorType dist2
Kokkos::ArithTraits< impl_scalar_type > implATS
Kokkos::ArithTraits< scalar_type > ATS
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< coords_type > ghostedCoordsMV
Kokkos::ArithTraits< impl_scalar_type > implATS
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_)
material_matrix_type material_matrix
Teuchos::RCP< coords_type > coordsMV
typename local_matrix_type::value_type scalar_type
Teuchos::RCP< diag_vec_type > diagVec
Teuchos::RCP< diag_vec_type > diagVec
typename matrix_type::local_matrix_type local_matrix_type
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
Kokkos::ArithTraits< magnitudeType > magATS