Ifpack2_BlockTriDiContainer_decl.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//
// Copyright 2009 NTESS and the Ifpack2 contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef IFPACK2_BLOCKTRIDICONTAINER_DECL_HPP
#define IFPACK2_BLOCKTRIDICONTAINER_DECL_HPP


#include "Ifpack2_config.h"
#include "Ifpack2_Container.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_RowMatrix.hpp"
#include "Tpetra_BlockCrsMatrix_decl.hpp"
#include <type_traits>
#include <string>

namespace Ifpack2 {


namespace BlockTriDiContainerDetails {
struct ImplNotAvailTag {};
struct ImplSimdTag {};
struct ImplSacadoTag {};

template <typename T>
struct ImplTag { typedef ImplNotAvailTag type; };
template <>
struct ImplTag<float> { typedef ImplSimdTag type; };
template <>
struct ImplTag<double> { typedef ImplSimdTag type; };
template <>
struct ImplTag<std::complex<float> > { typedef ImplSimdTag type; };
template <>
struct ImplTag<std::complex<double> > { typedef ImplSimdTag type; };

template <typename MatrixType>
struct ImplObject;
}  // namespace BlockTriDiContainerDetails

template <typename MatrixType,
          typename ImplTagType = typename BlockTriDiContainerDetails::ImplTag<typename MatrixType::scalar_type>::type>
class BlockTriDiContainer;

template <typename MatrixType>
class BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
  : public Container<MatrixType> {

 private:
  typedef MatrixType matrix_type;

  typedef typename MatrixType::scalar_type scalar_type;
#if KOKKOS_VERSION >= 40799
  typedef typename KokkosKernels::ArithTraits<scalar_type>::magnitudeType magnitude_type;
#else
  typedef typename Kokkos::ArithTraits<scalar_type>::magnitudeType magnitude_type;
#endif
  typedef typename Container<MatrixType>::local_ordinal_type local_ordinal_type;
  typedef typename Container<MatrixType>::global_ordinal_type global_ordinal_type;
  typedef typename Container<MatrixType>::node_type node_type;

  typedef typename Container<MatrixType>::mv_type mv_type;
  typedef typename Container<MatrixType>::map_type map_type;
  typedef typename Container<MatrixType>::vector_type vector_type;
  typedef typename Container<MatrixType>::import_type import_type;

  typedef typename Container<MatrixType>::HostView host_view_type;
  typedef typename Container<MatrixType>::ConstHostView const_host_view_type;
  typedef host_view_type HostView;
  typedef const_host_view_type ConstHostView;
  // typedef Tpetra::MultiVector<local_scalar_type, local_ordinal_type, global_ordinal_type, node_type> local_mv_type;
  // typedef typename Kokkos::View<local_scalar_type**, Kokkos::HostSpace> HostViewLocal;

  typedef Tpetra::CrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> crs_matrix_type;
  typedef Tpetra::BlockCrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> block_crs_matrix_type;

  const Teuchos::Array<Teuchos::Array<local_ordinal_type> > partitions_;

  typedef typename Container<MatrixType>::row_matrix_type row_matrix_type;

  static_assert(std::is_same<MatrixType, row_matrix_type>::value,
                "Ifpack2::BlockTriDiContainer: MatrixType must be a Tpetra::RowMatrix specialization.");
 public:


  BlockTriDiContainer(const Teuchos::RCP<const row_matrix_type>& matrix,
                      const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
                      const Teuchos::RCP<const import_type>& importer,
                      bool pointIndexed);

  BlockTriDiContainer(const Teuchos::RCP<const row_matrix_type>& matrix,
                      const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
                      const int n_subparts_per_part = 1,
                      bool overlapCommAndComp       = false,
                      bool useSequentialMethod      = false,
                      const int block_size          = -1,
                      const bool explicitConversion = false);

  ~BlockTriDiContainer() override;

  struct ComputeParameters {
#if KOKKOS_VERSION >= 40799
    magnitude_type addRadiallyToDiagonal = KokkosKernels::ArithTraits<magnitude_type>::zero();
#else
    magnitude_type addRadiallyToDiagonal = Kokkos::ArithTraits<magnitude_type>::zero();
#endif
  };

  struct ApplyParameters {
    bool zeroStartingSolution = false;
#if KOKKOS_VERSION >= 40799
    scalar_type dampingFactor = KokkosKernels::ArithTraits<scalar_type>::one();
#else
    scalar_type dampingFactor            = Kokkos::ArithTraits<scalar_type>::one();
#endif
    int maxNumSweeps = 1;
#if KOKKOS_VERSION >= 40799
    magnitude_type tolerance = KokkosKernels::ArithTraits<magnitude_type>::zero();
#else
    magnitude_type tolerance             = Kokkos::ArithTraits<magnitude_type>::zero();
#endif
    int checkToleranceEvery = 1;
  };



  void setParameters(const Teuchos::ParameterList& List) override;

  void clearBlocks() override;



  void initialize() override;

  void compute() override;

  // \brief Compute <tt>Y := D^{-1} (X - R*Y)</tt>.
  void applyInverseJacobi(const mv_type& X, mv_type& Y,
                          scalar_type dampingFactor,
                          bool zeroStartingSolution = false,
                          int numSweeps             = 1) const override;

  ComputeParameters createDefaultComputeParameters() const;

  void compute(const ComputeParameters& input);

  ApplyParameters createDefaultApplyParameters() const;

  int applyInverseJacobi(const mv_type& X, mv_type& Y,
                         const ApplyParameters& input) const;

  const magnitude_type getNorms0() const;

  const magnitude_type getNormsFinal() const;

  void
  apply(const_host_view_type X,
        host_view_type Y,
        int blockIndex,
        Teuchos::ETransp mode = Teuchos::NO_TRANS,
        scalar_type alpha     = Teuchos::ScalarTraits<scalar_type>::one(),
        scalar_type beta      = Teuchos::ScalarTraits<scalar_type>::zero()) const override;

  void
  weightedApply(const_host_view_type X,
                host_view_type Y,
                const_host_view_type W,
                int blockIndex,
                Teuchos::ETransp mode = Teuchos::NO_TRANS,
                scalar_type alpha     = Teuchos::ScalarTraits<scalar_type>::one(),
                scalar_type beta      = Teuchos::ScalarTraits<scalar_type>::zero()) const override;



  std::ostream& print(std::ostream& os) const override;



  std::string description() const override;

  void
  describe(Teuchos::FancyOStream& out,
           const Teuchos::EVerbosityLevel verbLevel =
               Teuchos::Describable::verbLevel_default) const override;


  static std::string getName();

 private:
  BlockTriDiContainer(const BlockTriDiContainer<MatrixType>& rhs);

  // hide details of impl using ImplObj; finally I understand why AMB did that way.
  Teuchos::RCP<BlockTriDiContainerDetails::ImplObject<MatrixType> > impl_;
  int n_subparts_per_part_;
  int block_size_ = -1;

  // initialize distributed and local objects
  void initInternal(const Teuchos::RCP<const row_matrix_type>& matrix,
                    const Teuchos::RCP<const import_type>& importer,
                    const bool overlapCommAndComp,
                    const bool useSeqMethod,
                    const int block_size          = -1,
                    const bool explicitConversion = false);

  void clearInternal();

  // Decide whether the fused block Jacobi path can and should be used.
  bool shouldUseFusedBlockJacobi(
      const Teuchos::RCP<const row_matrix_type>& matrix,
      const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
      bool useSeqMethod);
};

template <typename MatrixType>
class BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplNotAvailTag>
  : public Container<MatrixType> {
 private:
  typedef typename MatrixType::scalar_type scalar_type;
#if KOKKOS_VERSION >= 40799
  typedef typename KokkosKernels::ArithTraits<scalar_type>::magnitudeType magnitude_type;
#else
  typedef typename Kokkos::ArithTraits<scalar_type>::magnitudeType magnitude_type;
#endif
  typedef typename Container<MatrixType>::local_ordinal_type local_ordinal_type;
  typedef typename Container<MatrixType>::global_ordinal_type global_ordinal_type;

  typedef typename Container<MatrixType>::mv_type mv_type;
  typedef typename Container<MatrixType>::import_type import_type;

  typedef typename Container<MatrixType>::HostView host_view_type;
  typedef typename Container<MatrixType>::ConstHostView const_host_view_type;
  typedef typename Container<MatrixType>::row_matrix_type row_matrix_type;

  static_assert(std::is_same<MatrixType, row_matrix_type>::value,
                "Ifpack2::BlockTriDiContainer: MatrixType must be a Tpetra::RowMatrix specialization.");

 public:
  BlockTriDiContainer(const Teuchos::RCP<const row_matrix_type>& matrix,
                      const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
                      const Teuchos::RCP<const import_type>& importer,
                      bool pointIndexed)
    : Container<MatrixType>(matrix, partitions, pointIndexed) {
    TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Error: BlockTriDiContainer is not available for this scalar_type");
  }

  void setParameters(const Teuchos::ParameterList& List) override {}
  void clearBlocks() override {}

  void initialize() override {}
  void compute() override {}
  void applyInverseJacobi(const mv_type& X, mv_type& Y,
                          scalar_type dampingFactor,
                          bool zeroStartingSolution = false,
                          int numSweeps             = 1) const override {}

  void
  apply(const_host_view_type X,
        host_view_type Y,
        int blockIndex,
        Teuchos::ETransp mode = Teuchos::NO_TRANS,
        scalar_type alpha     = Teuchos::ScalarTraits<scalar_type>::one(),
        scalar_type beta      = Teuchos::ScalarTraits<scalar_type>::zero()) const override {}

  void
  weightedApply(const_host_view_type X,
                host_view_type Y,
                const_host_view_type W,
                int blockIndex,
                Teuchos::ETransp mode = Teuchos::NO_TRANS,
                scalar_type alpha     = Teuchos::ScalarTraits<scalar_type>::one(),
                scalar_type beta      = Teuchos::ScalarTraits<scalar_type>::zero()) const override {}

  std::ostream& print(std::ostream& os) const override {
    return os << "Ifpack2::BlockTriDiContainer::ImplNotAvailTag";
  }

  std::string description() const override {
    return "Ifpack2::BlockTriDiContainer::ImplNotAvailTag";
  }

  void
  describe(Teuchos::FancyOStream& out,
           const Teuchos::EVerbosityLevel verbLevel =
               Teuchos::Describable::verbLevel_default) const override {
    out << "Ifpack2::BlockTriDiContainer::ImplNotAvailTag";
  }

  static std::string getName() {
    return "Ifpack2::BlockTriDiContainer::ImplNotAvailTag";
  }
};

}  // namespace Ifpack2

#endif  // IFPACK2_BLOCKTRIDICONTAINER_DECL_HPP