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


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

#include "KokkosCompat_ClassicNodeAPI_Wrapper.hpp"
#include "Kokkos_DualView.hpp"
#include "Kokkos_StaticCrsGraph.hpp"

#include "Teuchos_Describable.hpp"
#include "Teuchos_ParameterListAcceptorDefaultBase.hpp"


namespace Tpetra {

#ifndef DOXYGEN_SHOULD_SKIP_THIS
  //
  // Dear users: These are just forward declarations.  Please skip
  // over them and go down to the CrsMatrix class declaration.  Thank
  // you.
  //
  template <class LO, class GO, class N, const bool isClassic>
  class CrsGraph;

  // forward declaration (needed for "friend" inside CrsGraph)
  template <class S, class LO, class GO, class N, const bool isClassic>
  class CrsMatrix;

  namespace Experimental {
    // forward declaration (needed for "friend" inside CrsGraph)
    template<class S, class LO, class GO, class N>
    class BlockCrsMatrix;
  } // namespace Experimental

  namespace Details {
    // Forward declaration of an implementation detail of CrsGraph::clone.
    template<class OutputCrsGraphType, class InputCrsGraphType>
    class CrsGraphCopier {
    public:
      static Teuchos::RCP<OutputCrsGraphType>
      clone (const InputCrsGraphType& graphIn,
             const Teuchos::RCP<typename OutputCrsGraphType::node_type> nodeOut,
             const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
    };
  } // namespace Details
#endif // DOXYGEN_SHOULD_SKIP_THIS

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

  enum ELocalGlobal {
    LocalIndices,
    GlobalIndices
  };

  namespace Details {
    enum EStorageStatus {
      STORAGE_2D, //<! 2-D storage
      STORAGE_1D_UNPACKED, //<! 1-D "unpacked" storage
      STORAGE_1D_PACKED, //<! 1-D "packed" storage
      STORAGE_UB //<! Invalid value; upper bound on enum values
    };
  } // namespace Details

