MueLu_UtilitiesBase_decl.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//        MueLu: A package for multigrid based preconditioning
//
// Copyright 2012 NTESS and the MueLu contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef MUELU_UTILITIESBASE_DECL_HPP
#define MUELU_UTILITIESBASE_DECL_HPP

#include <string>

#include "MueLu_ConfigDefs.hpp"

#include <Teuchos_DefaultComm.hpp>
#include <Teuchos_ScalarTraits.hpp>
#include <Teuchos_ParameterList.hpp>

#include "Kokkos_ArithTraits.hpp"

#include <Xpetra_BlockedCrsMatrix_fwd.hpp>
#include <Xpetra_BlockedMap_fwd.hpp>
#include <Xpetra_BlockedVector_fwd.hpp>
#include <Xpetra_CrsGraphFactory_fwd.hpp>
#include <Xpetra_CrsGraph_fwd.hpp>
#include <Xpetra_CrsMatrix_fwd.hpp>
#include <Xpetra_CrsMatrixWrap_fwd.hpp>
#include <Xpetra_Import_fwd.hpp>
#include <Xpetra_ImportFactory_fwd.hpp>
#include <Xpetra_Map_fwd.hpp>
#include <Xpetra_MapFactory_fwd.hpp>
#include <Xpetra_Matrix_fwd.hpp>
#include <Xpetra_MultiVector_fwd.hpp>
#include <Xpetra_MultiVectorFactory_fwd.hpp>
#include <Xpetra_Operator_fwd.hpp>
#include <Xpetra_Vector_fwd.hpp>
#include <Xpetra_VectorFactory_fwd.hpp>

namespace MueLu {

// MPI helpers
#define MueLu_sumAll(rcpComm, in, out) \
  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out))
#define MueLu_minAll(rcpComm, in, out) \
  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out))
#define MueLu_maxAll(rcpComm, in, out) \
  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out))

template <class Scalar,
          class LocalOrdinal  = DefaultLocalOrdinal,
          class GlobalOrdinal = DefaultGlobalOrdinal,
          class Node          = DefaultNode>
class UtilitiesBase {
 public:
#undef MUELU_UTILITIESBASE_SHORT
#include "MueLu_UseShortNames.hpp"
 public:
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;

  static RCP<Matrix> Crs2Op(RCP<CrsMatrix> Op);

  static RCP<CrsMatrixWrap> GetThresholdedMatrix(const RCP<Matrix>& Ain, const Magnitude threshold, const bool keepDiagonal = true, const GlobalOrdinal expectedNNZperRow = -1);

  static RCP<Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>> GetThresholdedGraph(const RCP<Matrix>& A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow = -1);

  static Teuchos::ArrayRCP<Scalar> GetMatrixDiagonal_arcp(const Matrix& A);

  static RCP<Vector> GetMatrixDiagonal(const Matrix& A);

