Tpetra_BlockCrsMatrix_decl.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//          Tpetra: Templated Linear Algebra Services Package
//                 Copyright (2008) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// 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.
//
// ************************************************************************
// @HEADER
 
#ifndef TPETRA_BLOCKCRSMATRIX_DECL_HPP
#define TPETRA_BLOCKCRSMATRIX_DECL_HPP
 
 
#include "Tpetra_CrsGraph.hpp"
#include "Tpetra_RowMatrix.hpp"
#include "Tpetra_BlockMultiVector_decl.hpp"
#include "Tpetra_CrsMatrix_decl.hpp"
 
namespace Tpetra {
 
template<class Scalar,
         class LO,
         class GO,
         class Node>
class BlockCrsMatrix :
  virtual public ::Tpetra::RowMatrix<Scalar, LO, GO, Node>,
  virtual public ::Tpetra::DistObject<char, LO, GO, Node>
{
private:
  using dist_object_type = ::Tpetra::DistObject<char, LO, GO, Node>;
  using BMV = BlockMultiVector<Scalar, LO, GO, Node>;
  using STS = Teuchos::ScalarTraits<Scalar>;
 
protected:
  typedef char packet_type;
 
public:
 
 
  using scalar_type = Scalar;
 
  using impl_scalar_type = typename BMV::impl_scalar_type;
 
  typedef LO local_ordinal_type;
  typedef GO global_ordinal_type;
  typedef Node node_type;
 
  typedef typename Node::device_type device_type;
  typedef typename device_type::execution_space execution_space;
  typedef typename device_type::memory_space memory_space;
 
  typedef ::Tpetra::Map<LO, GO, node_type> map_type;
  typedef ::Tpetra::MultiVector<Scalar, LO, GO, node_type> mv_type;
  typedef ::Tpetra::CrsGraph<LO, GO, node_type> crs_graph_type;
 
  typedef Kokkos::View<impl_scalar_type**,
                       Kokkos::LayoutRight,
                       device_type,
                       Kokkos::MemoryTraits<Kokkos::Unmanaged> >
          little_block_type;
  typedef typename little_block_type::HostMirror little_block_host_type;
 
  typedef Kokkos::View<const impl_scalar_type**,
                       Kokkos::LayoutRight,
                       device_type,
                       Kokkos::MemoryTraits<Kokkos::Unmanaged> >
          const_little_block_type;
  typedef typename BMV::little_vec_type little_vec_type;
  typedef typename BMV::little_host_vec_type little_host_vec_type;
 
  typedef typename BMV::const_little_vec_type const_little_vec_type;
  typedef typename BMV::const_little_host_vec_type const_host_little_vec_type;
 
  using row_matrix_type = RowMatrix<Scalar, LO, GO, node_type>;
  using local_inds_device_view_type =
        typename row_matrix_type::local_inds_device_view_type;
  using local_inds_host_view_type =
        typename row_matrix_type::local_inds_host_view_type;
  using nonconst_local_inds_host_view_type =
        typename row_matrix_type::nonconst_local_inds_host_view_type;
 
  using global_inds_device_view_type =
        typename row_matrix_type::global_inds_device_view_type;
  using global_inds_host_view_type =
        typename row_matrix_type::global_inds_host_view_type;
  using nonconst_global_inds_host_view_type =
        typename row_matrix_type::nonconst_global_inds_host_view_type;
 
  using values_device_view_type =
        typename row_matrix_type::values_device_view_type;
  using values_host_view_type =
        typename row_matrix_type::values_host_view_type;
  using nonconst_values_host_view_type =
        typename row_matrix_type::nonconst_values_host_view_type;
 
 
 
  BlockCrsMatrix ();
 
  BlockCrsMatrix (const crs_graph_type& graph, const LO blockSize);
 
  BlockCrsMatrix (const crs_graph_type& graph,
                  const map_type& domainPointMap,
                  const map_type& rangePointMap,
                  const LO blockSize);
 
  virtual ~BlockCrsMatrix () {}
 
 
 
  Teuchos::RCP<const map_type> getDomainMap () const override;
 
  Teuchos::RCP<const map_type> getRangeMap () const override;
 
  Teuchos::RCP<const map_type> getRowMap () const override;
 
  Teuchos::RCP<const map_type> getColMap () const override;
 
  global_size_t getGlobalNumRows() const override;
 
  size_t getNodeNumRows() const override;
 
  size_t getNodeMaxNumRowEntries() const override;
 
  void
  apply (const mv_type& X,
         mv_type& Y,
         Teuchos::ETransp mode = Teuchos::NO_TRANS,
         Scalar alpha = Teuchos::ScalarTraits<Scalar>::one (),
         Scalar beta = Teuchos::ScalarTraits<Scalar>::zero ()) const override;
 
  bool hasTransposeApply () const override {
    // FIXME (mfh 04 May 2014) Transpose and conjugate transpose modes
    // are not implemented yet.  Fill in applyBlockTrans() to fix this.
    return false;
  }
 
  void setAllToScalar (const Scalar& alpha);
 
 
 
  std::string description () const override;
 
  void
  describe (Teuchos::FancyOStream& out,
            const Teuchos::EVerbosityLevel verbLevel) const override;
 
 
 
  LO getBlockSize () const { return blockSize_; }
 
  virtual Teuchos::RCP<const ::Tpetra::RowGraph<LO,GO,Node> > getGraph () const override;
 
  const crs_graph_type & getCrsGraph () const { return graph_; }
 
  void
  applyBlock (const BlockMultiVector<Scalar, LO, GO, Node>& X,
              BlockMultiVector<Scalar, LO, GO, Node>& Y,
              Teuchos::ETransp mode = Teuchos::NO_TRANS,
              const Scalar alpha = Teuchos::ScalarTraits<Scalar>::one (),
              const Scalar beta = Teuchos::ScalarTraits<Scalar>::zero ());
 
  LO
  replaceLocalValues (const LO localRowInd,
                      const LO colInds[],
                      const Scalar vals[],
                      const LO numColInds) const;
 
  LO
  sumIntoLocalValues (const LO localRowInd,
                      const LO colInds[],
                      const Scalar vals[],
                      const LO numColInds) const;
 
 
#ifdef TPETRA_ENABLE_DEPRECATED_CODE
  LO
  getLocalRowView (const LO localRowInd,
                   const LO*& colInds,
                   Scalar*& vals,
                   LO& numInds) const;
 
 
  void
  getLocalRowView (LO LocalRow,
                   Teuchos::ArrayView<const LO> &indices,
                   Teuchos::ArrayView<const Scalar> &values) const override;
#endif // TPETRA_ENABLE_DEPRECATED_CODE
  void
  getLocalRowView (LO LocalRow,
                   local_inds_host_view_type &indices,
                   values_host_view_type &values) const override;
 
  void
  getLocalRowViewNonConst (LO LocalRow,
                           local_inds_host_view_type &indices,
                           nonconst_values_host_view_type &values) const;
  
  virtual void
  getLocalRowCopy (LO LocalRow,
                   nonconst_local_inds_host_view_type &Indices,
                   nonconst_values_host_view_type &Values,
                   size_t& NumEntries) const override;
#ifdef TPETRA_ENABLE_DEPRECATED_CODE
  virtual void
  getLocalRowCopy (LO LocalRow,
                   const Teuchos::ArrayView<LO> &Indices,
                   const Teuchos::ArrayView<Scalar> &Values,
                   size_t &NumEntries) const override;
#endif
 
#ifdef TPETRA_ENABLE_DEPRECATED_CODE
  little_block_type
  getLocalBlock (const LO localRowInd, const LO localColInd) const;
#endif
 
  little_block_type
  getLocalBlockDeviceNonConst (const LO localRowInd, const LO localColInd) const;
 
  little_block_host_type
  getLocalBlockHostNonConst (const LO localRowInd, const LO localColInd) const;
 
 
  LO
  getLocalRowOffsets (const LO localRowInd,
                      ptrdiff_t offsets[],
                      const LO colInds[],
                      const LO numColInds) const;
 
  LO
  replaceLocalValuesByOffsets (const LO localRowInd,
                               const ptrdiff_t offsets[],
                               const Scalar vals[],
                               const LO numOffsets) const;
 
  LO
  absMaxLocalValuesByOffsets (const LO localRowInd,
                              const ptrdiff_t offsets[],
                              const Scalar vals[],
                              const LO numOffsets) const;
 
  LO
  sumIntoLocalValuesByOffsets (const LO localRowInd,
                               const ptrdiff_t offsets[],
                               const Scalar vals[],
                               const LO numOffsets) const;
 
  size_t getNumEntriesInLocalRow (const LO localRowInd) const override;
 
  bool localError () const {
    return *localError_;
  }
 
  std::string errorMessages () const {
    return (*errs_).is_null () ? std::string ("") : (*errs_)->str ();
  }
 
  void
  getLocalDiagOffsets (const Kokkos::View<size_t*, device_type,
                         Kokkos::MemoryUnmanaged>& offsets) const;
 
 
  void
  getLocalDiagCopy (const Kokkos::View<impl_scalar_type***, device_type,
                                       Kokkos::MemoryUnmanaged>& diag,
                    const Kokkos::View<const size_t*, device_type,
                                       Kokkos::MemoryUnmanaged>& offsets) const;
 
#ifdef TPETRA_ENABLE_DEPRECATED_CODE
  void
  getLocalDiagCopy (const Kokkos::View<impl_scalar_type***, device_type,
                                       Kokkos::MemoryUnmanaged>& diag,
                    const Teuchos::ArrayView<const size_t>& offsets) const;
#endif
 
protected:
  LO
  absMaxLocalValues (const LO localRowInd,
                     const LO colInds[],
                     const Scalar vals[],
                     const LO numColInds) const;
 
 
 
  using buffer_device_type = typename DistObject<Scalar, LO, GO,
                                                 Node>::buffer_device_type;
 
  virtual bool checkSizes (const ::Tpetra::SrcDistObject& source) override;
 
  virtual void
  copyAndPermute
  (const SrcDistObject& sourceObj,
   const size_t numSameIDs,
   const Kokkos::DualView<const local_ordinal_type*,
     buffer_device_type>& permuteToLIDs,
   const Kokkos::DualView<const local_ordinal_type*,
     buffer_device_type>& permuteFromLIDs,
   const CombineMode CM) override;
 
  virtual void
  packAndPrepare
  (const SrcDistObject& sourceObj,
   const Kokkos::DualView<const local_ordinal_type*,
     buffer_device_type>& exportLIDs,
   Kokkos::DualView<packet_type*,
     buffer_device_type>& exports,
   Kokkos::DualView<size_t*,
     buffer_device_type> numPacketsPerLID,
   size_t& constantNumPackets) override;
 
  virtual void
  unpackAndCombine
  (const Kokkos::DualView<const local_ordinal_type*,
     buffer_device_type>& importLIDs,
   Kokkos::DualView<packet_type*,
     buffer_device_type> imports,
   Kokkos::DualView<size_t*,
     buffer_device_type> numPacketsPerLID,
   const size_t constantNumPackets,
   const CombineMode combineMode) override;
 
private:
  crs_graph_type graph_;
  Teuchos::RCP<crs_graph_type> graphRCP_;
  map_type rowMeshMap_;
  map_type domainPointMap_;
  map_type rangePointMap_;
  LO blockSize_;
 
  using graph_row_offset_host_type = typename crs_graph_type::local_graph_device_type::row_map_type::HostMirror;
  graph_row_offset_host_type ptrHost_;
 
  using graph_column_indices_host_type =   typename crs_graph_type::local_graph_device_type::entries_type::HostMirror;
  graph_column_indices_host_type indHost_;
 
  using impl_scalar_type_dualview = Kokkos::DualView<impl_scalar_type*, device_type>;
  using impl_scalar_type_wrapped_dualview = Details::WrappedDualView<impl_scalar_type_dualview>;
  mutable impl_scalar_type_wrapped_dualview val_;
 
  Teuchos::RCP<Teuchos::RCP<BMV> > X_colMap_;
  Teuchos::RCP<Teuchos::RCP<BMV> > Y_rowMap_;
 
  Teuchos::RCP<Teuchos::RCP<typename crs_graph_type::import_type> > pointImporter_;
 
  LO offsetPerBlock_;
 
  Teuchos::RCP<bool> localError_;
 
  Teuchos::RCP<Teuchos::RCP<std::ostringstream> > errs_;
 
  std::ostream& markLocalErrorAndGetStream ();
 
  // //! Clear the local error state and stream.
  // void clearLocalErrorStateAndStream ();
 
  template<class Device>
  struct is_cuda {
#if defined(KOKKOS_ENABLE_CUDA)
    // CudaHostPinnedSpace::execution_space ==
    // HostSpace::execution_space.  That's OK; it's host memory, that
    // just happens to be Cuda accessible.  But what if somebody gives
    // us Device<Cuda, CudaHostPinnedSpace>?  It looks like they mean
    // to run on device then, so we should sync to device.
    static constexpr bool value =
      std::is_same<typename Device::execution_space, Kokkos::Cuda>::value;
    // Gonna badly fake this here for other execspaces
#elif defined(KOKKOS_ENABLE_HIP)
    static constexpr bool value =
      std::is_same<typename Device::execution_space, Kokkos::Experimental::HIP>::value;
#elif defined(KOKKOS_ENABLE_SYCL)
    static constexpr bool value =
      std::is_same<typename Device::execution_space, Kokkos::Experimental::SYCL>::value;
#else
    static constexpr bool value = false;
#endif
  };
 
public:
#ifdef TPETRA_ENABLE_DEPRECATED_CODE
    // KK: sync modify syntax will not work
    //     the interface is deprecated bu the functionalities are removed
 
 
  inline void modify_host()
  {
    //throw std::logic_error("do not use");
  }
 
  inline void modify_device()
  {
    //throw std::logic_error("do not use");
  }
 
  template<class MemorySpace>
  void modify ()
  {
    //throw std::logic_error("do not use");
  }
 
  inline bool need_sync_host() const
  {
    //throw std::logic_error("do not use");
    return false; 
  }
 
  inline bool need_sync_device() const
  {
    //throw std::logic_error("do not use");
    return false; 
  }
 
  template<class MemorySpace>
  bool need_sync () const
  {
    //throw std::logic_error("do not use");
    return false;
  }
 
  inline void sync_host()
  {
    //throw std::logic_error("do not use");
  }
 
  inline void sync_device()
  {
    //throw std::logic_error("do not use");
  }
 
  template<class MemorySpace>
  void sync ()
  {
    //throw std::logic_error("do not use");
  }
#endif
 
 
    typename impl_scalar_type_dualview::t_host::const_type
    getValuesHost() const;
 
    typename impl_scalar_type_dualview::t_dev::const_type
    getValuesDevice() const;
 
#ifdef TPETRA_ENABLE_DEPRECATED_CODE
  template<class MemorySpace>
  typename std::conditional<is_cuda<MemorySpace>::value,
                            typename impl_scalar_type_dualview::t_dev,
                            typename impl_scalar_type_dualview::t_host>::type
  getValues () const
  {
    // Unlike std::conditional, if_c has a select method.
    return Kokkos::Impl::if_c<
      is_cuda<MemorySpace>::value,
      typename impl_scalar_type_dualview::t_dev,
      typename impl_scalar_type_dualview::t_host
      >::select (this->getValuesDeviceNonConst (), this->getValuesHostNonConst ());
  }
#endif
 
    typename impl_scalar_type_dualview::t_host
    getValuesHostNonConst() const;
 
    typename impl_scalar_type_dualview::t_dev
    getValuesDeviceNonConst() const;
 
    typename impl_scalar_type_dualview::t_host::const_type
    getValuesHost (const LO& lclRow) const;
 
    typename impl_scalar_type_dualview::t_dev::const_type
    getValuesDevice (const LO& lclRow) const;
 
    typename impl_scalar_type_dualview::t_host
    getValuesHostNonConst (const LO& lclRow);
 
    typename impl_scalar_type_dualview::t_dev
    getValuesDeviceNonConst (const LO& lclRow);
 
private:
 
  void
  applyBlockTrans (const BlockMultiVector<Scalar, LO, GO, Node>& X,
                   BlockMultiVector<Scalar, LO, GO, Node>& Y,
                   const Teuchos::ETransp mode,
                   const Scalar alpha,
                   const Scalar beta);
 
  void
  applyBlockNoTrans (const BlockMultiVector<Scalar, LO, GO, Node>& X,
                     BlockMultiVector<Scalar, LO, GO, Node>& Y,
                     const Scalar alpha,
                     const Scalar beta);
 
  void
  localApplyBlockNoTrans (const BlockMultiVector<Scalar, LO, GO, Node>& X,
                          BlockMultiVector<Scalar, LO, GO, Node>& Y,
                          const Scalar alpha,
                          const Scalar beta);
 
  LO
  findRelOffsetOfColumnIndex (const LO localRowIndex,
                              const LO colIndexToFind,
                              const LO hint = 0) const;
 
  LO offsetPerBlock () const;
 
  const_little_block_type
  getConstLocalBlockFromInput (const impl_scalar_type* val, const size_t pointOffset) const;
 
  little_block_type
  getNonConstLocalBlockFromInput (impl_scalar_type* val, const size_t pointOffset) const;
 
  little_block_host_type
  getNonConstLocalBlockFromInputHost (impl_scalar_type* val, const size_t pointOffset) const;
 
 
 
public:
  virtual Teuchos::RCP<const Teuchos::Comm<int> > getComm() const override;
 
 
  virtual global_size_t getGlobalNumCols() const override;
 
  virtual size_t getNodeNumCols() const override;
 
  virtual GO getIndexBase() const override;
 
  virtual global_size_t getGlobalNumEntries() const override;
 
  virtual size_t getNodeNumEntries() const override;
 
  virtual size_t getNumEntriesInGlobalRow (GO globalRow) const override;
 
  virtual size_t getGlobalMaxNumRowEntries () const override;
 
  virtual bool hasColMap () const override;
 
  virtual bool isLocallyIndexed () const override;
 
  virtual bool isGloballyIndexed () const override;
 
  virtual bool isFillComplete () const override;
 
  virtual bool supportsRowViews () const override;
 
 
 
 
  virtual void
  getGlobalRowCopy (GO GlobalRow,
                    nonconst_global_inds_host_view_type &Indices,
                    nonconst_values_host_view_type &Values,
                    size_t& NumEntries) const override;
#ifdef TPETRA_ENABLE_DEPRECATED_CODE
  virtual void
  getGlobalRowCopy (GO GlobalRow,
                    const Teuchos::ArrayView<GO> &Indices,
                    const Teuchos::ArrayView<Scalar> &Values,
                    size_t& NumEntries) const override;
#endif
 
#ifdef TPETRA_ENABLE_DEPRECATED_CODE
  virtual void
  getGlobalRowView (GO GlobalRow,
                    Teuchos::ArrayView<const GO>& indices,
                    Teuchos::ArrayView<const Scalar>& values) const override;
#endif // TPETRA_ENABLE_DEPRECATED_CODE
  virtual void
  getGlobalRowView (GO GlobalRow,
                    global_inds_host_view_type & indices,
                    values_host_view_type & values) const override;
 
  virtual void getLocalDiagCopy (::Tpetra::Vector<Scalar,LO,GO,Node>& diag) const override;
 
 
 
  virtual void leftScale (const ::Tpetra::Vector<Scalar, LO, GO, Node>& x) override;
 
  virtual void rightScale (const ::Tpetra::Vector<Scalar, LO, GO, Node>& x) override;
 
  virtual typename ::Tpetra::RowMatrix<Scalar, LO, GO, Node>::mag_type
  getFrobeniusNorm () const override;
};
 
} // namespace Tpetra
 
#endif // TPETRA_BLOCKCRSMATRIX_DECL_HPP