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 Kokkos::View<const impl_scalar_type**,
                       Kokkos::LayoutRight,
                       device_type,
                       Kokkos::MemoryTraits<Kokkos::Unmanaged> >
          const_little_block_type;
  typedef Kokkos::View<impl_scalar_type*,
                       Kokkos::LayoutRight,
                       device_type,
                       Kokkos::MemoryTraits<Kokkos::Unmanaged> >
          little_vec_type;
  typedef Kokkos::View<const impl_scalar_type*,
                       Kokkos::LayoutRight,
                       device_type,
                       Kokkos::MemoryTraits<Kokkos::Unmanaged> >
          const_little_vec_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;

  Teuchos::RCP<const map_type> getRangeMap () const;

  Teuchos::RCP<const map_type> getRowMap () const;

  Teuchos::RCP<const map_type> getColMap () const;

  global_size_t getGlobalNumRows() const;

  size_t getNodeNumRows() const;

  size_t getNodeMaxNumRowEntries() const;

  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;

  bool hasTransposeApply () const {
    // 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;

  void
  describe (Teuchos::FancyOStream& out,
            const Teuchos::EVerbosityLevel verbLevel) const;



  LO getBlockSize () const { return blockSize_; }

  virtual Teuchos::RCP<const ::Tpetra::RowGraph<LO,GO,Node> > getGraph () const;

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

  void
  gaussSeidelCopy (MultiVector<Scalar,LO,GO,Node> &X,
                   const ::Tpetra::MultiVector<Scalar,LO,GO,Node> &B,
                   const ::Tpetra::MultiVector<Scalar,LO,GO,Node> &D,
                   const Scalar& dampingFactor,
                   const ESweepDirection direction,
                   const int numSweeps,
                   const bool zeroInitialGuess) const;

  void
  reorderedGaussSeidelCopy (::Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
                            const ::Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
                            const ::Tpetra::MultiVector<Scalar,LO,GO,Node>& D,
                            const Teuchos::ArrayView<LO>& rowIndices,
                            const Scalar& dampingFactor,
                            const ESweepDirection direction,
                            const int numSweeps,
                            const bool zeroInitialGuess) const;

  void
  localGaussSeidel (const BlockMultiVector<Scalar, LO, GO, Node>& Residual,
                          BlockMultiVector<Scalar, LO, GO, Node>& Solution,
                    const Kokkos::View<impl_scalar_type***, device_type,
                          Kokkos::MemoryUnmanaged>& D_inv,
                    const Scalar& omega,
                    const ESweepDirection direction) const;

  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;

  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;

  void
  getLocalRowCopy (LO LocalRow,
                   const Teuchos::ArrayView<LO> &Indices,
                   const Teuchos::ArrayView<Scalar> &Values,
                   size_t &NumEntries) const;

  little_block_type
  getLocalBlock (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
  sumIntoLocalValuesByOffsets (const LO localRowInd,
                               const ptrdiff_t offsets[],
                               const Scalar vals[],
                               const LO numOffsets) const;

  size_t getNumEntriesInLocalRow (const LO localRowInd) const;

  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;

  void
  getLocalDiagCopy (const Kokkos::View<impl_scalar_type***, device_type,
                                       Kokkos::MemoryUnmanaged>& diag,
                    const Teuchos::ArrayView<const size_t>& offsets) const;


protected:
  LO
  absMaxLocalValues (const LO localRowInd,
                     const LO colInds[],
                     const Scalar vals[],
                     const LO numColInds) const;

  LO
  absMaxLocalValuesByOffsets (const LO localRowInd,
                              const ptrdiff_t offsets[],
                              const Scalar vals[],
                              const LO numOffsets) const;



  using buffer_device_type = typename DistObject<Scalar, LO, GO,
                                                 Node>::buffer_device_type;

  virtual bool checkSizes (const ::Tpetra::SrcDistObject& source);

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

  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,
   Distributor& /* distor */);

  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,
   Distributor& /* distor */,
   const CombineMode combineMode);

private:
  crs_graph_type graph_;
  Teuchos::RCP<crs_graph_type> graphRCP_;
  map_type rowMeshMap_;
  map_type domainPointMap_;
  map_type rangePointMap_;
  LO blockSize_;

  typename crs_graph_type::local_graph_type::row_map_type::HostMirror ptrHost_;

  typename crs_graph_type::local_graph_type::entries_type::HostMirror indHost_;

  typename Kokkos::DualView<impl_scalar_type*, device_type> 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;