  // static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero());
  static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps() * 100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero(), const bool doLumped = false);

  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);

  static Teuchos::RCP<Xpetra::Vector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>> GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A);

  static Teuchos::RCP<Xpetra::Vector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>> GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, const Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& BlockNumber);

  static Teuchos::RCP<Vector> GetInverse(Teuchos::RCP<const Vector> v, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps() * 100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero());

  static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A);

  static RCP<Vector> GetMatrixOverlappedDeletedRowsum(const Matrix& A);

  static RCP<Xpetra::Vector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>>
  GetMatrixOverlappedAbsDeletedRowsum(const Matrix& A);

  // TODO: should NOT return an Array. Definition must be changed to:
  // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
  // or
  // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS);

  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector& Resid);

  static RCP<MultiVector> Residual(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS);

  static void Residual(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector& Resid);

  static Scalar PowerMethod(const Matrix& A, bool scaleByDiag = true,
                            LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123);

  static Scalar PowerMethod(const Matrix& A, const RCP<Vector>& diagInvVec,
                            LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123);

  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os);

  static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>>& v, LocalOrdinal i0, LocalOrdinal i1);

  static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const Teuchos::ArrayView<double>& weight, const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>>& v, LocalOrdinal i0, LocalOrdinal i1);

  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);

  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);
  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::ArrayRCP<const bool> DetectDirichletRowsExt(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, bool& bHasZeroDiagonal, const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero());

  static void FindNonZeros(const Teuchos::ArrayRCP<const Scalar> vals,
                           Teuchos::ArrayRCP<bool> nonzeros);

  static void FindNonZeros(const typename Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::dual_view_type::t_dev_const_um vals,
                           Kokkos::View<bool*, typename Node::device_type> nonzeros);

  static void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
                                            const Teuchos::ArrayRCP<bool>& dirichletRows,
                                            Teuchos::ArrayRCP<bool> dirichletCols,
                                            Teuchos::ArrayRCP<bool> dirichletDomain);

  static void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
                                            const Kokkos::View<bool*, typename Node::device_type>& dirichletRows,
                                            Kokkos::View<bool*, typename Node::device_type> dirichletCols,
                                            Kokkos::View<bool*, typename Node::device_type> dirichletDomain);

  static void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, const Magnitude rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows);

  static void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, const Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& BlockNumber, const Magnitude rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows);

  static void ApplyRowSumCriterion(const Matrix& A,
                                   const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
                                   Kokkos::View<bool*, typename NO::device_type::memory_space>& dirichletRows);

  static void ApplyRowSumCriterionHost(const Matrix& A,
                                       const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
                                       Kokkos::View<bool*, Kokkos::HostSpace>& dirichletRows);

  static void ApplyRowSumCriterion(const Matrix& A,
                                   const Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& BlockNumber,
                                   const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
                                   Kokkos::View<bool*, typename NO::device_type::memory_space>& dirichletRows);

  static void ApplyRowSumCriterionHost(const Matrix& A,
                                       const Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& BlockNumber,
                                       const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
                                       Kokkos::View<bool*, Kokkos::HostSpace>& dirichletRows);

  static Teuchos::ArrayRCP<const bool> DetectDirichletCols(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
                                                           const Teuchos::ArrayRCP<const bool>& dirichletRows);

  static Kokkos::View<bool*, typename NO::device_type> DetectDirichletCols(const Matrix& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows);

  static Scalar Frobenius(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& B);

  static void SetRandomSeed(const Teuchos::Comm<int>& comm);

  // Finds the OAZ Dirichlet rows for this matrix
  // so far only used in IntrepidPCoarsenFactory
  // TODO check whether we can use DetectDirichletRows instead
  static void FindDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
                                std::vector<LocalOrdinal>& dirichletRows, bool count_twos_as_dirichlet = false);

  // Applies Ones-and-Zeros to matrix rows
  // Takes a vector of row indices
  static void ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
                                   const std::vector<LocalOrdinal>& dirichletRows);

  // Applies Ones-and-Zeros to matrix rows
  // Takes a Boolean array.
  static void ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
                                   const Teuchos::ArrayRCP<const bool>& dirichletRows);

  static void ApplyOAZToMatrixRows(RCP<Matrix>& A, const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows);

  // Zeros out rows
  // Takes a vector containg Dirichlet row indices
  static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
                                const std::vector<LocalOrdinal>& dirichletRows,
                                Scalar replaceWith = Teuchos::ScalarTraits<Scalar>::zero());

  // Zeros out rows
  // Takes a Boolean ArrayRCP
  static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
                                const Teuchos::ArrayRCP<const bool>& dirichletRows,
                                Scalar replaceWith = Teuchos::ScalarTraits<Scalar>::zero());

  // Zeros out rows
  // Takes a Boolean ArrayRCP
  static void ZeroDirichletRows(Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& X,
                                const Teuchos::ArrayRCP<const bool>& dirichletRows,
                                Scalar replaceWith = Teuchos::ScalarTraits<Scalar>::zero());

  static void ZeroDirichletRows(RCP<Matrix>& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows, SC replaceWith = Teuchos::ScalarTraits<SC>::zero());

  static void ZeroDirichletRows(RCP<MultiVector>& X, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows, SC replaceWith = Teuchos::ScalarTraits<SC>::zero());

  // Zeros out columns
  // Takes a Boolean vector
  static void ZeroDirichletCols(Teuchos::RCP<Matrix>& A,
                                const Teuchos::ArrayRCP<const bool>& dirichletCols,
                                Scalar replaceWith = Teuchos::ScalarTraits<Scalar>::zero());

  static void ZeroDirichletCols(RCP<Matrix>& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletCols, SC replaceWith = Teuchos::ScalarTraits<SC>::zero());

  // Finds the OAZ Dirichlet rows for this matrix
  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);

  // You can use this to de-normalize a tenative prolongator, for instance
  static RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> ReplaceNonZerosWithOnes(const RCP<Matrix>& original);

  // This routine takes a BlockedMap and an Importer (assuming that the BlockedMap matches the source of the importer) and generates a BlockedMap corresponding
  // 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)
  // This is largely intended to be used in repartitioning of blocked matrices
  static RCP<const Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>> GeneratedBlockedTargetMap(const Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>& sourceBlockedMap,
                                                                                                    const Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node>& Importer);

  // Checks to see if the first chunk of the colMap is also the row map.  This simiplifies a bunch of
  // operation in coarsening
  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>> CuthillMcKee(const Matrix& Op);

  static RCP<Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>> ReverseCuthillMcKee(const Matrix& Op);

};  // class UtilitiesBase

// Useful Kokkos conversions
template <class View, unsigned AppendValue>
struct AppendTrait {
  // static_assert(false, "Error: NOT a Kokkos::View");
};

template <class MT, unsigned T>
struct CombineMemoryTraits {
  // empty
};

template <unsigned U, unsigned T>
struct CombineMemoryTraits<Kokkos::MemoryTraits<U>, T> {
  typedef Kokkos::MemoryTraits<U | T> type;
};

template <class DataType, unsigned T, class... Pack>
struct AppendTrait<Kokkos::View<DataType, Pack...>, T> {
  typedef Kokkos::View<DataType, Pack...> view_type;
  using type = Kokkos::View<DataType, typename view_type::array_layout, typename view_type::device_type, typename CombineMemoryTraits<typename view_type::memory_traits, T>::type>;
};


}  // namespace MueLu

#define MUELU_UTILITIESBASE_SHORT
#endif  // MUELU_UTILITIESBASE_DECL_HPP