  template <class LocalOrdinal = Details::DefaultTypes::local_ordinal_type,
            class GlobalOrdinal = Details::DefaultTypes::global_ordinal_type,
            class Node = Details::DefaultTypes::node_type,
            const bool classic = Node::classic>
  class CrsGraph :
    public RowGraph<LocalOrdinal, GlobalOrdinal, Node>,
    public DistObject<GlobalOrdinal,
                      LocalOrdinal,
                      GlobalOrdinal,
                      Node>,
    public Teuchos::ParameterListAcceptorDefaultBase
  {
    static_assert (! classic, "The 'classic' version of Tpetra was deprecated long ago, and has been removed.");

    template <class S, class LO, class GO, class N, const bool isClassic>
    friend class CrsMatrix;
    template <class LO2, class GO2, class N2, const bool isClassic>
    friend class CrsGraph;
    template <class S, class LO, class GO, class N>
    friend class ::Tpetra::Experimental::BlockCrsMatrix;

    typedef DistObject<GlobalOrdinal, LocalOrdinal, GlobalOrdinal, Node> dist_object_type;

  public:
    typedef LocalOrdinal local_ordinal_type;
    typedef GlobalOrdinal global_ordinal_type;
    typedef Node node_type;

    typedef typename Node::device_type device_type;
    typedef typename device_type::execution_space execution_space;

    typedef Kokkos::StaticCrsGraph<LocalOrdinal,
                                   Kokkos::LayoutLeft,
                                   execution_space> local_graph_type;
    typedef local_graph_type LocalStaticCrsGraphType TPETRA_DEPRECATED;

    typedef typename local_graph_type::row_map_type t_RowPtrs TPETRA_DEPRECATED;
    typedef typename local_graph_type::row_map_type::non_const_type t_RowPtrsNC TPETRA_DEPRECATED;
    typedef typename local_graph_type::entries_type::non_const_type t_LocalOrdinal_1D TPETRA_DEPRECATED;

    typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
    typedef Tpetra::Import<LocalOrdinal, GlobalOrdinal, Node> import_type;
    typedef Tpetra::Export<LocalOrdinal, GlobalOrdinal, Node> export_type;



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

    CrsGraph (const Teuchos::RCP<const map_type>& rowMap,
              const Kokkos::DualView<const size_t*, execution_space>& numEntPerRow,
              const ProfileType pftype = DynamicProfile,
              const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

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

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

    CrsGraph (const Teuchos::RCP<const map_type>& rowMap,
              const Teuchos::RCP<const map_type>& colMap,
              const Kokkos::DualView<const size_t*, execution_space>& numEntPerRow,
              ProfileType pftype = DynamicProfile,
              const Teuchos::RCP<Teuchos::ParameterList>& params = null);

    CrsGraph (const Teuchos::RCP<const map_type>& rowMap,
              const Teuchos::RCP<const map_type>& colMap,
              const Teuchos::ArrayRCP<const size_t>& numEntPerRow,
              ProfileType pftype = DynamicProfile,
              const Teuchos::RCP<Teuchos::ParameterList>& params = null);

    CrsGraph (const Teuchos::RCP<const map_type>& rowMap,
              const Teuchos::RCP<const map_type>& colMap,
              const typename local_graph_type::row_map_type& rowPointers,
              const typename local_graph_type::entries_type::non_const_type& columnIndices,
              const Teuchos::RCP<Teuchos::ParameterList>& params = null);

    CrsGraph (const Teuchos::RCP<const map_type>& rowMap,
              const Teuchos::RCP<const map_type>& colMap,
              const Teuchos::ArrayRCP<size_t> & rowPointers,
              const Teuchos::ArrayRCP<LocalOrdinal> & columnIndices,
              const Teuchos::RCP<Teuchos::ParameterList>& params = null);

    CrsGraph (const Teuchos::RCP<const map_type>& rowMap,
              const Teuchos::RCP<const map_type>& colMap,
              const local_graph_type& lclGraph,
              const Teuchos::RCP<Teuchos::ParameterList>& params);

    template<class Node2>
    Teuchos::RCP<CrsGraph<LocalOrdinal, GlobalOrdinal, Node2, Node2::classic> >
    clone (const Teuchos::RCP<Node2>& node2,
           const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const
    {
      typedef CrsGraph<LocalOrdinal, GlobalOrdinal, Node2, Node2::classic> output_crs_graph_type;
      typedef CrsGraph<LocalOrdinal, GlobalOrdinal, Node, classic> input_crs_graph_type;
      typedef Details::CrsGraphCopier<output_crs_graph_type, input_crs_graph_type> copier_type;
      return copier_type::clone (*this, node2, params);
    }

    virtual ~CrsGraph();



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

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



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


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


    void removeLocalIndices (LocalOrdinal localRow);



    void globalAssemble ();

    void resumeFill (const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

    void
    fillComplete (const Teuchos::RCP<const map_type> &domainMap,
                  const Teuchos::RCP<const map_type> &rangeMap,
                  const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

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

    void
    expertStaticFillComplete (const Teuchos::RCP<const map_type> & domainMap,
                              const Teuchos::RCP<const map_type> & rangeMap,
                              const Teuchos::RCP<const import_type> &importer=Teuchos::null,
                              const Teuchos::RCP<const export_type> &exporter=Teuchos::null,
                              const Teuchos::RCP<Teuchos::ParameterList> &params=Teuchos::null);


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

    Teuchos::RCP<node_type> getNode() const;

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

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

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

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

    Teuchos::RCP<const import_type> getImporter () const;

    Teuchos::RCP<const export_type> 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 Teuchos::ArrayView<GlobalOrdinal>& Indices,
                      size_t& NumIndices) const;

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


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


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



    std::string description() const;

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



    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<GlobalOrdinal> &exports,
                    const Teuchos::ArrayView<size_t> & numPacketsPerLID,
                    size_t& constantNumPackets,
                    Distributor &distor);

    virtual void
    pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
          Teuchos::Array<GlobalOrdinal>& exports,
          const Teuchos::ArrayView<size_t>& numPacketsPerLID,
          size_t& constantNumPackets,
          Distributor& distor) const;

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


    void
    getNumEntriesPerLocalRowUpperBound (Teuchos::ArrayRCP<const size_t>& boundPerLocalRow,
                                        size_t& boundForAllLocalRows,
                                        bool& boundSameForAllLocalRows) const;

    void
    setAllIndices (const typename local_graph_type::row_map_type& rowPointers,
                   const typename local_graph_type::entries_type::non_const_type& columnIndices);

    void
    setAllIndices (const Teuchos::ArrayRCP<size_t> & rowPointers,
                   const Teuchos::ArrayRCP<LocalOrdinal> & columnIndices);

    Teuchos::ArrayRCP<const size_t> getNodeRowPtrs () const;


    Teuchos::ArrayRCP<const LocalOrdinal> getNodePackedIndices() const;

    void replaceColMap (const Teuchos::RCP<const map_type>& newColMap);

    void
    reindexColumns (const Teuchos::RCP<const map_type>& newColMap,
                    const Teuchos::RCP<const import_type>& newImport = Teuchos::null,
                    const bool sortIndicesInEachRow = true);

    void
    replaceDomainMapAndImporter (const Teuchos::RCP<const map_type>& newDomainMap,
                                 const Teuchos::RCP<const import_type>& newImporter);

    virtual void
    removeEmptyProcessesInPlace (const Teuchos::RCP<const map_type>& newMap);

    template<class ViewType, class OffsetViewType >
    struct pack_functor {
      typedef typename ViewType::execution_space execution_space;
      ViewType src;
      ViewType dest;
      OffsetViewType src_offset;
      OffsetViewType dest_offset;
      typedef typename OffsetViewType::non_const_value_type ScalarIndx;

      pack_functor(ViewType dest_, ViewType src_, OffsetViewType dest_offset_, OffsetViewType src_offset_):
        src(src_),dest(dest_),src_offset(src_offset_),dest_offset(dest_offset_) {};

      KOKKOS_INLINE_FUNCTION
      void operator() (size_t row) const {
        ScalarIndx i = src_offset(row);
        ScalarIndx j = dest_offset(row);
        const ScalarIndx k = dest_offset(row+1);
        for(;j<k;j++,i++) {
          dest(j) = src(i);
        }
      }
    };

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

    bool indicesAreAllocated () const;
    void allocateIndices (const ELocalGlobal lg);

    template <class T>
    Teuchos::ArrayRCP<Teuchos::Array<T> > allocateValues2D () const;

    template <class T>
    RowInfo updateLocalAllocAndValues (const RowInfo rowInfo,
                                       const size_t newAllocSize,
                                       Teuchos::Array<T>& rowVals)
    {
#ifdef HAVE_TPETRA_DEBUG
      TEUCHOS_TEST_FOR_EXCEPT( ! isLocallyIndexed () );
      TEUCHOS_TEST_FOR_EXCEPT( ! indicesAreAllocated() );
      TEUCHOS_TEST_FOR_EXCEPT( newAllocSize == 0 );
      TEUCHOS_TEST_FOR_EXCEPT( newAllocSize < rowInfo.allocSize );
      TEUCHOS_TEST_FOR_EXCEPT( ! rowMap_->isNodeLocalElement (rowInfo.localRow) );
#endif // HAVE_TPETRA_DEBUG

      // Teuchos::ArrayRCP::resize automatically copies over values on reallocation.
      lclInds2D_[rowInfo.localRow].resize (newAllocSize);
      rowVals.resize (newAllocSize);
      nodeNumAllocated_ += (newAllocSize - rowInfo.allocSize);

      RowInfo rowInfoOut = rowInfo;
      rowInfoOut.allocSize = newAllocSize;
      return rowInfoOut;
    }

    template <class T>
    RowInfo
    updateGlobalAllocAndValues (const RowInfo rowInfo,
                                const size_t newAllocSize,
                                Teuchos::Array<T>& rowVals)
    {
#ifdef HAVE_TPETRA_DEBUG
      TEUCHOS_TEST_FOR_EXCEPT( ! isGloballyIndexed () );
      TEUCHOS_TEST_FOR_EXCEPT( ! indicesAreAllocated () );
      TEUCHOS_TEST_FOR_EXCEPT( newAllocSize == 0 );
      TEUCHOS_TEST_FOR_EXCEPT( newAllocSize < rowInfo.allocSize );
      TEUCHOS_TEST_FOR_EXCEPT( ! rowMap_->isNodeLocalElement (rowInfo.localRow) );
#endif // HAVE_TPETRA_DEBUG

      // Teuchos::ArrayRCP::resize automatically copies over values on reallocation.
      gblInds2D_[rowInfo.localRow].resize (newAllocSize);
      rowVals.resize (newAllocSize);
      nodeNumAllocated_ += (newAllocSize - rowInfo.allocSize);

      RowInfo rowInfoOut = rowInfo;
      rowInfoOut.allocSize = newAllocSize;
      return rowInfoOut;
    }



    void makeColMap ();
    void makeIndicesLocal ();
    void makeImportExport ();



    template<ELocalGlobal lg>
    size_t filterIndices (const SLocalGlobalNCViews& inds) const
    {
      using Teuchos::ArrayView;
      static_assert (lg == GlobalIndices || lg == LocalIndices,
                     "Tpetra::CrsGraph::filterIndices: The template parameter "
                     "lg must be either GlobalIndices or LocalIndicies.");

      const map_type& cmap = *colMap_;
      size_t numFiltered = 0;
#ifdef HAVE_TPETRA_DEBUG
      size_t numFiltered_debug = 0;
#endif
      if (lg == GlobalIndices) {
        ArrayView<GlobalOrdinal> ginds = inds.ginds;
        typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin();
        typename ArrayView<GlobalOrdinal>::iterator cptr = ginds.begin();
        while (cptr != ginds.end()) {
          if (cmap.isNodeGlobalElement(*cptr)) {
            *fend++ = *cptr;
#ifdef HAVE_TPETRA_DEBUG
            ++numFiltered_debug;
#endif
          }
          ++cptr;
        }
        numFiltered = fend - ginds.begin();
      }
      else if (lg == LocalIndices) {
        ArrayView<LocalOrdinal> linds = inds.linds;
        typename ArrayView<LocalOrdinal>::iterator fend = linds.begin();
        typename ArrayView<LocalOrdinal>::iterator cptr = linds.begin();
        while (cptr != linds.end()) {
          if (cmap.isNodeLocalElement(*cptr)) {
            *fend++ = *cptr;
#ifdef HAVE_TPETRA_DEBUG
            ++numFiltered_debug;
#endif
          }
          ++cptr;
        }
        numFiltered = fend - linds.begin();
      }
#ifdef HAVE_TPETRA_DEBUG
      TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
#endif
      return numFiltered;
    }


    template<class T>
    size_t
    filterGlobalIndicesAndValues (const Teuchos::ArrayView<GlobalOrdinal>& ginds,
                                  const Teuchos::ArrayView<T>& vals) const
    {
      using Teuchos::ArrayView;
      const map_type& cmap = *colMap_;
      size_t numFiltered = 0;
      typename ArrayView<T>::iterator fvalsend = vals.begin();
      typename ArrayView<T>::iterator valscptr = vals.begin();
#ifdef HAVE_TPETRA_DEBUG
      size_t numFiltered_debug = 0;
#endif
      typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin();
      typename ArrayView<GlobalOrdinal>::iterator cptr = ginds.begin();
      while (cptr != ginds.end()) {
        if (cmap.isNodeGlobalElement (*cptr)) {
          *fend++ = *cptr;
          *fvalsend++ = *valscptr;
#ifdef HAVE_TPETRA_DEBUG
          ++numFiltered_debug;
#endif
        }
        ++cptr;
        ++valscptr;
      }
      numFiltered = fend - ginds.begin();
#ifdef HAVE_TPETRA_DEBUG
      TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
      TEUCHOS_TEST_FOR_EXCEPT( valscptr != vals.end() );
      const size_t numFilteredActual =
        static_cast<size_t> (fvalsend - vals.begin ());
      TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFilteredActual );
#endif // HAVE_TPETRA_DEBUG
      return numFiltered;
    }

    template<class T>
    size_t
    filterLocalIndicesAndValues (const Teuchos::ArrayView<LocalOrdinal>& linds,
                                 const Teuchos::ArrayView<T>& vals) const
    {
      using Teuchos::ArrayView;
      const map_type& cmap = *colMap_;
      size_t numFiltered = 0;
      typename ArrayView<T>::iterator fvalsend = vals.begin();
      typename ArrayView<T>::iterator valscptr = vals.begin();
#ifdef HAVE_TPETRA_DEBUG
      size_t numFiltered_debug = 0;
#endif
      typename ArrayView<LocalOrdinal>::iterator fend = linds.begin();
      typename ArrayView<LocalOrdinal>::iterator cptr = linds.begin();
      while (cptr != linds.end()) {
        if (cmap.isNodeLocalElement (*cptr)) {
          *fend++ = *cptr;
          *fvalsend++ = *valscptr;
#ifdef HAVE_TPETRA_DEBUG
          ++numFiltered_debug;
#endif
        }
        ++cptr;
        ++valscptr;
      }
      numFiltered = fend - linds.begin();
#ifdef HAVE_TPETRA_DEBUG
      TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
      TEUCHOS_TEST_FOR_EXCEPT( valscptr != vals.end() );
      const size_t numFilteredActual =
        Teuchos::as<size_t> (fvalsend - vals.begin ());
      TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFilteredActual );
#endif
      return numFiltered;
    }

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

