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 "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 RCP<const map_type>& rowMap,
              size_t maxNumEntriesPerRow,
              ProfileType pftype = DynamicProfile,
              const RCP<ParameterList>& params = null);

    CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap,
              const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
              ProfileType pftype = DynamicProfile,
              const RCP<ParameterList>& params = 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 RCP<ParameterList>& params = 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 RCP<ParameterList>& params = null);

    virtual ~CrsGraph();




    void setParameterList (const RCP<ParameterList>& params);

    RCP<const 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(const RCP<ParameterList> &params = null);

      void fillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap, const RCP<ParameterList> &params = null);

      void fillComplete(const RCP<ParameterList> &params = null);




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


      ArrayRCP<const LocalOrdinal> getNodePackedIndices() 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:
      typedef typename LocalMatOps::template graph<LocalOrdinal,Node>::graph_type local_graph_type;

      // 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>
      ArrayRCP<T> allocateValues1D () const;
      template <class T>
      ArrayRCP<ArrayRCP<T> > allocateValues2D () 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 fillLocalGraph(const RCP<ParameterList> &params);
      const RCP<const local_graph_type> getLocalGraph() const;
      const RCP<local_graph_type> 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()
      RCP<local_graph_type> 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 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 = rowPtrs_[R] and e = rowPtrs_[R+1]
      // only the first numRowEntries_[R] of these are valid
      // both of these are null for 2D (Dynamic) allocations
      // rowPtrs_ has length N+1, while numRowEntries_ has length N
      // we may delete this to save memory on fillComplete, if "Delete Row Pointers" is speified
      ArrayRCP<size_t> rowPtrs_;
      //
      // 2D/Dynamic structures.
      //
      ArrayRCP<ArrayRCP< LocalOrdinal> > lclInds2D_;
      ArrayRCP<ArrayRCP<GlobalOrdinal> > gblInds2D_;

      // This is deleted after fillCOmplete if "Optimize Storage" is set to \c true
      ArrayRCP<size_t>       numRowEntries_;

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

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

      bool haveRowInfo_;
      inline bool                     hasRowInfo() const {
#ifdef HAVE_TPETRA_DEBUG
        bool actuallyHasRowInfo = true;
        if (indicesAreAllocated() && getProfileType() == StaticProfile && rowPtrs_ == null) actuallyHasRowInfo = false;
        TEUCHOS_TEST_FOR_EXCEPTION(
            actuallyHasRowInfo != haveRowInfo_,
            std::logic_error, "Internal logic error. Please contact Tpetra team."
        )
#endif
        return haveRowInfo_;
      }

    bool insertGlobalIndicesWarnedEfficiency_;
    bool insertLocalIndicesWarnedEfficiency_;

  }; // class CrsGraph

} // namespace Tpetra

#endif