Tpetra_CrsGraph_decl.hpp

// @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_CRSGRAPH_DECL_HPP
#define TPETRA_CRSGRAPH_DECL_HPP

#include <Teuchos_Describable.hpp>
#include <Teuchos_CompileTimeAssert.hpp>

#include <Kokkos_DefaultNode.hpp>
#include <Kokkos_DefaultKernels.hpp>
#include <Kokkos_CrsGraph.hpp>
#include <Kokkos_NodeHelpers.hpp>

#include "Tpetra_ConfigDefs.hpp"
#include "Tpetra_RowGraph.hpp"
#include "Tpetra_DistObject.hpp"
#include "Tpetra_Util.hpp"

namespace Tpetra {

#ifndef DOXYGEN_SHOULD_SKIP_THIS
  // forward declaration
  template <class S, class LO, class GO, class N, class SpMatOps>
  class CrsMatrix;
#endif

  struct RowInfo {
    size_t localRow;
    size_t allocSize;
    size_t numEntries;
    size_t offset1D;
  };

  enum ELocalGlobal {
    LocalIndices,
    GlobalIndices
  };


  template <class LocalOrdinal, 
            class GlobalOrdinal = LocalOrdinal, 
            class Node = Kokkos::DefaultNode::DefaultNodeType,
            class LocalMatOps = typename Kokkos::DefaultKernels<void,LocalOrdinal,Node>::SparseOps >
  class CrsGraph : 
    public RowGraph<LocalOrdinal,GlobalOrdinal,Node>,
    public DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>,
    public Teuchos::ParameterListAcceptorDefaultBase
  {
    template <class S, class LO, class GO, class N, class SpMatOps>
    friend class CrsMatrix;

  public: 
    typedef LocalOrdinal                         local_ordinal_type;
    typedef GlobalOrdinal                        global_ordinal_type;
    typedef Node                                 node_type;
    typedef Map<LocalOrdinal,GlobalOrdinal,Node> map_type;



    CrsGraph (const Teuchos::RCP<const map_type>& rowMap, 
              size_t maxNumEntriesPerRow, 
              ProfileType pftype = DynamicProfile,
              const Teuchos::RCP<Teuchos::ParameterList>& plist = Teuchos::null);

    CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap, 
              const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc, 
              ProfileType pftype = DynamicProfile,
              const Teuchos::RCP<Teuchos::ParameterList>& plist = Teuchos::null);

    CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap, 
              const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& colMap, 
              size_t maxNumEntriesPerRow, 
              ProfileType pftype = DynamicProfile,
              const Teuchos::RCP<Teuchos::ParameterList>& plist = Teuchos::null);

    CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap, 
              const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& colMap, 
              const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, 
              ProfileType pftype = DynamicProfile,
              const Teuchos::RCP<Teuchos::ParameterList>& plist = Teuchos::null);

    virtual ~CrsGraph();




    void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist);

    Teuchos::RCP<const Teuchos::ParameterList>
    getValidParameters () const;





      void insertGlobalIndices(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &indices);


      void insertLocalIndices(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &indices);


      void removeLocalIndices(LocalOrdinal localRow);




      void globalAssemble();

      void resumeFill();

      void fillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap, OptimizeOption os = DoOptimizeStorage);

      void fillComplete(OptimizeOption os = DoOptimizeStorage);




      const RCP<const Comm<int> > & getComm() const;

      RCP<Node> getNode() const;

      const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getRowMap() const;

      const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getColMap() const;

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

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

      RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > getImporter() const;

      RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > getExporter() const;


      global_size_t getGlobalNumRows() const;


      global_size_t getGlobalNumCols() const;

      size_t getNodeNumRows() const;


      size_t getNodeNumCols() const;

      GlobalOrdinal getIndexBase() const;


      global_size_t getGlobalNumEntries() const;

      size_t getNodeNumEntries() const;


      size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const;


      size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const;


      size_t getNodeAllocationSize() const;


      size_t getNumAllocatedEntriesInGlobalRow(GlobalOrdinal globalRow) const;


      size_t getNumAllocatedEntriesInLocalRow(LocalOrdinal localRow) const;


      global_size_t getGlobalNumDiags() const;


      size_t getNodeNumDiags() const;


      size_t getGlobalMaxNumRowEntries() const;


      size_t getNodeMaxNumRowEntries() const;

      bool hasColMap() const; 


      bool isLowerTriangular() const;


      bool isUpperTriangular() const;

      bool isLocallyIndexed() const;

      bool isGloballyIndexed() const;

      bool isFillComplete() const;

      bool isFillActive() const;


      bool isSorted() const;


      bool isStorageOptimized() const;

      ProfileType getProfileType() const;


      void getGlobalRowCopy(GlobalOrdinal GlobalRow, 
                            const ArrayView<GlobalOrdinal> &Indices, 
                            size_t &NumIndices
                            ) const;


      void getLocalRowCopy(LocalOrdinal LocalRow, 
                           const ArrayView<LocalOrdinal> &indices, 
                           size_t &NumIndices
                           ) const;


      void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal> &Indices) const;


      void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal> &indices) const;




      std::string description() const;

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




      bool checkSizes(const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>& source);

      void copyAndPermute(const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node> & source,
                          size_t numSameIDs,
                          const ArrayView<const LocalOrdinal> &permuteToLIDs,
                          const ArrayView<const LocalOrdinal> &permuteFromLIDs);

      void packAndPrepare(const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node> & source,
                          const ArrayView<const LocalOrdinal> &exportLIDs,
                          Array<GlobalOrdinal> &exports,
                          const ArrayView<size_t> & numPacketsPerLID,
                          size_t& constantNumPackets,
                          Distributor &distor);

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




      ArrayRCP<const size_t>       getNodeRowBegs() const;


      ArrayRCP<const LocalOrdinal> getNodePackedIndices() const;




      TPETRA_DEPRECATED void optimizeStorage();

      TPETRA_DEPRECATED ArrayRCP<const GlobalOrdinal> getGlobalRowView(GlobalOrdinal GlobalRow) const;

      TPETRA_DEPRECATED ArrayRCP<const LocalOrdinal> getLocalRowView(LocalOrdinal LocalRow) const;


    private:
      // We forbid copy construction by declaring this method private
      // and not implementing it.
      CrsGraph (const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> &Source);

      // We forbid assignment (operator=) by declaring this method
      // private and not implementing it.
      CrsGraph<LocalOrdinal,GlobalOrdinal,Node>& 
      operator= (const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> &rhs);

    protected:
      // these structs are conveniences, to cut down on the number of
      // arguments to some of the methods below.
      struct SLocalGlobalViews {
        ArrayView<const GlobalOrdinal> ginds;
        ArrayView<const LocalOrdinal>  linds;
      };
      struct SLocalGlobalNCViews {
        ArrayView<GlobalOrdinal>       ginds;
        ArrayView<LocalOrdinal>        linds;
      };
      //
      // Allocation
      //
      bool indicesAreAllocated() const;

      void allocateIndices (ELocalGlobal lg);

      template<class T>
      void 
      allocateValues (ArrayRCP<T> &values1D, ArrayRCP<ArrayRCP<T> > &values2D) const;

      template<ELocalGlobal lg>
      RowInfo 
      updateAlloc (RowInfo rowinfo, size_t allocSize);
      template<ELocalGlobal lg, class T> 
      RowInfo 
      updateAllocAndValues (RowInfo rowinfo, size_t allocSize, ArrayRCP<T> &rowVals);
      //
      // Local versus global indices
      //
      void computeIndexState();
      void makeColMap();
      void makeIndicesLocal();
      void makeImportExport();
      //
      // insert/suminto/replace
      //
      template<ELocalGlobal lg>
      size_t 
      filterIndices (const SLocalGlobalNCViews &inds) const;

      template<ELocalGlobal lg, class T>
      size_t 
      filterIndicesAndValues (const SLocalGlobalNCViews &inds, const ArrayView<T> &vals) const;

      template<ELocalGlobal lg, ELocalGlobal I> 
      size_t
      insertIndices (RowInfo rowInfo, const SLocalGlobalViews &newInds);

      template<ELocalGlobal lg, ELocalGlobal I, class IterO, class IterN> 
      void
      insertIndicesAndValues (RowInfo rowInfo, const SLocalGlobalViews &newInds, IterO rowVals, IterN newVals);

      template<ELocalGlobal lg, class IterO, class IterN, class BinaryFunction> 
      void 
      transformValues (RowInfo rowInfo, const SLocalGlobalViews &inds, IterO rowVals, IterN newVals, BinaryFunction f) const;
      //
      // Sorting and merging
      //
      bool                       isMerged() const;
      void                       setSorted(bool sorted);
      void                       setMerged(bool merged);
      void                       sortAllIndices();
      void                       sortRowIndices(RowInfo rowinfo);
      template <class Scalar> void sortRowIndicesAndValues(RowInfo rowinfo, ArrayView<Scalar> values);
      void                                             mergeAllIndices();
      void                                             mergeRowIndices(RowInfo rowinfo);
      template <class Iter, class BinaryFunction> void mergeRowIndicesAndValues(RowInfo rowinfo, Iter rowValueIter, BinaryFunction f);
      // 
      void setDomainRangeMaps(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap);
      void staticAssertions() const;
      // global consts
      void clearGlobalConstants();
      void computeGlobalConstants();
      // graph data accessors
      RowInfo                         getRowInfo(size_t myRow) const;
      ArrayView<const LocalOrdinal>   getLocalView(RowInfo rowinfo) const;
      ArrayView<LocalOrdinal>         getLocalViewNonConst(RowInfo rowinfo);
      ArrayView<const GlobalOrdinal>  getGlobalView(RowInfo rowinfo) const;
      ArrayView<GlobalOrdinal>        getGlobalViewNonConst(RowInfo rowinfo);
      size_t                          findLocalIndex(RowInfo rowinfo, LocalOrdinal ind) const;
      size_t                          findGlobalIndex(RowInfo rowinfo, GlobalOrdinal ind) const;
      // local Kokkos objects
      void pushToLocalGraph();
      void pullFromLocalGraph();
      void fillLocalGraph(OptimizeOption os);
      const Kokkos::CrsGraph<LocalOrdinal,Node,LocalMatOps> & getLocalGraph() const;
      Kokkos::CrsGraph<LocalOrdinal,Node,LocalMatOps>       & getLocalGraphNonConst();
      // debugging
      void checkInternalState() const;

      // Tpetra support objects
      RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap_, colMap_, rangeMap_, domainMap_;
      RCP<Import<LocalOrdinal,GlobalOrdinal,Node> > importer_;
      RCP<Export<LocalOrdinal,GlobalOrdinal,Node> > exporter_;

      // local data, stored in a Kokkos::CrsGraph. only initialized after fillComplete()
      Kokkos::CrsGraph<LocalOrdinal,Node,LocalMatOps> lclGraph_;

      // Local and Global Counts
      // nodeNumEntries_ and nodeNumAllocated_ are required to be always consistent
      // nodeMaxNumEntries_, nodeNumDiags_ and the global quantities are computed during fillComplete() and only valid when isFillComplete()
      global_size_t globalNumEntries_, globalNumDiags_, globalMaxNumRowEntries_;
      size_t          nodeNumEntries_,   nodeNumDiags_,   nodeMaxNumRowEntries_, nodeNumAllocated_;

      // allocate static or dynamic?
      ProfileType pftype_;
      // requested allocation sizes; we have to preserve these, because we perform late-allocation
      // number of non-zeros to allocate per row; set to null after they are allocated.
      ArrayRCP<const size_t> numAllocPerRow_;
      // number of non-zeros to allocate for all row; either this or numAllocPerRow_ is used, but not both.
      size_t numAllocForAllRows_;

      // graph indices. before allocation, all are null. 
      // after allocation, except during makeIndicesLocal(), one of local or global is null.
      // we will never have 1D and 2D structures being non-null
      // this is host memory
      // 1D == StaticAllocation, 2D == DynamicAllocation
      //
      // 1D/Static structures
      //
      ArrayRCP< LocalOrdinal>                     lclInds1D_;
      ArrayRCP<GlobalOrdinal>                     gblInds1D_;
      // offset to the beg and end entries of each row. only used for 1D (Static) allocation.
      // i.e., indices for row R are lclInds1D_[i] for i in [b,e) where b = rowBegs_[R] and e = rowEnds_[R]
      // for a packed (optimized) allocation, we will typically have &rowBegs_[R]+1 == &rowEnds_[R]
      // these are null for 2D (Dynamic) allocations
      // the allocation size is computed by looking at the difference between rowBegs_[r+1] and rowBegs_[r]
      // rowBegs_ therefore has length N+1, while rowEnds_ has length N
      ArrayRCP<size_t> rowBegs_, rowEnds_;
      //
      // 2D/Dynamic structures. 
      //
      ArrayRCP<ArrayRCP< LocalOrdinal> > lclInds2D_;
      ArrayRCP<ArrayRCP<GlobalOrdinal> > gblInds2D_;
      ArrayRCP<size_t>       numEntriesPerRow_;

      // TODO: these might be useful in the future
      // ArrayRCP< typedef ArrayRCP<const GlobalOrdinal>::iterator > gRowPtrs_;
      // ArrayRCP< typedef ArrayRCP<GlobalOrdinal>::iterator > gRowPtrsNC_;
      // ArrayRCP< typedef ArrayRCP<const LocalOrdinal>::iterator > lRowPtrs_;
      // ArrayRCP< typedef ArrayRCP<LocalOrdinal>::iterator > lRowPtrsNC_;

      bool indicesAreAllocated_,
           indicesAreLocal_,
           indicesAreGlobal_,
           fillComplete_, 
           lowerTriangular_,
           upperTriangular_,
           indicesAreSorted_,
           noRedundancies_,
           haveGlobalConstants_;

      // non-local data
      std::map<GlobalOrdinal, std::deque<GlobalOrdinal> > nonlocals_;

  }; // class CrsGraph

} // namespace Tpetra

#endif