    template<class Scalar>
    void
    insertIndicesAndValues (const RowInfo& rowInfo,
                            const SLocalGlobalViews& newInds,
                            const Teuchos::ArrayView<Scalar>& oldRowVals,
                            const Teuchos::ArrayView<const Scalar>& newRowVals,
                            const ELocalGlobal lg,
                            const ELocalGlobal I);
    void
    insertGlobalIndicesImpl (const LocalOrdinal myRow,
                             const Teuchos::ArrayView<const GlobalOrdinal> &indices);
    void
    insertLocalIndicesImpl (const LocalOrdinal myRow,
                            const Teuchos::ArrayView<const LocalOrdinal> &indices);
    void
    insertLocalIndicesFiltered (const LocalOrdinal localRow,
                                const Teuchos::ArrayView<const LocalOrdinal> &indices);

    void
    insertGlobalIndicesFiltered (const GlobalOrdinal localRow,
                                 const Teuchos::ArrayView<const GlobalOrdinal> &indices);

    template<class Scalar, class BinaryFunction>
    LocalOrdinal
    transformLocalValues (RowInfo rowInfo,
                          const Teuchos::ArrayView<Scalar>& rowVals,
                          const Teuchos::ArrayView<const LocalOrdinal>& inds,
                          const Teuchos::ArrayView<const Scalar>& newVals,
                          BinaryFunction f) const
    {
      typedef typename Teuchos::ArrayView<Scalar>::size_type size_type;
      const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid ();
      const size_type numElts = inds.size ();
      size_t hint = 0; // Guess for the current index k into rowVals

      // Get a view of the column indices in the row.  This amortizes
      // the cost of getting the view over all the entries of inds.
      Teuchos::ArrayView<const LocalOrdinal> colInds = getLocalView (rowInfo);

      LocalOrdinal numValid = 0; // number of valid local column indices
      for (size_type j = 0; j < numElts; ++j) {
        const size_t k = findLocalIndex (rowInfo, inds[j], colInds, hint);
        if (k != STINV) {
          rowVals[k] = f (rowVals[k], newVals[j]); // use binary function f
          hint = k+1;
          ++numValid;
        }
      }
      return numValid;
    }

