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) {
80 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
87 using ATS = Kokkos::ArithTraits<scalar_type>;
89 using implATS = Kokkos::ArithTraits<impl_scalar_type>;
91 using magATS = Kokkos::ArithTraits<magnitudeType>;
114 if (!importer.is_null()) {
133 KOKKOS_INLINE_FUNCTION
143 for (
size_t j = 0; j <
coords.extent(1); ++j) {
151 return Kokkos::max(d_row, d_col);
155 template <
class local_ordinal_type,
class material_vector_type,
class material_matrix_type>
162 TensorInversion(material_vector_type& material_vector_, material_matrix_type& material_matrix_)
166 KOKKOS_FORCEINLINE_FUNCTION
173 auto matrix_material = Kokkos::subview(
material_matrix, i, Kokkos::ALL(), Kokkos::ALL());
174 KokkosBatched::SerialLU<KokkosBatched::Algo::LU::Unblocked>::invoke(matrix_material);
178 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
185 using ATS = Kokkos::ArithTraits<scalar_type>;
187 using implATS = Kokkos::ArithTraits<impl_scalar_type>;
189 using magATS = Kokkos::ArithTraits<magnitudeType>;
215 if (!importer.is_null()) {
227 if (!importer.is_null()) {
231 ghostedMaterial = material_;
234 using execution_space =
typename Node::execution_space;
235 using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
238 auto lclMaterial = ghostedMaterial->getLocalViewDevice(Xpetra::Access::ReadOnly);
242 Kokkos::parallel_for(
"MueLu:TensorMaterialDistanceFunctor::inversion", range_type(0, lclMaterial.extent(0)), functor);
246 KOKKOS_INLINE_FUNCTION
255 auto matrix_row_material = Kokkos::subview(
material, row, Kokkos::ALL(), Kokkos::ALL());
256 auto dist = Kokkos::subview(
lcl_dist, row, Kokkos::ALL());
258 for (
size_t j = 0; j <
coords.extent(1); ++j) {
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);
265 for (
size_t j = 0; j <
coords.extent(1); ++j) {
273 auto matrix_col_material = Kokkos::subview(
material, col, Kokkos::ALL(), Kokkos::ALL());
274 auto dist = Kokkos::subview(
lcl_dist, row, Kokkos::ALL());
276 for (
size_t j = 0; j <
coords.extent(1); ++j) {
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);
283 for (
size_t j = 0; j <
coords.extent(1); ++j) {
288 return Kokkos::max(implATS::magnitude(d_row), implATS::magnitude(d_col));
295 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType>
298 DistanceFunctorType& distFunctor) {
299 using scalar_type =
Scalar;
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>;
312 auto lclA = A.getLocalMatrixDevice();
313 auto lclDiag = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
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;
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);
328 d = distFunctor.distance2(row, col);
329 d2 += implATS::one() / d;
330 haveAddedToDiag =
true;
336 lclDiag(row, 0) = !haveAddedToDiag ? implATS::squareroot(implATS::rmax()) : d2;
340 if (!importer.is_null()) {
349 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType>
352 DistanceFunctorType& distFunctor) {
353 using scalar_type =
Scalar;
356 using node_type =
Node;
357 using ATS = Kokkos::ArithTraits<scalar_type>;
358 using impl_scalar_type =
typename ATS::val_type;
359 using implATS = Kokkos::ArithTraits<impl_scalar_type>;
360 using magnitudeType =
typename implATS::magnitudeType;
361 using execution_space =
typename Node::execution_space;
362 using range_type = Kokkos::RangePolicy<LocalOrdinal, execution_space>;
363 using mATS = Kokkos::ArithTraits<magnitudeType>;
367 auto lclA = A.getLocalMatrixDevice();
368 auto lclDiag = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
370 Kokkos::parallel_for(
371 "MueLu:CoalesceDropF:Build:scalar_filter:laplacian_diag",
372 range_type(0, lclA.numRows()),
373 KOKKOS_LAMBDA(
const local_ordinal_type& row) {
374 auto rowView = lclA.rowConst(row);
375 auto length = rowView.length;
377 impl_scalar_type mymax = implATS::zero();
380 for (local_ordinal_type colID = 0; colID < length; colID++) {
381 auto col = rowView.colidx(colID);
383 d = distFunctor.distance2(row, col);
384 d2 = implATS::one() / d;
385 if (implATS::magnitude(mymax) < -implATS::magnitude(d2))
386 mymax = -implATS::magnitude(d2);
389 lclDiag(row, 0) = mymax;
393 if (!importer.is_null()) {
412 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType, Misc::StrengthMeasure measure>
421 using diag_view_type =
typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
425 using ATS = Kokkos::ArithTraits<scalar_type>;
428 using mATS = Kokkos::ArithTraits<magnitudeType>;
441 :
A(A_.getLocalMatrixDevice())
447 auto lclDiag2d =
diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
448 diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
451 auto lclDiag2d =
diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
452 diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
456 KOKKOS_FORCEINLINE_FUNCTION
458 auto row =
A.rowConst(rlid);
459 const size_t offset =
A.graph.row_map(rlid);
461 auto clid = row.colidx(k);
465 val =
one /
dist2.distance2(rlid, clid);
471 auto aiiajj = ATS::magnitude(
diag(rlid)) * ATS::magnitude(
diag(clid));
472 auto aij2 = ATS::magnitude(val) * ATS::magnitude(val);
477 auto neg_aij = -ATS::real(val);
478 auto max_neg_aik =
eps * ATS::real(
diag(rlid));
479 results(offset + k) = Kokkos::max((neg_aij <= max_neg_aik) ?
DROP :
KEEP,
482 auto aiiajj = ATS::magnitude(
diag(rlid)) * ATS::magnitude(
diag(clid));
483 const bool is_nonpositive = ATS::real(val) <= mATS::zero();
484 magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val);
495 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType, Misc::StrengthMeasure measure>
505 using diag_view_type =
typename Kokkos::DualView<const scalar_type*, Kokkos::LayoutStride, typename Node::device_type, Kokkos::MemoryUnmanaged>::t_dev;
509 using ATS = Kokkos::ArithTraits<scalar_type>;
512 using mATS = Kokkos::ArithTraits<magnitudeType>;
527 :
A(A_.getLocalMatrixDevice())
535 auto lclDiag2d =
diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
536 diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
539 auto lclDiag2d =
diagVec->getLocalViewDevice(Xpetra::Access::ReadOnly);
540 diag = Kokkos::subview(lclDiag2d, Kokkos::ALL(), 0);
544 KOKKOS_FORCEINLINE_FUNCTION
547 auto row =
A.rowConst(rlid);
548 const size_t offset =
A.graph.row_map(rlid);
550 auto clid = row.colidx(k);
554 if (brlid != bclid) {
555 val =
one /
dist2.distance2(brlid, bclid);
561 auto aiiajj = ATS::magnitude(
diag(brlid)) * ATS::magnitude(
diag(bclid));
562 auto aij2 = ATS::magnitude(val) * ATS::magnitude(val);
567 auto neg_aij = -ATS::real(val);
568 auto max_neg_aik =
eps * ATS::real(
diag(brlid));
569 results(offset + k) = Kokkos::max((neg_aij <= max_neg_aik) ?
DROP :
KEEP,
572 auto aiiajj = ATS::magnitude(
diag(brlid)) * ATS::magnitude(
diag(bclid));
573 const bool is_nonpositive = ATS::real(val) <= mATS::zero();
574 magnitudeType aij2 = ATS::magnitude(val) * ATS::magnitude(val);
586 template <Misc::StrengthMeasure measure,
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType>
589 DistanceFunctorType& dist2_,
595 template <Misc::StrengthMeasure measure,
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class DistanceFunctorType>
599 DistanceFunctorType& dist2_,
603 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
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
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
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
Kokkos::ArithTraits< impl_scalar_type > implATS
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_material_type material
Kokkos::View< DecisionType *, memory_space > results_view
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 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 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
Kokkos::View< const bool *, memory_space > boundary_nodes_view
typename ATS::val_type impl_scalar_type
Kokkos::ArithTraits< impl_scalar_type > implATS
KOKKOS_FORCEINLINE_FUNCTION magnitudeType distance2(const local_ordinal_type row, const local_ordinal_type col) const
typename implATS::magnitudeType magnitudeType
Kokkos::ArithTraits< scalar_type > ATS
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
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 local_matrix_type::ordinal_type local_ordinal_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
Teuchos::RCP< coords_type > ghostedCoordsMV
LocalOrdinal local_ordinal_type
DistanceFunctorType dist2
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
Teuchos::RCP< diag_vec_type > diagVec
Teuchos::RCP< diag_vec_type > diagVec
typename matrix_type::local_matrix_type local_matrix_type
Kokkos::ArithTraits< magnitudeType > magATS