Ifpack2_BlockTriDiContainer_decl.hpp

Go to the documentation of this file.

/*@HEADER
// ***********************************************************************
//
//       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2009) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
//@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);

    ~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_;
    
    // 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);

    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