#else
    static constexpr bool value = false;
#endif // defined(KOKKOS_ENABLE_CUDA)
  };

public:


  inline void modify_host()
  {
    val_.modify_host();
  }

  inline void modify_device()
  {
    val_.modify_device();
  }

  template<class MemorySpace>
  void modify ()
  {
    if (is_cuda<MemorySpace>::value) {
      this->modify_device ();
    }
    else {
      this->modify_host ();
    }
  }

  inline bool need_sync_host() const
  {
    return val_.need_sync_host();
  }

  inline bool need_sync_device() const
  {
    return val_.need_sync_device();
  }

  template<class MemorySpace>
  bool need_sync () const
  {
    if (is_cuda<MemorySpace>::value) {
      return this->need_sync_device ();
    }
    else {
      return this->need_sync_host ();
    }
  }

  inline void sync_host()
  {
    val_.sync_host();
  }

  inline void sync_device()
  {
    val_.sync_device();
  }

  template<class MemorySpace>
  void sync ()
  {
    if (is_cuda<MemorySpace>::value) {
      this->sync_device ();
    }
    else {
      this->sync_host ();
    }
  }

  // \brief Get the host view of the matrix's values
  typename Kokkos::DualView<impl_scalar_type*, device_type>::t_host getValuesHost () const {
    return val_.view_host();
  }

  // \brief Get the device view of the matrix's values
  typename Kokkos::DualView<impl_scalar_type*, device_type>::t_dev getValuesDevice () const {
    return val_.view_device();
  }

  template<class MemorySpace>
  typename std::conditional<is_cuda<MemorySpace>::value,
                            typename Kokkos::DualView<impl_scalar_type*, device_type>::t_dev,
                            typename Kokkos::DualView<impl_scalar_type*, device_type>::t_host>::type
  getValues ()
  {
    // Unlike std::conditional, if_c has a select method.
    return Kokkos::Impl::if_c<
        is_cuda<MemorySpace>::value,
        typename Kokkos::DualView<impl_scalar_type*, device_type>::t_dev,
        typename Kokkos::DualView<impl_scalar_type*, device_type>::t_host
      >::select (this->getValuesDevice (), this->getValuesHost ());
  }


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;

  const_little_block_type
  getConstLocalBlockFromAbsOffset (const size_t absBlockOffset) const;

  little_block_type
  getNonConstLocalBlockFromAbsOffset (const size_t absBlockOffset) const;

  const_little_block_type
  getConstLocalBlockFromRelOffset (const LO lclMeshRow,
                                   const size_t relMeshOffset) const;

public:
  virtual Teuchos::RCP<const Teuchos::Comm<int> > getComm() const;


  virtual global_size_t getGlobalNumCols() const;

  virtual size_t getNodeNumCols() const;

  virtual GO getIndexBase() const;

  virtual global_size_t getGlobalNumEntries() const;

  virtual size_t getNodeNumEntries() const;

  virtual size_t getNumEntriesInGlobalRow (GO globalRow) const;

  virtual size_t getGlobalMaxNumRowEntries () const;

  virtual bool hasColMap () const;

  virtual bool isLocallyIndexed () const;

  virtual bool isGloballyIndexed () const;

  virtual bool isFillComplete () const;

  virtual bool supportsRowViews () const;




  virtual void
  getGlobalRowCopy (GO GlobalRow,
                    const Teuchos::ArrayView<GO> &Indices,
                    const Teuchos::ArrayView<Scalar> &Values,
                    size_t& NumEntries) const;

  virtual void
  getGlobalRowView (GO GlobalRow,
                    Teuchos::ArrayView<const GO>& indices,
                    Teuchos::ArrayView<const Scalar>& values) const;

  virtual void getLocalDiagCopy (::Tpetra::Vector<Scalar,LO,GO,Node>& diag) const;



  virtual void leftScale (const ::Tpetra::Vector<Scalar, LO, GO, Node>& x);

  virtual void rightScale (const ::Tpetra::Vector<Scalar, LO, GO, Node>& x);

  virtual typename ::Tpetra::RowMatrix<Scalar, LO, GO, Node>::mag_type
  getFrobeniusNorm () const;
};

} // namespace Tpetra

#endif // TPETRA_BLOCKCRSMATRIX_DECL_HPP