MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_UtilitiesBase_decl.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_UTILITIESBASE_DECL_HPP
11 #define MUELU_UTILITIESBASE_DECL_HPP
12 
13 #include <string>
14 
15 #include "MueLu_ConfigDefs.hpp"
16 
17 #include <Teuchos_DefaultComm.hpp>
18 #include <Teuchos_ScalarTraits.hpp>
20 
21 #if KOKKOS_VERSION >= 40799
22 #include "KokkosKernels_ArithTraits.hpp"
23 #else
24 #include "Kokkos_ArithTraits.hpp"
25 #endif
26 
31 #include <Xpetra_CrsGraph_fwd.hpp>
32 #include <Xpetra_CrsMatrix_fwd.hpp>
34 #include <Xpetra_Import_fwd.hpp>
36 #include <Xpetra_Map_fwd.hpp>
38 #include <Xpetra_Matrix_fwd.hpp>
41 #include <Xpetra_Operator_fwd.hpp>
42 #include <Xpetra_Vector_fwd.hpp>
44 
45 namespace MueLu {
46 
47 // MPI helpers
48 #define MueLu_sumAll(rcpComm, in, out) \
49  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out))
50 #define MueLu_minAll(rcpComm, in, out) \
51  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out))
52 #define MueLu_maxAll(rcpComm, in, out) \
53  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out))
54 
62 template <class Scalar,
65  class Node = DefaultNode>
66 class UtilitiesBase {
67  public:
68 #undef MUELU_UTILITIESBASE_SHORT
69 #include "MueLu_UseShortNames.hpp"
70  public:
72 
74 
81  static RCP<CrsMatrixWrap> GetThresholdedMatrix(const RCP<Matrix>& Ain, const Magnitude threshold, const bool keepDiagonal = true, const GlobalOrdinal expectedNNZperRow = -1);
82 
89  static RCP<Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>> GetThresholdedGraph(const RCP<Matrix>& A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow = -1);
90 
97  static Teuchos::ArrayRCP<Scalar> GetMatrixDiagonal_arcp(const Matrix& A);
98 
105  static RCP<Vector> GetMatrixDiagonal(const Matrix& A);
106 
113  // static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero());
114  static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps() * 100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero(), const bool doLumped = false);
115 
122  static Teuchos::RCP<Vector> GetLumpedMatrixDiagonal(Matrix const& A, const bool doReciprocal = false,
124  Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero(),
125  const bool replaceSingleEntryRowWithZero = false,
126  const bool useAverageAbsDiagVal = false);
127 
135 
137 
146 
154  static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A);
155 
164  static RCP<Vector> GetMatrixOverlappedDeletedRowsum(const Matrix& A);
165 
167  GetMatrixOverlappedAbsDeletedRowsum(const Matrix& A);
168 
176  static GlobalOrdinal CountNegativeDiagonalEntries(const Matrix& A);
177 
178  // TODO: should NOT return an Array. Definition must be changed to:
179  // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
180  // or
181  // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
182  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS);
183 
184  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector& Resid);
185 
186  static RCP<MultiVector> Residual(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS);
187 
188  static void Residual(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector& Resid);
189 
201  static Scalar PowerMethod(const Matrix& A, bool scaleByDiag = true,
202  LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123);
203 
215  static Scalar PowerMethod(const Matrix& A, const RCP<Vector>& diagInvVec,
216  LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123);
217 
218  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os);
219 
225 
231 
245 
253  static Kokkos::View<bool*, typename NO::device_type::memory_space> DetectDirichletRows_kokkos(const Matrix& A, const Magnitude& tol = Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<SC>::magnitudeType>::zero(), const bool count_twos_as_dirichlet = false);
254  static Kokkos::View<bool*, typename Kokkos::HostSpace> DetectDirichletRows_kokkos_host(const Matrix& A, const Magnitude& tol = Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<SC>::magnitudeType>::zero(), const bool count_twos_as_dirichlet = false);
255 
269 
276 
284  static void FindNonZeros(const Teuchos::ArrayRCP<const Scalar> vals,
285  Teuchos::ArrayRCP<bool> nonzeros);
286 
294  Kokkos::View<bool*, typename Node::device_type> nonzeros);
295 
305  const Teuchos::ArrayRCP<bool>& dirichletRows,
306  Teuchos::ArrayRCP<bool> dirichletCols,
307  Teuchos::ArrayRCP<bool> dirichletDomain);
308 
317  const Kokkos::View<bool*, typename Node::device_type>& dirichletRows,
318  Kokkos::View<bool*, typename Node::device_type> dirichletCols,
319  Kokkos::View<bool*, typename Node::device_type> dirichletDomain);
320 
333 
335 
336  static void ApplyRowSumCriterion(const Matrix& A,
337  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
338  Kokkos::View<bool*, typename NO::device_type::memory_space>& dirichletRows);
339 
340  static void ApplyRowSumCriterionHost(const Matrix& A,
341  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
342  Kokkos::View<bool*, Kokkos::HostSpace>& dirichletRows);
343 
344  static void ApplyRowSumCriterion(const Matrix& A,
346  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
347  Kokkos::View<bool*, typename NO::device_type::memory_space>& dirichletRows);
348 
349  static void ApplyRowSumCriterionHost(const Matrix& A,
351  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
352  Kokkos::View<bool*, Kokkos::HostSpace>& dirichletRows);
353 
365  const Teuchos::ArrayRCP<const bool>& dirichletRows);
366 
377  static Kokkos::View<bool*, typename NO::device_type> DetectDirichletCols(const Matrix& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows);
378 
384 
394  static void SetRandomSeed(const Teuchos::Comm<int>& comm);
395 
396  // Finds the OAZ Dirichlet rows for this matrix
397  // so far only used in IntrepidPCoarsenFactory
398  // TODO check whether we can use DetectDirichletRows instead
400  std::vector<LocalOrdinal>& dirichletRows, bool count_twos_as_dirichlet = false);
401 
402  // Applies Ones-and-Zeros to matrix rows
403  // Takes a vector of row indices
405  const std::vector<LocalOrdinal>& dirichletRows);
406 
407  // Applies Ones-and-Zeros to matrix rows
408  // Takes a Boolean array.
410  const Teuchos::ArrayRCP<const bool>& dirichletRows);
411 
412  static void ApplyOAZToMatrixRows(RCP<Matrix>& A, const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows);
413 
414  // Zeros out rows
415  // Takes a vector containg Dirichlet row indices
417  const std::vector<LocalOrdinal>& dirichletRows,
419 
420  // Zeros out rows
421  // Takes a Boolean ArrayRCP
423  const Teuchos::ArrayRCP<const bool>& dirichletRows,
425 
426  // Zeros out rows
427  // Takes a Boolean ArrayRCP
429  const Teuchos::ArrayRCP<const bool>& dirichletRows,
431 
432  static void ZeroDirichletRows(RCP<Matrix>& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows, SC replaceWith = Teuchos::ScalarTraits<SC>::zero());
433 
434  static void ZeroDirichletRows(RCP<MultiVector>& X, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows, SC replaceWith = Teuchos::ScalarTraits<SC>::zero());
435 
436  // Zeros out columns
437  // Takes a Boolean vector
439  const Teuchos::ArrayRCP<const bool>& dirichletCols,
441 
442  static void ZeroDirichletCols(RCP<Matrix>& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletCols, SC replaceWith = Teuchos::ScalarTraits<SC>::zero());
443 
444  // Finds the OAZ Dirichlet rows for this matrix
448 
450  // You can use this to de-normalize a tenative prolongator, for instance
452 
453  // This routine takes a BlockedMap and an Importer (assuming that the BlockedMap matches the source of the importer) and generates a BlockedMap corresponding
454  // to the Importer's target map. We assume that the targetMap is unique (which, is not a strict requirement of an Importer, but is here and no, we don't check)
455  // This is largely intended to be used in repartitioning of blocked matrices
458 
459  // Checks to see if the first chunk of the colMap is also the row map. This simiplifies a bunch of
460  // operation in coarsening
462 
467 
471 
472 }; // class UtilitiesBase
473 
474 // Useful Kokkos conversions
475 template <class View, unsigned AppendValue>
476 struct AppendTrait {
477  // static_assert(false, "Error: NOT a Kokkos::View");
478 };
479 
480 template <class MT, unsigned T>
482  // empty
483 };
484 
485 template <unsigned U, unsigned T>
486 struct CombineMemoryTraits<Kokkos::MemoryTraits<U>, T> {
487  typedef Kokkos::MemoryTraits<U | T> type;
488 };
489 
490 template <class DataType, unsigned T, class... Pack>
491 struct AppendTrait<Kokkos::View<DataType, Pack...>, T> {
492  typedef Kokkos::View<DataType, Pack...> view_type;
493  using type = Kokkos::View<DataType, typename view_type::array_layout, typename view_type::device_type, typename CombineMemoryTraits<typename view_type::memory_traits, T>::type>;
494 };
495 
497 
498 } // namespace MueLu
499 
500 #define MUELU_UTILITIESBASE_SHORT
501 #endif // MUELU_UTILITIESBASE_DECL_HPP
static void FindNonZeros(const Teuchos::ArrayRCP< const Scalar > vals, Teuchos::ArrayRCP< bool > nonzeros)
Find non-zero values in an ArrayRCP Compares the value to 2 * machine epsilon.
KOKKOS_FORCEINLINE_FUNCTION bool isDirichletRow(typename CrsMatrix::ordinal_type rowId, KokkosSparse::SparseRowViewConst< CrsMatrix > &row, const typename Kokkos::ArithTraits< typename CrsMatrix::value_type >::magnitudeType &tol, const bool count_twos_as_dirichlet)
static RCP< Vector > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, const std::vector< LocalOrdinal > &dirichletRows)
static Teuchos::RCP< Vector > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Return vector containing: max_{i=k}(-a_ik), for each for i in the matrix.
static void ZeroDirichletCols(Teuchos::RCP< Matrix > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< bool > &dirichletRows, Teuchos::ArrayRCP< bool > dirichletCols, Teuchos::ArrayRCP< bool > dirichletDomain)
Detects Dirichlet columns &amp; domains from a list of Dirichlet rows.
MueLu::DefaultLocalOrdinal LocalOrdinal
static bool MapsAreNested(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &colMap)
static RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReverseCuthillMcKee(const Matrix &Op)
static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node >> &isDirichletRow, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node >> &isDirichletCol)
static void ApplyRowSumCriterionHost(const Matrix &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, Kokkos::HostSpace > &dirichletRows)
Tpetra::KokkosClassic::DefaultNode::DefaultNodeType DefaultNode
static Scalar PowerMethod(const Matrix &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Power method.
static GlobalOrdinal CountNegativeDiagonalEntries(const Matrix &A)
Counts the number of negative diagonal entries.
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
Apply Rowsum Criterion.
static RCP< Vector > GetMatrixOverlappedDeletedRowsum(const Matrix &A)
Extract Overlapped Matrix Deleted Rowsum.
static RCP< Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > GetThresholdedGraph(const RCP< Matrix > &A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a graph.
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReplaceNonZerosWithOnes(const RCP< Matrix > &original)
Creates a copy of a matrix where the non-zero entries are replaced by ones.
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
MueLu::DefaultNode Node
static RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > CuthillMcKee(const Matrix &Op)
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
static Teuchos::ArrayRCP< const bool > DetectDirichletCols(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
Detect Dirichlet columns based on Dirichlet rows.
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
Kokkos::View< DataType, typename view_type::array_layout, typename view_type::device_type, typename CombineMemoryTraits< typename view_type::memory_traits, T >::type > type
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
static RCP< CrsMatrixWrap > GetThresholdedMatrix(const RCP< Matrix > &Ain, const Magnitude threshold, const bool keepDiagonal=true, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a matrix.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Magnitude >::zero(), bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
static RCP< Xpetra::Vector< Magnitude, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedAbsDeletedRowsum(const Matrix &A)
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal_arcp(const Matrix &A)
Extract Matrix Diagonal.
static RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > GeneratedBlockedTargetMap(const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > &sourceBlockedMap, const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer)
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
Teuchos::ScalarTraits< Scalar >::magnitudeType Magnitude
static RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps()*100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool doLumped=false)
Extract Matrix Diagonal.
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Detect Dirichlet rows (extended version)
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
static void EnforceInitialCondition(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &InitialGuess, const Magnitude &tol=Teuchos::ScalarTraits< Magnitude >::zero(), const bool count_twos_as_dirichlet=false)
Detect Dirichlet rows and copy values from RHS multivector to InitialGuess for Dirichlet rows...
Scalar SC
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)
Frobenius inner product of two matrices.
static Teuchos::RCP< Vector > GetInverse(Teuchos::RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps()*100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Return vector containing inverse of input vector.
static Kokkos::View< bool *, typename Kokkos::HostSpace > DetectDirichletRows_kokkos_host(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar >> &v, LocalOrdinal i0, LocalOrdinal i1)
Squared distance between two rows in a multivector.
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::magnitude(Teuchos::ScalarTraits< Scalar >::zero()), Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
Extract Matrix Diagonal of lumped matrix.
static Kokkos::View< bool *, typename NO::device_type::memory_space > DetectDirichletRows_kokkos(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
static RCP< Matrix > Crs2Op(RCP< CrsMatrix > Op)