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 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_UTILITIESBASE_DECL_HPP
47 #define MUELU_UTILITIESBASE_DECL_HPP
48 
49 #include <string>
50 
51 #include "MueLu_ConfigDefs.hpp"
52 
53 #include <Teuchos_DefaultComm.hpp>
54 #include <Teuchos_ScalarTraits.hpp>
56 
57 #include "Kokkos_ArithTraits.hpp"
58 
63 #include <Xpetra_CrsGraph_fwd.hpp>
64 #include <Xpetra_CrsMatrix_fwd.hpp>
66 #include <Xpetra_Import_fwd.hpp>
68 #include <Xpetra_Map_fwd.hpp>
70 #include <Xpetra_Matrix_fwd.hpp>
73 #include <Xpetra_Operator_fwd.hpp>
74 #include <Xpetra_Vector_fwd.hpp>
76 
77 namespace MueLu {
78 
79 // MPI helpers
80 #define MueLu_sumAll(rcpComm, in, out) \
81  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out))
82 #define MueLu_minAll(rcpComm, in, out) \
83  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out))
84 #define MueLu_maxAll(rcpComm, in, out) \
85  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out))
86 
94 template <class Scalar,
97  class Node = DefaultNode>
98 class UtilitiesBase {
99  public:
100 #undef MUELU_UTILITIESBASE_SHORT
101 #include "MueLu_UseShortNames.hpp"
102  public:
104 
105  static RCP<Matrix> Crs2Op(RCP<CrsMatrix> Op);
106 
113  static RCP<CrsMatrixWrap> GetThresholdedMatrix(const RCP<Matrix>& Ain, const Scalar threshold, const bool keepDiagonal = true, const GlobalOrdinal expectedNNZperRow = -1);
114 
121  static RCP<Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>> GetThresholdedGraph(const RCP<Matrix>& A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow = -1);
122 
129  static Teuchos::ArrayRCP<Scalar> GetMatrixDiagonal_arcp(const Matrix& A);
130 
137  static RCP<Vector> GetMatrixDiagonal(const Matrix& A);
138 
145  // static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero());
146  static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps() * 100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero(), const bool doLumped = false);
147 
154  static Teuchos::RCP<Vector> GetLumpedMatrixDiagonal(Matrix const& A, const bool doReciprocal = false,
156  Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero(),
157  const bool replaceSingleEntryRowWithZero = false,
158  const bool useAverageAbsDiagVal = false);
159 
167 
169 
178 
186  static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A);
187 
196  static RCP<Vector> GetMatrixOverlappedDeletedRowsum(const Matrix& A);
197 
199  GetMatrixOverlappedAbsDeletedRowsum(const Matrix& A);
200 
201  // TODO: should NOT return an Array. Definition must be changed to:
202  // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
203  // or
204  // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
205  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS);
206 
207  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector& Resid);
208 
209  static RCP<MultiVector> Residual(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS);
210 
211  static void Residual(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector& Resid);
212 
224  static Scalar PowerMethod(const Matrix& A, bool scaleByDiag = true,
225  LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123);
226 
238  static Scalar PowerMethod(const Matrix& A, const RCP<Vector>& diagInvVec,
239  LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123);
240 
241  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os);
242 
248 
254 
268 
276  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);
277  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);
278 
292 
300  static void FindNonZeros(const Teuchos::ArrayRCP<const Scalar> vals,
301  Teuchos::ArrayRCP<bool> nonzeros);
302 
310  Kokkos::View<bool*, typename Node::device_type> nonzeros);
311 
321  const Teuchos::ArrayRCP<bool>& dirichletRows,
322  Teuchos::ArrayRCP<bool> dirichletCols,
323  Teuchos::ArrayRCP<bool> dirichletDomain);
324 
333  const Kokkos::View<bool*, typename Node::device_type>& dirichletRows,
334  Kokkos::View<bool*, typename Node::device_type> dirichletCols,
335  Kokkos::View<bool*, typename Node::device_type> dirichletDomain);
336 
349 
351 
352  static void ApplyRowSumCriterion(const Matrix& A,
353  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
354  Kokkos::View<bool*, typename NO::device_type::memory_space>& dirichletRows);
355 
356  static void ApplyRowSumCriterionHost(const Matrix& A,
357  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
358  Kokkos::View<bool*, Kokkos::HostSpace>& dirichletRows);
359 
360  static void ApplyRowSumCriterion(const Matrix& A,
362  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
363  Kokkos::View<bool*, typename NO::device_type::memory_space>& dirichletRows);
364 
365  static void ApplyRowSumCriterionHost(const Matrix& A,
367  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
368  Kokkos::View<bool*, Kokkos::HostSpace>& dirichletRows);
369 
381  const Teuchos::ArrayRCP<const bool>& dirichletRows);
382 
393  static Kokkos::View<bool*, typename NO::device_type> DetectDirichletCols(const Matrix& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows);
394 
400 
410  static void SetRandomSeed(const Teuchos::Comm<int>& comm);
411 
412  // Finds the OAZ Dirichlet rows for this matrix
413  // so far only used in IntrepidPCoarsenFactory
414  // TODO check whether we can use DetectDirichletRows instead
416  std::vector<LocalOrdinal>& dirichletRows, bool count_twos_as_dirichlet = false);
417 
418  // Applies Ones-and-Zeros to matrix rows
419  // Takes a vector of row indices
421  const std::vector<LocalOrdinal>& dirichletRows);
422 
423  // Applies Ones-and-Zeros to matrix rows
424  // Takes a Boolean array.
426  const Teuchos::ArrayRCP<const bool>& dirichletRows);
427 
428  static void ApplyOAZToMatrixRows(RCP<Matrix>& A, const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows);
429 
430  // Zeros out rows
431  // Takes a vector containg Dirichlet row indices
433  const std::vector<LocalOrdinal>& dirichletRows,
435 
436  // Zeros out rows
437  // Takes a Boolean ArrayRCP
439  const Teuchos::ArrayRCP<const bool>& dirichletRows,
441 
442  // Zeros out rows
443  // Takes a Boolean ArrayRCP
445  const Teuchos::ArrayRCP<const bool>& dirichletRows,
447 
448  static void ZeroDirichletRows(RCP<Matrix>& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows, SC replaceWith = Teuchos::ScalarTraits<SC>::zero());
449 
450  static void ZeroDirichletRows(RCP<MultiVector>& X, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows, SC replaceWith = Teuchos::ScalarTraits<SC>::zero());
451 
452  // Zeros out columns
453  // Takes a Boolean vector
455  const Teuchos::ArrayRCP<const bool>& dirichletCols,
457 
458  static void ZeroDirichletCols(RCP<Matrix>& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletCols, SC replaceWith = Teuchos::ScalarTraits<SC>::zero());
459 
460  // Finds the OAZ Dirichlet rows for this matrix
464 
466  // You can use this to de-normalize a tenative prolongator, for instance
468 
469  // This routine takes a BlockedMap and an Importer (assuming that the BlockedMap matches the source of the importer) and generates a BlockedMap corresponding
470  // 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)
471  // This is largely intended to be used in repartitioning of blocked matrices
474 
475  // Checks to see if the first chunk of the colMap is also the row map. This simiplifies a bunch of
476  // operation in coarsening
478 
483 
487 
488 }; // class UtilitiesBase
489 
490 // Useful Kokkos conversions
491 template <class View, unsigned AppendValue>
492 struct AppendTrait {
493  // static_assert(false, "Error: NOT a Kokkos::View");
494 };
495 
496 template <class MT, unsigned T>
498  // empty
499 };
500 
501 template <unsigned U, unsigned T>
502 struct CombineMemoryTraits<Kokkos::MemoryTraits<U>, T> {
503  typedef Kokkos::MemoryTraits<U | T> type;
504 };
505 
506 template <class DataType, unsigned T, class... Pack>
507 struct AppendTrait<Kokkos::View<DataType, Pack...>, T> {
508  typedef Kokkos::View<DataType, Pack...> view_type;
509  using type = Kokkos::View<DataType, typename view_type::array_layout, typename view_type::device_type, typename CombineMemoryTraits<typename view_type::memory_traits, T>::type>;
510 };
511 
513 
514 } // namespace MueLu
515 
516 #define MUELU_UTILITIESBASE_SHORT
517 #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.
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 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.
static RCP< CrsMatrixWrap > GetThresholdedMatrix(const RCP< Matrix > &Ain, const Scalar threshold, const bool keepDiagonal=true, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a matrix.
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 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 Teuchos::ArrayRCP< Magnitude > 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 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 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)
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)