Tpetra_VbrMatrix_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_VBRMATRIX_DECL_HPP
#define TPETRA_VBRMATRIX_DECL_HPP


#include <Tpetra_ConfigDefs.hpp>
#include <Tpetra_DistObject_decl.hpp>
#include <Tpetra_Operator.hpp>
#include <Tpetra_VbrUtils.hpp>
#include <Kokkos_DefaultNode.hpp>
#include <Kokkos_DefaultKernels.hpp>

// Forward declarations, to avoid including files that we don't really
// need just for declaring VbrMatrix.  The #ifndef ... #endif prevents
// Doxygen from generating spurious documentation for the forward
// declarations.
#ifndef DOXYGEN_SHOULD_SKIP_THIS
namespace Tpetra {
  template<class LocalOrdinal, class GlobalOrdinal, class Node>
  class BlockMap;

  template<class LocalOrdinal, class GlobalOrdinal, class Node>
  class BlockCrsGraph;

  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  class Vector;
} // namespace Tpetra

namespace KokkosClassic {
  template<class Scalar, class LocalOrdinal, class Node>
  class VbrMatrix;
} // namespace KokkosClassic

namespace Teuchos {
  template<class OrdinalType, class ScalarType>
  class SerialDenseMatrix;
} // namespace Teuchos
#endif // DOXYGEN_SHOULD_SKIP_THIS

namespace Tpetra {


template<class Scalar = Operator<>::scalar_type,
         class LocalOrdinal = typename Operator<Scalar>::local_ordinal_type,
         class GlobalOrdinal = typename Operator<Scalar, LocalOrdinal>::global_ordinal_type,
         class Node = typename Operator<Scalar, LocalOrdinal, GlobalOrdinal>::node_type>
class VbrMatrix :
  public Tpetra::DistObject<char, LocalOrdinal, GlobalOrdinal, Node>,
  public Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
public:
  typedef Scalar        scalar_type;
  typedef LocalOrdinal  local_ordinal_type;
  typedef GlobalOrdinal global_ordinal_type;
  typedef Node          node_type;
  typedef typename KokkosClassic::DefaultKernels<Scalar, LocalOrdinal, Node>::BlockSparseOps sparse_ops_type;




  VbrMatrix (const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& blkRowMap,
             size_t maxNumEntriesPerRow,
             ProfileType pftype = DynamicProfile);


  VbrMatrix (const Teuchos::RCP<const BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& blkGraph);

  virtual ~VbrMatrix();





  template <class DomainScalar, class RangeScalar>
  void
  multiply (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X,
            MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
            Teuchos::ETransp trans,
            RangeScalar alpha,
            RangeScalar beta) const;


  template <class DomainScalar, class RangeScalar>
  void
  solve (const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
         MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X,
         Teuchos::ETransp trans) const;





  Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getDomainMap() const;


  Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getRangeMap() const;


  void
  apply (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
         MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
         Teuchos::ETransp trans = Teuchos::NO_TRANS,
         Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
         Scalar beta = Teuchos::ScalarTraits<Scalar>::zero()) const;


  void
  applyInverse (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & Y,
                MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
                Teuchos::ETransp trans) const;


  bool hasTransposeApply() const;




  Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > getBlockRowMap() const;

  Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > getBlockColMap() const;

  Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > getBlockDomainMap() const;

  Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > getBlockRangeMap() const;

  Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getPointRowMap() const;

  Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getPointColMap() const;

  bool isFillComplete() const;





  void putScalar(Scalar s);


