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;
  }
  
  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;
    typedef typename Kokkos::ArithTraits<scalar_type>::magnitudeType magnitude_type;
    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 {
      magnitude_type addRadiallyToDiagonal = Kokkos::ArithTraits<magnitude_type>::zero();
    };

    struct ApplyParameters {
      bool zeroStartingSolution = false;
      scalar_type dampingFactor = Kokkos::ArithTraits<scalar_type>::one();
      int maxNumSweeps = 1;
      magnitude_type tolerance = Kokkos::ArithTraits<magnitude_type>::zero();
      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();
  };
  
  template <typename MatrixType>
  class BlockTriDiContainer<MatrixType,BlockTriDiContainerDetails::ImplNotAvailTag> 
    : public Container<MatrixType> {
  private:
    typedef typename MatrixType::scalar_type scalar_type;
    typedef typename Kokkos::ArithTraits<scalar_type>::magnitudeType magnitude_type;
    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