Tpetra_Experimental_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.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
// @HEADER

#ifndef TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_DECL_HPP
#define TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_DECL_HPP


#include "Tpetra_CrsGraph.hpp"
#include "Tpetra_RowMatrix.hpp"
#include "Tpetra_Experimental_BlockMultiVector.hpp"
#include "Tpetra_CrsMatrix_decl.hpp"

namespace Tpetra {
namespace Experimental {
namespace Classes {

template<class Scalar = ::Tpetra::Details::DefaultTypes::scalar_type,
         class LO = ::Tpetra::Details::DefaultTypes::local_ordinal_type,
         class GO = ::Tpetra::Details::DefaultTypes::global_ordinal_type,
         class Node = ::Tpetra::Details::DefaultTypes::node_type>
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 TPETRA_DEPRECATED
  getLocalDiagOffsets (Teuchos::ArrayRCP<size_t>& 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;

  void TPETRA_DEPRECATED
  getLocalDiagCopy (BlockCrsMatrix<Scalar,LO,GO,Node>& 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;



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

  virtual void
  copyAndPermute (const ::Tpetra::SrcDistObject& source,
                  size_t numSameIDs,
                  const Teuchos::ArrayView<const LO>& permuteToLIDs,
                  const Teuchos::ArrayView<const LO>& permuteFromLIDs);

  virtual void
  packAndPrepare (const ::Tpetra::SrcDistObject& source,
                  const Teuchos::ArrayView<const LO>& exportLIDs,
                  Teuchos::Array<packet_type>& exports,
                  const Teuchos::ArrayView<size_t>& numPacketsPerLID,
                  size_t& constantNumPackets,
                  ::Tpetra::Distributor& distor);

  virtual void
  unpackAndCombine (const Teuchos::ArrayView<const LO> &importLIDs,
                    const Teuchos::ArrayView<const packet_type> &imports,
                    const Teuchos::ArrayView<size_t> &numPacketsPerLID,
                    size_t constantNumPackets,
                    ::Tpetra::Distributor& distor,
                    ::Tpetra::CombineMode CM);

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

public:


  template<class MemorySpace>
  void modify ()
  {
    // It's legit to use a memory space, execution space, or
    // Kokkos::Device specialization as the template parameter of
    // Kokkos::DualView::modify.  That method just extracts the
    // memory_space typedef of its template parameter anyway.
    // However, insisting on a memory space avoids unnecessary
    // instantiations.
    val_.template modify<typename MemorySpace::memory_space> ();
  }

  template<class MemorySpace>
  bool need_sync () const
  {
    // It's legit to use a memory space, execution space, or
    // Kokkos::Device specialization as the template parameter of
    // Kokkos::DualView::need_sync.  That method just extracts the
    // memory_space typedef of its template parameter anyway.
    // However, insisting on a memory space avoids unnecessary
    // instantiations.
    return val_.template need_sync<typename MemorySpace::memory_space> ();
  }

  template<class MemorySpace>
  void sync ()
  {
    // It's legit to use a memory space, execution space, or
    // Kokkos::Device specialization as the template parameter of
    // Kokkos::DualView::sync.  That method just extracts the
    // memory_space typedef of its template parameter anyway.
    // However, insisting on a memory space avoids unnecessary
    // instantiations.
    val_.template sync<typename MemorySpace::memory_space> ();
  }

  template<class MemorySpace>
  auto getValues () -> decltype (val_.template view<typename MemorySpace::memory_space> ())
  {
    return val_.template view<typename MemorySpace::memory_space> ();
  }


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 Teuchos::RCP<Node> getNode() 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 global_size_t TPETRA_DEPRECATED getGlobalNumDiags() const;

  virtual size_t TPETRA_DEPRECATED getNodeNumDiags() const;

  virtual size_t getGlobalMaxNumRowEntries() const;

  virtual bool hasColMap() const;

  virtual bool TPETRA_DEPRECATED isLowerTriangular () const;

  virtual bool TPETRA_DEPRECATED isUpperTriangular () 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 Classes
} // namespace Experimental
} // namespace Tpetra

#endif // TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_DECL_HPP