    template<class Scalar, class BinaryFunction>
    LocalOrdinal
    transformGlobalValues (RowInfo rowInfo,
                           const Teuchos::ArrayView<Scalar>& rowVals,
                           const Teuchos::ArrayView<const GlobalOrdinal>& inds,
                           const Teuchos::ArrayView<const Scalar>& newVals,
                           BinaryFunction f) const
    {
      typedef typename Teuchos::ArrayView<Scalar>::size_type size_type;
      const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid ();
      const size_type numElts = inds.size ();
      size_t hint = 0; // guess at the index's relative offset in the row

      LocalOrdinal numValid = 0; // number of valid local column indices
      for (size_type j = 0; j < numElts; ++j) {
        const size_t k = findGlobalIndex (rowInfo, inds[j], hint);
        if (k != STINV) {
          rowVals[k] = f (rowVals[k], newVals[j]); // use binary function f
          hint = k+1;
          numValid++;
        }
      }
      return numValid;
    }



    bool isMerged () const;

    void setLocallyModified ();

    void sortAllIndices ();

    void sortRowIndices (const RowInfo rowinfo);

    template <class Scalar>
    void sortRowIndicesAndValues (const RowInfo rowinfo,
                                  const Teuchos::ArrayView<Scalar>& values);

    void mergeAllIndices ();