  void setGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, const Teuchos::SerialDenseMatrix<int,Scalar>& blockEntry);


  void setLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, const Teuchos::SerialDenseMatrix<int,Scalar>& blockEntry);


  void sumIntoGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, const Teuchos::SerialDenseMatrix<int,Scalar>& blockEntry);


  void sumIntoLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, const Teuchos::SerialDenseMatrix<int,Scalar>& blockEntry);


  void setGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);


  void setLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);


  void sumIntoGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);


  void sumIntoLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);





  void fillComplete(const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& blockDomainMap, const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& blockRangeMap);


  void fillComplete();

  void globalAssemble();





  void getGlobalBlockRowView(GlobalOrdinal globalBlockRow,
                             LocalOrdinal& numPtRows,
                             Teuchos::ArrayView<const GlobalOrdinal>& blockCols,
                             Teuchos::Array<LocalOrdinal>& ptColsPerBlockCol,
                             Teuchos::Array<Teuchos::ArrayRCP<const Scalar> >& blockEntries) const;


  void getLocalBlockRowView(LocalOrdinal localBlockRow,
                            LocalOrdinal& numPtRows,
                            Teuchos::ArrayView<const LocalOrdinal>& blockCols,
                            Teuchos::Array<LocalOrdinal>& ptColsPerBlockCol,
                            Teuchos::ArrayRCP<const Scalar>& blockEntries) const;


  void getGlobalBlockEntryView(GlobalOrdinal globalBlockRow,
                               GlobalOrdinal globalBlockCol,
                               LocalOrdinal& numPtRows,
                               LocalOrdinal& numPtCols,
                               Teuchos::ArrayRCP<const Scalar>& blockEntry) const;


  void getGlobalBlockEntryViewNonConst(GlobalOrdinal globalBlockRow,
                                       GlobalOrdinal globalBlockCol,
                                       LocalOrdinal& numPtRows,
                                       LocalOrdinal& numPtCols,
                                       Teuchos::ArrayRCP<Scalar>& blockEntry);


  void getLocalBlockEntryView(LocalOrdinal localBlockRow,
                              LocalOrdinal localBlockCol,
                              LocalOrdinal& numPtRows,
                              LocalOrdinal& numPtCols,
                              Teuchos::ArrayRCP<const Scalar>& blockEntry) const;


  void getLocalBlockEntryViewNonConst(LocalOrdinal localBlockRow,
                                      LocalOrdinal localBlockCol,
                                      LocalOrdinal& numPtRows,
                                      LocalOrdinal& numPtCols,
                                      Teuchos::ArrayRCP<Scalar>& blockEntry);


  void getLocalDiagCopy (Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& diag) const;

  Teuchos::RCP<const BlockCrsGraph<LocalOrdinal, GlobalOrdinal, Node> >
  getBlockCrsGraph () {
    return constBlkGraph_;
  }




  virtual bool checkSizes (const SrcDistObject& source);

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

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

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




  std::string description() const;

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

 private:
  Teuchos::RCP<Node> getNode() const;

  void updateImport(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X) const;
  void updateExport(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const;

  void createImporterExporter();
  void optimizeStorage();
  void fillLocalMatrix();
  void fillLocalMatVec();

  //private data members:

  //We hold two graph pointers, one const and the other non-const.
  //If a BlockCrsGraph is provided at construction, it is const and VbrMatrix
  //never changes it.
  //If a BlockCrsGraph is not provided at construction, VbrMatrix creates one
  //internally and fills it as the matrix is filled, up until fillComplete()
  //is called.
  //
  //blkGraph_ is either the internally created graph, or is null.
  //constBlkGraph_ is either a pointer to blkGraph_, or a pointer to the
  //graph provided at construction time.
  //
  Teuchos::RCP<BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node> > blkGraph_;
  Teuchos::RCP<const BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node> > constBlkGraph_;

  typedef KokkosClassic::VbrMatrix<Scalar, LocalOrdinal, Node> local_matrix_type;

  // We use a pointer so that the _decl header file only needs a
  // forward declaration, not an include.  We don't really need an RCP
  // here, since the local matrix object never gets shared outside the
  // class.  std::unique_ptr (C++11) would be more appropriate.
  Teuchos::RCP<local_matrix_type> lclMatrix_;

  //A variable-block-row matrix is represented by 6 arrays
  //in packed (contiguous storage) form. For a description of these
  //arrays, see the text at the bottom of this file.
  //(2 of those arrays, rptr and cptr, are represented by arrays in the
  //getBlockRowMap() and getBlockColMap() objects, and
  //another two of those arrays, bptr and bindx, are represented by arrays in
  //the BlockCrsGraph object.)
  //This is noted in the comments for rptr,cptr,bptr,bindx below.
  //
  //These arrays are handled as if they may point to memory that resides on
  //a separate device (e.g., a GPU). In other words, when the contents of these
  //arrays are manipulated, we use views or buffers obtained from the Node
  //object.
  Teuchos::ArrayRCP<Scalar> pbuf_values1D_;
  Teuchos::ArrayRCP<LocalOrdinal> pbuf_indx_;

  sparse_ops_type lclMatOps_;
  Teuchos::RCP<Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer_;
  Teuchos::RCP<Tpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter_;
  mutable Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > importedVec_;
  mutable Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > exportedVec_;

  typedef typename std::map<GlobalOrdinal,Teuchos::ArrayRCP<Scalar> > RowGlobalCols;

  //We use an array-of-maps to represent the variable-block-row matrix in
  //un-packed '2D' form.
  //
  //This unpacked data is assumed to be resident in CPU (host) memory.
  //It doesn't make sense to copy this data back and forth to a separate
  //compute device (e.g., a GPU), since we don't support doing matrix-vector
  //products until after fillComplete is called, at which time contiguous
  //arrays are allocated on the device and matrix data is copied into them.
  Teuchos::RCP<Teuchos::Array<RowGlobalCols> > data_2D_;

  VbrUtils::VbrData<LocalOrdinal,GlobalOrdinal,Scalar> nonlocal_data_;

  bool is_fill_completed_;
  bool is_storage_optimized_;
};//class VbrMatrix

}//namespace Tpetra