    void mergeRowIndices (RowInfo rowinfo);

    template<class Scalar>
    void
    mergeRowIndicesAndValues (RowInfo rowinfo,
                              const Teuchos::ArrayView<Scalar>& rowValues);

    void
    setDomainRangeMaps (const Teuchos::RCP<const map_type>& domainMap,
                        const Teuchos::RCP<const map_type>& rangeMap);

    void staticAssertions() const;
    void clearGlobalConstants();
    void computeGlobalConstants();

    RowInfo getRowInfo (const size_t myRow) const;

    Teuchos::ArrayView<const LocalOrdinal>
    getLocalView (const RowInfo rowinfo) const;

    Teuchos::ArrayView<LocalOrdinal>
    getLocalViewNonConst (const RowInfo rowinfo);

    Teuchos::ArrayView<const GlobalOrdinal>
    getGlobalView (const RowInfo rowinfo) const;

    Teuchos::ArrayView<GlobalOrdinal>
    getGlobalViewNonConst (const RowInfo rowinfo);

    size_t
    findLocalIndex (RowInfo rowinfo,
                    LocalOrdinal ind,
                    size_t hint = 0) const;

    size_t
    findLocalIndex (RowInfo rowinfo,
                    LocalOrdinal ind,
                    Teuchos::ArrayView<const LocalOrdinal> colInds,
                    size_t hint = 0) const;

    size_t findGlobalIndex (RowInfo rowinfo, GlobalOrdinal ind, size_t hint = 0) const;

    local_graph_type getLocalGraph () const;

    TPETRA_DEPRECATED local_graph_type getLocalGraph_Kokkos () const;

    void fillLocalGraph (const Teuchos::RCP<Teuchos::ParameterList>& params);

    bool hasRowInfo () const;

    void checkInternalState () const;

    Teuchos::RCP<const map_type> rowMap_;
    Teuchos::RCP<const map_type> colMap_;
    Teuchos::RCP<const map_type> rangeMap_;
    Teuchos::RCP<const map_type> domainMap_;

    Teuchos::RCP<const import_type> importer_;

    Teuchos::RCP<const export_type> exporter_;

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

    ProfileType pftype_;

    Kokkos::DualView<const size_t*, Kokkos::LayoutLeft, execution_space> k_numAllocPerRow_;

    size_t numAllocForAllRows_;



    typename local_graph_type::entries_type::non_const_type k_lclInds1D_;

    typedef Kokkos::View<GlobalOrdinal*, execution_space> t_GlobalOrdinal_1D;

    t_GlobalOrdinal_1D k_gblInds1D_;

    typename local_graph_type::row_map_type::const_type k_rowPtrs_;



    Teuchos::ArrayRCP<Teuchos::Array<LocalOrdinal> > lclInds2D_;

    Teuchos::ArrayRCP<Teuchos::Array<GlobalOrdinal> > gblInds2D_;

    typedef Kokkos::DualView<size_t*, Kokkos::LayoutLeft, execution_space> t_numRowEntries_;

    t_numRowEntries_ k_numRowEntries_;


    Details::EStorageStatus storageStatus_;

    bool indicesAreAllocated_;
    bool indicesAreLocal_;
    bool indicesAreGlobal_;
    bool fillComplete_;

    bool lowerTriangular_;
    bool upperTriangular_;
    bool indicesAreSorted_;
    bool noRedundancies_;
    bool haveLocalConstants_;
    bool haveGlobalConstants_;

    std::map<GlobalOrdinal, std::vector<GlobalOrdinal> > nonlocals_;

    bool sortGhostsAssociatedWithEachProcessor_;

  }; // class CrsGraph

  template <class LocalOrdinal, class GlobalOrdinal, class Node, const bool classic = Node::classic>
  Teuchos::RCP<CrsGraph<LocalOrdinal, GlobalOrdinal, Node, classic> >
  createCrsGraph (const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &map,
                  size_t maxNumEntriesPerRow = 0,
                  const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
  {
    using Teuchos::rcp;
    typedef CrsGraph<LocalOrdinal, GlobalOrdinal, Node, classic> graph_type;
    return rcp (new graph_type (map, maxNumEntriesPerRow, DynamicProfile, params));
  }

  namespace Details {

    template<class LocalOrdinal,
             class GlobalOrdinal,
             class OutputNodeType,
             class InputNodeType>
    class CrsGraphCopier<CrsGraph<LocalOrdinal, GlobalOrdinal, OutputNodeType>,
                         CrsGraph<LocalOrdinal, GlobalOrdinal, InputNodeType> > {
    public:
      typedef CrsGraph<LocalOrdinal, GlobalOrdinal, InputNodeType> input_crs_graph_type;
      typedef CrsGraph<LocalOrdinal, GlobalOrdinal, OutputNodeType> output_crs_graph_type;

      static Teuchos::RCP<output_crs_graph_type>
      clone (const input_crs_graph_type& graphIn,
             const Teuchos::RCP<OutputNodeType> &nodeOut,
             const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
      {
        using Teuchos::arcp;
        using Teuchos::ArrayRCP;
        using Teuchos::ArrayView;
        using Teuchos::null;
        using Teuchos::outArg;
        using Teuchos::ParameterList;
        using Teuchos::parameterList;
        using Teuchos::RCP;
        using Teuchos::rcp;
        using Teuchos::REDUCE_MIN;
        using Teuchos::reduceAll;
        using Teuchos::sublist;
        using std::cerr;
        using std::endl;
        typedef LocalOrdinal LO;
        typedef GlobalOrdinal GO;
        typedef typename ArrayView<const GO>::size_type size_type;
        typedef ::Tpetra::Map<LO, GO, InputNodeType> input_map_type;
        typedef ::Tpetra::Map<LO, GO, OutputNodeType> output_map_type;
        const char prefix[] = "Tpetra::Details::CrsGraphCopier::clone: ";

        // Set parameters' default values.
        bool debug = false;
        bool fillCompleteClone = true;
        bool useLocalIndices = graphIn.hasColMap ();
        ProfileType pftype = StaticProfile;
        // If the user provided a ParameterList, get values from there.
        if (! params.is_null ()) {
          fillCompleteClone = params->get ("fillComplete clone", fillCompleteClone);
          useLocalIndices = params->get ("Locally indexed clone", useLocalIndices);
          if (params->get ("Static profile clone", true) == false) {
            pftype = DynamicProfile;
          }
          debug = params->get ("Debug", debug);
        }

        const Teuchos::Comm<int>& comm = * (graphIn.getRowMap ()->getComm ());
        const int myRank = comm.getRank ();

        TEUCHOS_TEST_FOR_EXCEPTION(
                                   ! graphIn.hasColMap () && useLocalIndices, std::runtime_error,
                                   prefix << "You asked clone() to use local indices (by setting the "
                                   "\"Locally indexed clone\" parameter to true), but the source graph "
                                   "does not yet have a column Map, so this is impossible.");

        if (debug) {
          std::ostringstream os;
          os << "Process " << myRank << ": Cloning row Map" << endl;
          cerr << os.str ();
        }

        RCP<const output_map_type> clonedRowMap =
          graphIn.getRowMap ()->template clone<OutputNodeType> (nodeOut);

        // Invoke the output graph's constructor, using the input graph's
        // upper bounds on the number of entries in each local row.
        RCP<output_crs_graph_type> clonedGraph; // returned by this function
        {
          ArrayRCP<const size_t> numEntriesPerRow;
          size_t numEntriesForAll = 0;
          bool boundSameForAllLocalRows = true;

          if (debug) {
            std::ostringstream os;
            os << "Process " << myRank << ": Getting per-row bounds" << endl;
            cerr << os.str ();
          }
          graphIn.getNumEntriesPerLocalRowUpperBound (numEntriesPerRow,
                                                      numEntriesForAll,
                                                      boundSameForAllLocalRows);
          if (debug) {
            std::ostringstream os;
            os << "Process " << myRank << ": numEntriesForAll = "
               << numEntriesForAll << endl;
            cerr << os.str ();
          }

          if (debug) {
            std::ostringstream os;
            os << "Process " << myRank << ": graphIn.getNodeMaxNumRowEntries() = "
               << graphIn.getNodeMaxNumRowEntries () << endl;
            cerr << os.str ();
          }

          RCP<ParameterList> graphparams;
          if (params.is_null ()) {
            graphparams = parameterList ("CrsGraph");
          } else {
            graphparams = sublist (params, "CrsGraph");
          }
          if (useLocalIndices) {
            RCP<const output_map_type> clonedColMap =
              graphIn.getColMap ()->template clone<OutputNodeType> (nodeOut);
            if (boundSameForAllLocalRows) {
              clonedGraph = rcp (new output_crs_graph_type (clonedRowMap, clonedColMap,
                                                            numEntriesForAll, pftype,
                                                            graphparams));
            } else {
              clonedGraph = rcp (new output_crs_graph_type (clonedRowMap, clonedColMap,
                                                            numEntriesPerRow, pftype,
                                                            graphparams));
            }
          } else {
            if (boundSameForAllLocalRows) {
              clonedGraph = rcp (new output_crs_graph_type (clonedRowMap,
                                                            numEntriesForAll, pftype,
                                                            graphparams));
            } else {
              clonedGraph = rcp (new output_crs_graph_type (clonedRowMap,
                                                            numEntriesPerRow,
                                                            pftype, graphparams));
            }
          }

          if (debug) {
            std::ostringstream os;
            os << "Process " << myRank << ": Invoked output graph's constructor" << endl;
            cerr << os.str ();
          }

          // done with these
          numEntriesPerRow = null;
          numEntriesForAll = 0;
        }

        const input_map_type& inputRowMap = * (graphIn.getRowMap ());
        const size_type numRows =
          static_cast<size_type> (inputRowMap.getNodeNumElements ());

        bool failed = false;

        if (useLocalIndices) {
          const LO localMinLID = inputRowMap.getMinLocalIndex ();
          const LO localMaxLID = inputRowMap.getMaxLocalIndex ();

          if (graphIn.isLocallyIndexed ()) {
            if (numRows != 0) {
              try {
                ArrayView<const LO> linds;
                for (LO lrow = localMinLID; lrow <= localMaxLID; ++lrow) {
                  graphIn.getLocalRowView (lrow, linds);
                  if (linds.size () != 0) {
                    clonedGraph->insertLocalIndices (lrow, linds);
                  }
                }
              }
              catch (std::exception& e) {
                std::ostringstream os;
                os << "Process " << myRank << ": copying (reading local by view, "
                  "writing local) indices into the output graph threw an "
                  "exception: " << e.what () << endl;
                cerr << os.str ();
                failed = true;
              }
            }
          }
          else { // graphIn.isGloballyIndexed()
            TEUCHOS_TEST_FOR_EXCEPTION(
                                       ! graphIn.hasColMap () && useLocalIndices, std::invalid_argument,
                                       prefix << "You asked clone() to use local indices (by setting the "
                                       "\"Locally indexed clone\" parameter to true), but the source graph "
                                       "does not yet have a column Map, so this is impossible.");

            // The input graph has a column Map, but is globally indexed.
            // That's a bit weird, but we'll run with it.  In this case,
            // getLocalRowView won't work, but getLocalRowCopy should
            // still work; it will just have to convert from global to
            // local indices internally.

            try {
              // Make space for getLocalRowCopy to put column indices.
              //
              // This is only a hint; we may have to resize in the loop
              // below.  getNodeMaxNumRowEntries() may return nonsense if
              // fill is active.  The key bool in CrsGraph is
              // haveLocalConstants_.
              size_t myMaxNumRowEntries =
                graphIn.isFillActive () ? static_cast<size_t> (0) :
                graphIn.getNodeMaxNumRowEntries ();

              Array<LO> linds (myMaxNumRowEntries);

              // Copy each row into the new graph, using local indices.
              for (LO lrow = localMinLID; lrow <= localMaxLID; ++lrow) {
                size_t theNumEntries = graphIn.getNumEntriesInLocalRow (lrow);
                if (theNumEntries > myMaxNumRowEntries) {
                  myMaxNumRowEntries = theNumEntries;
                  linds.resize (myMaxNumRowEntries);
                }
                graphIn.getLocalRowCopy (lrow, linds (), theNumEntries);
                if (theNumEntries != 0) {
                  clonedGraph->insertLocalIndices (lrow, linds (0, theNumEntries));
                }
              }
            }
            catch (std::exception& e) {
              std::ostringstream os;
              os << "Process " << myRank << ": copying (reading local by copy, "
                "writing local) indices into the output graph threw an exception: "
                 << e.what () << endl;
              cerr << os.str ();
              failed = true;
            }
          }
        }
        else { /* useGlobalIndices */
          if (numRows != 0) {
            const GlobalOrdinal localMinGID = inputRowMap.getMinGlobalIndex ();
            const GlobalOrdinal localMaxGID = inputRowMap.getMaxGlobalIndex ();
            const bool inputRowMapIsContiguous = inputRowMap.isContiguous ();

            if (graphIn.isGloballyIndexed ()) {
              ArrayView<const GlobalOrdinal> ginds;

              if (inputRowMapIsContiguous) {
                try {
                  for (GO grow = localMinGID; grow <= localMaxGID; ++grow) {
                    graphIn.getGlobalRowView (grow, ginds);
                    if (ginds.size () != 0) {
                      clonedGraph->insertGlobalIndices (grow, ginds);
                    }
                  }
                }
                catch (std::exception& e) {
                  std::ostringstream os;
                  os << "Process " << myRank << ": copying (reading global by view, "
                    "writing global) indices into the output graph threw an "
                    "exception: " << e.what () << endl;
                  cerr << os.str ();
                  failed = true;
                }
              }
              else { // input row Map is not contiguous
                try {
                  ArrayView<const GO> inputRowMapGIDs = inputRowMap.getNodeElementList ();
                  for (size_type k = 0; k < numRows; ++k) {
                    const GO grow = inputRowMapGIDs[k];
                    graphIn.getGlobalRowView (grow, ginds);
                    if (ginds.size () != 0) {
                      clonedGraph->insertGlobalIndices (grow, ginds);
                    }
                  }
                }
                catch (std::exception& e) {
                  std::ostringstream os;
                  os << "Process " << myRank << ": copying (reading global by view, "
                    "writing global) indices into the output graph threw an "
                    "exception: " << e.what () << endl;
                  cerr << os.str ();
                  failed = true;
                }
              }
            }
            else { // graphIn.isLocallyIndexed()
              // Make space for getGlobalRowCopy to put column indices.
              //
              // This is only a hint; we may have to resize in the loop
              // below.  getNodeMaxNumRowEntries() may return nonsense if
              // fill is active.  The key bool in CrsGraph is
              // haveLocalConstants_.
              size_t myMaxNumRowEntries =
                graphIn.isFillActive () ? static_cast<size_t> (0) :
                graphIn.getNodeMaxNumRowEntries ();

              Array<GO> ginds (myMaxNumRowEntries);

              if (inputRowMapIsContiguous) {
                try {
                  for (GO grow = localMinGID; grow <= localMaxGID; ++grow) {
                    size_t theNumEntries = graphIn.getNumEntriesInGlobalRow (grow);
                    if (theNumEntries > myMaxNumRowEntries) {
                      myMaxNumRowEntries = theNumEntries;
                      ginds.resize (myMaxNumRowEntries);
                    }
                    graphIn.getGlobalRowCopy (grow, ginds (), theNumEntries);
                    if (theNumEntries != 0) {
                      clonedGraph->insertGlobalIndices (grow, ginds (0, theNumEntries));
                    }
                  }
                }
                catch (std::exception& e) {
                  std::ostringstream os;
                  os << "Process " << myRank << ": copying (reading global by copy, "
                    "writing global) indices into the output graph threw an "
                    "exception: " << e.what () << endl;
                  cerr << os.str ();
                  failed = true;
                }
              }
              else { // input row Map is not contiguous
                try {
                  ArrayView<const GO> inputRowMapGIDs = inputRowMap.getNodeElementList ();
                  for (size_type k = 0; k < numRows; ++k) {
                    const GO grow = inputRowMapGIDs[k];

                    size_t theNumEntries = graphIn.getNumEntriesInGlobalRow (grow);
                    if (theNumEntries > myMaxNumRowEntries) {
                      myMaxNumRowEntries = theNumEntries;
                      ginds.resize (myMaxNumRowEntries);
                    }
                    graphIn.getGlobalRowCopy (grow, ginds (), theNumEntries);
                    if (theNumEntries != 0) {
                      clonedGraph->insertGlobalIndices (grow, ginds (0, theNumEntries));
                    }
                  }
                }
                catch (std::exception& e) {
                  std::ostringstream os;
                  os << "Process " << myRank << ": copying (reading global by copy, "
                    "writing global) indices into the output graph threw an "
                    "exception: " << e.what () << endl;
                  cerr << os.str ();
                  failed = true;
                }
              }
            }
          } // numRows != 0
        }

        if (debug) {
          std::ostringstream os;
          os << "Process " << myRank << ": copied entries" << endl;
          cerr << os.str ();
        }

        if (fillCompleteClone) {
          RCP<ParameterList> fillparams = params.is_null () ?
            parameterList ("fillComplete") :
            sublist (params, "fillComplete");
          try {
            RCP<const output_map_type> clonedRangeMap;
            RCP<const output_map_type> clonedDomainMap;
            if (! graphIn.getRangeMap ().is_null () &&
                graphIn.getRangeMap () != graphIn.getRowMap ()) {
              clonedRangeMap =
                graphIn.getRangeMap ()->template clone<OutputNodeType> (nodeOut);
            }
            else {
              clonedRangeMap = clonedRowMap;
            }
            if (! graphIn.getDomainMap ().is_null ()
                && graphIn.getDomainMap () != graphIn.getRowMap ()) {
              clonedDomainMap =
                graphIn.getDomainMap ()->template clone<OutputNodeType> (nodeOut);
            }
            else {
              clonedDomainMap = clonedRowMap;
            }

            if (debug) {
              std::ostringstream os;
              os << "Process " << myRank << ": About to call fillComplete on "
                "cloned graph" << endl;
              cerr << os.str ();
            }
            clonedGraph->fillComplete (clonedDomainMap, clonedRangeMap, fillparams);
          }
          catch (std::exception &e) {
            failed = true;
            std::ostringstream os;
            os << prefix << "Process " << myRank << ": Caught the following "
              "exception while calling fillComplete() on clone of type"
               << endl << Teuchos::typeName (*clonedGraph) << endl;
            cerr << os.str ();
          }
        }

        int lclSuccess = failed ? 0 : 1;
        int gblSuccess = 1;
        reduceAll<int, int> (comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
        TEUCHOS_TEST_FOR_EXCEPTION(
                                   gblSuccess != 1, std::logic_error, prefix <<
                                   "Clone failed on at least one process.");

        if (debug) {
          std::ostringstream os;
          os << "Process " << myRank << ": Done with CrsGraph::clone" << endl;
          cerr << os.str ();
        }
        return clonedGraph;
      }
    };

  } // namespace Details
} // namespace Tpetra

#endif // TPETRA_CRSGRAPH_DECL_HPP