//----------------------------------------------------------------------------
// Description of arrays representing the VBR format:
//
// (For more description, see this URL (valid as of 5/26/2010):
// http://docs.sun.com/source/817-0086-10/prog-sparse-support.html)
// ...and of course more can be found using google...
// The old Aztec manual was a great resource for this but I can't
// find a copy of that these days...
//
//
// Here is a brief description of the 6 arrays that are required to
// represent a VBR matrix in packed (contiguous-memory-storage) format:
//
// rptr: length num_block_rows + 1
//       rptr[i]: the pt-row corresponding to the i-th block-row
//       Note: rptr is getBlockRowMap()->getNodeFirstPointInBlocks().
//
// cptr: length num_distinct_block_cols + 1
//       cptr[j]: the pt-col corresponding to the j-th block-col
//       Note: cptr is getBlockColMap()->getNodeFirstPointInBlocks().
//
// bptr: length num_block_rows + 1
//       bptr[i]: location in bindx of the first nonzero block-entry
//                of the i-th block-row
//       Note: bptr is blkGraph_->getNodeRowOffsets();
//
// bindx: length num-nonzero-block-entries
//        bindx[j]: block-col-index of j-th block-entry
//        Note: bindx is blkGraph_->getNodePackedIndices();
//
// indx: length num-nonzero-block-entries + 1
//       indx[j]: location in vals of the beginning of the j-th
//       block-entry
//
// vals: length num-nonzero-scalar-entries
//
//
// Some example look-ups:
//
// int nbr = num_block_rows;
// int total_num_block_nonzeros = bptr[nbr];
// int total_num_scalar_nonzeros = indx[num_block_nonzeros];
//
// //get arrays for i-th block-row:
// int* bindx_i = &bindx[bptr[i]];
// double* vals_i = &val[indx[bptr[i]]];
// int num_block_nonzeros_in_row_i = bptr[i+1]-bptr[i];
//
//----------------------------------------------------------------------------

#endif //TPETRA_VBRMATRIX_DECL_HPP