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_RowGraph.hpp"
#include "Tpetra_DistObject.hpp"
#include "Tpetra_Exceptions.hpp"
#include "Tpetra_Util.hpp" // need this here for sort2

#include "Kokkos_Sparse_findRelOffset.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

  namespace { // (anonymous)

    template<class ViewType>
    struct UnmanagedView {
      static_assert (Kokkos::is_view<ViewType>::value,
                     "ViewType must be a Kokkos::View specialization.");
      // FIXME (mfh 02 Dec 2015) Right now, this strips away other
      // memory traits.  Christian will add an "AllTraits" enum which is
      // the enum value of MemoryTraits<T>, that will help us fix this.
      typedef Kokkos::View<typename ViewType::data_type,
                           typename ViewType::array_layout,
                           typename ViewType::device_type,
                           Kokkos::MemoryUnmanaged> type;
    };

    template<class T, class BinaryFunction>
    T atomic_binary_function_update (volatile T* const dest, const T& inputVal, BinaryFunction f)
    {
      T oldVal = *dest;
      T assume;

      // NOTE (mfh 30 Nov 2015) I do NOT need a fence here for IBM
      // POWER architectures, because 'newval' depends on 'assume',
      // which depends on 'oldVal', which depends on '*dest'.  This
      // sets up a chain of read dependencies that should ensure
      // correct behavior given a sane memory model.
      do {
        assume = oldVal;
        T newVal = f (assume, inputVal);
        oldVal = Kokkos::atomic_compare_exchange (dest, assume, newVal);
      } while (assume != oldVal);

      return oldVal;
    }

  } // namespace (anonymous)
#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 = ::Tpetra::Details::DefaultTypes::local_ordinal_type,
            class GlobalOrdinal = ::Tpetra::Details::DefaultTypes::global_ordinal_type,
            class Node = ::Tpetra::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 = Teuchos::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 = Teuchos::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 = Teuchos::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 = Teuchos::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 = Teuchos::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 Teuchos::ParameterList> getValidParameters () const;



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

    void
    insertGlobalIndices (const GlobalOrdinal globalRow,
                         const LocalOrdinal numEnt,
                         const GlobalOrdinal inds[]);


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

    void
    insertLocalIndices (const LocalOrdinal localRow,
                        const LocalOrdinal numEnt,
                        const LocalOrdinal inds[]);


    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 Teuchos::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
    getLocalDiagOffsets (const Kokkos::View<size_t*, device_type, Kokkos::MemoryUnmanaged>& offsets) const;

    void
    getLocalDiagOffsets (Teuchos::ArrayRCP<size_t>& offsets) const;

    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
    {
      using Teuchos::arcp;
      using Teuchos::Array;
      using Teuchos::ArrayRCP;
      using Teuchos::null;
      const char tfecfFuncName[] = "allocateValues2D: ";

      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (! indicesAreAllocated (), std::runtime_error,
         "Graph indices must be allocated before values.");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (getProfileType () != DynamicProfile, std::runtime_error,
         "Graph indices must be allocated in a dynamic profile.");

      ArrayRCP<Array<T> > values2D;
      values2D = arcp<Array<T> > (getNodeNumRows ());
      if (lclInds2D_ != null) {
        const size_t numRows = lclInds2D_.size ();
        for (size_t r = 0; r < numRows; ++r) {
          values2D[r].resize (lclInds2D_[r].size ());
        }
      }
      else if (gblInds2D_ != null) {
        const size_t numRows = gblInds2D_.size ();
        for (size_t r = 0; r < numRows; ++r) {
          values2D[r].resize (gblInds2D_[r].size ());
        }
      }
      return values2D;
    }

    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)
    {
#ifdef HAVE_TPETRA_DEBUG
      const char tfecfFuncName[] = "insertIndicesAndValues: ";
#endif // HAVE_TPETRA_DEBUG

#ifdef HAVE_TPETRA_DEBUG
      size_t numNewInds = 0;
      try {
        numNewInds = insertIndices (rowInfo, newInds, lg, I);
      } catch (std::exception& e) {
        TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
          (true, std::runtime_error, "insertIndices threw an exception: "
           << e.what ());
      }
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (numNewInds > static_cast<size_t> (oldRowVals.size ()),
         std::runtime_error, "numNewInds (" << numNewInds << ") > "
         "oldRowVals.size() (" << oldRowVals.size () << ".");
#else
      const size_t numNewInds = insertIndices (rowInfo, newInds, lg, I);
#endif // HAVE_TPETRA_DEBUG

      typedef typename Teuchos::ArrayView<Scalar>::size_type size_type;

#ifdef HAVE_TPETRA_DEBUG
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (rowInfo.numEntries + numNewInds > static_cast<size_t> (oldRowVals.size ()),
         std::runtime_error, "rowInfo.numEntries (" << rowInfo.numEntries << ")"
         " + numNewInds (" << numNewInds << ") > oldRowVals.size() ("
         << oldRowVals.size () << ").");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (static_cast<size_type> (numNewInds) > newRowVals.size (),
         std::runtime_error, "numNewInds (" << numNewInds << ") > "
         "newRowVals.size() (" << newRowVals.size () << ").");
#endif // HAVE_TPETRA_DEBUG

      size_type oldInd = static_cast<size_type> (rowInfo.numEntries);

#ifdef HAVE_TPETRA_DEBUG
      try {
#endif // HAVE_TPETRA_DEBUG
        //NOTE: The code in the else branch fails on GCC 4.9 and newer in the assignement oldRowVals[oldInd] = newRowVals[newInd];
        //We supply a workaround n as well as other code variants which produce or not produce the error
#if defined(__GNUC__) && defined(__GNUC_MINOR__) && defined(__GNUC_PATCHLEVEL__)
#define GCC_VERSION __GNUC__*100+__GNUC_MINOR__*10+__GNUC_PATCHLEVEL__
#if GCC_VERSION >= 490
#define GCC_WORKAROUND
#endif
#endif
#ifdef GCC_WORKAROUND
        size_type nNI = static_cast<size_type>(numNewInds);
        if (nNI > 0)
          memcpy(&oldRowVals[oldInd], &newRowVals[0], nNI*sizeof(Scalar));
        /*
        //Original Code Fails
        for (size_type newInd = 0; newInd < static_cast<size_type> (numNewInds);
        ++newInd, ++oldInd) {
        oldRowVals[oldInd] = newRowVals[newInd];
        }

        //char cast variant fails
        char* oldRowValPtr = (char*)&oldRowVals[oldInd];
        const char* newRowValPtr = (const char*) &newRowVals[0];

        for(size_type newInd = 0; newInd < (nNI * sizeof(Scalar)); newInd++) {
        oldRowValPtr[newInd] = newRowValPtr[newInd];
        }

        //Raw ptr variant fails
        Scalar* oldRowValPtr = &oldRowVals[oldInd];
        Scalar* newRowValPtr = const_cast<Scalar*>(&newRowVals[0]);

        for(size_type newInd = 0; newInd < nNI; newInd++) {
        oldRowValPtr[newInd] = newRowValPtr[newInd];
        }

        //memcpy works
        for (size_type newInd = 0; newInd < nNI; newInd++) {
        memcpy( &oldRowVals[oldInd+newInd], &newRowVals[newInd], sizeof(Scalar));
        }

        //just one loop index fails
        for (size_type newInd = 0; newInd < nNI; newInd++) {
        oldRowVals[oldInd+newInd] = newRowVals[newInd];
        }

        //inline increment fails
        for (size_type newInd = 0; newInd < numNewInds;) {
        oldRowVals[oldInd++] = newRowVals[newInd++];
        }

        */

#else // GCC Workaround above
        for (size_type newInd = 0; newInd < static_cast<size_type> (numNewInds);
             ++newInd, ++oldInd) {
          oldRowVals[oldInd] = newRowVals[newInd];
        }
#endif // GCC Workaround
#ifdef HAVE_TPETRA_DEBUG
      }
      catch (std::exception& e) {
        TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
          (true, std::runtime_error, "for loop for copying values threw an "
           "exception: " << e.what ());
      }
#endif // HAVE_TPETRA_DEBUG
    }

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

    static const bool useAtomicUpdatesByDefault =
#ifdef KOKKOS_HAVE_SERIAL
      ! std::is_same<execution_space, Kokkos::Serial>::value;
#else
      true;
#endif // KOKKOS_HAVE_SERIAL

    template<class OutputScalarViewType,
             class LocalIndicesViewType,
             class InputScalarViewType,
             class BinaryFunction>
    LocalOrdinal
    transformLocalValues (const RowInfo& rowInfo,
                          const typename UnmanagedView<OutputScalarViewType>::type& rowVals,
                          const typename UnmanagedView<LocalIndicesViewType>::type& inds,
                          const typename UnmanagedView<InputScalarViewType>::type& newVals,
                          BinaryFunction f,
                          const bool atomic = useAtomicUpdatesByDefault) const
    {
      // We use static_assert here to check the template parameters,
      // rather than std::enable_if (e.g., on the return value, to
      // enable compilation only if the template parameters match the
      // desired attributes).  This turns obscure link errors into
      // clear compilation errors.  It also makes the return value a
      // lot easier to see.
      static_assert (Kokkos::is_view<OutputScalarViewType>::value,
                     "Template parameter OutputScalarViewType must be "
                     "a Kokkos::View.");
      static_assert (Kokkos::is_view<LocalIndicesViewType>::value,
                     "Template parameter LocalIndicesViewType must be "
                     "a Kokkos::View.");
      static_assert (Kokkos::is_view<InputScalarViewType>::value,
                     "Template parameter InputScalarViewType must be a "
                     "Kokkos::View.");
      static_assert (static_cast<int> (OutputScalarViewType::rank) == 1,
                     "Template parameter OutputScalarViewType must "
                     "have rank 1.");
      static_assert (static_cast<int> (LocalIndicesViewType::rank) == 1,
                     "Template parameter LocalIndicesViewType must "
                     "have rank 1.");
      static_assert (static_cast<int> (InputScalarViewType::rank) == 1,
                     "Template parameter InputScalarViewType must have "
                     "rank 1.");
      static_assert (std::is_same<
                       typename OutputScalarViewType::non_const_value_type,
                       typename InputScalarViewType::non_const_value_type>::value,
                     "Template parameters OutputScalarViewType and "
                     "InputScalarViewType must contain values of the same "
                     "type.");
      static_assert (std::is_same<
                       typename LocalIndicesViewType::non_const_value_type,
                       local_ordinal_type>::value,
                     "Template parameter LocalIndicesViewType must "
                     "contain values of type local_ordinal_type.");

      typedef typename OutputScalarViewType::non_const_value_type ST;
      typedef LocalOrdinal LO;
      typedef GlobalOrdinal GO;

      if (newVals.dimension_0 () != inds.dimension_0 ()) {
        // The sizes of the input arrays must match.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }
      const LO numElts = static_cast<LO> (inds.dimension_0 ());
      const bool sorted = this->isSorted ();

      LO numValid = 0; // number of valid input column indices
      size_t hint = 0; // Guess for the current index k into rowVals

      if (isLocallyIndexed ()) {
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getLocalKokkosRowView (rowInfo);

        for (LO j = 0; j < numElts; ++j) {
          const LO lclColInd = inds(j);
          const size_t offset =
            KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                         lclColInd, hint, sorted);
          if (offset != rowInfo.numEntries) {
            if (atomic) {
              // NOTE (mfh 30 Nov 2015) The commented-out code is
              // wrong because another thread may have changed
              // rowVals(offset) between those two lines of code.
              //
              //const ST newVal = f (rowVals(offset), newVals(j));
              //Kokkos::atomic_assign (&rowVals(offset), newVal);

              volatile ST* const dest = &rowVals(offset);
              (void) atomic_binary_function_update (dest, newVals(j), f);
            }
            else {
              // use binary function f
              rowVals(offset) = f (rowVals(offset), newVals(j));
            }
            hint = offset + 1;
            ++numValid;
          }
        }
      }
      else if (isGloballyIndexed ()) {
        // NOTE (mfh 26 Nov 2015) Dereferencing an RCP or reading its
        // pointer does NOT change its reference count.  Thus, this
        // code is still thread safe.
        if (colMap_.is_null ()) {
          // NO input column indices are valid in this case.  Either
          // the column Map hasn't been set yet (so local indices
          // don't exist yet), or the calling process owns no graph
          // entries.
          return numValid;
        }
        const map_type& colMap = *colMap_;
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getGlobalKokkosRowView (rowInfo);

        const GO GINV = Teuchos::OrdinalTraits<GO>::invalid ();
        for (LO j = 0; j < numElts; ++j) {
          const GO gblColInd = colMap.getGlobalElement (inds(j));
          if (gblColInd != GINV) {
            const size_t offset =
              KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                           gblColInd, hint, sorted);
            if (offset != rowInfo.numEntries) {
              if (atomic) {
                // NOTE (mfh 30 Nov 2015) The commented-out code is
                // wrong because another thread may have changed
                // rowVals(offset) between those two lines of code.
                //
                //const ST newVal = f (rowVals(offset), newVals(j));
                //Kokkos::atomic_assign (&rowVals(offset), newVal);

                volatile ST* const dest = &rowVals(offset);
                (void) atomic_binary_function_update (dest, newVals(j), f);
              }
              else {
                // use binary function f
                rowVals(offset) = f (rowVals(offset), newVals(j));
              }
              hint = offset + 1;
              numValid++;
            }
          }
        }
      }
      // If the graph is neither locally nor globally indexed on the
      // calling process, that means the calling process has no graph
      // entries.  Thus, none of the input column indices are valid.

      return numValid;
    }

    template<class OutputScalarViewType,
             class LocalIndicesViewType,
             class InputScalarViewType>
    LocalOrdinal
    sumIntoLocalValues (const RowInfo& rowInfo,
                        const typename UnmanagedView<OutputScalarViewType>::type& rowVals,
                        const typename UnmanagedView<LocalIndicesViewType>::type& inds,
                        const typename UnmanagedView<InputScalarViewType>::type& newVals,
                        const bool atomic = useAtomicUpdatesByDefault) const
    {
      // We use static_assert here to check the template parameters,
      // rather than std::enable_if (e.g., on the return value, to
      // enable compilation only if the template parameters match the
      // desired attributes).  This turns obscure link errors into
      // clear compilation errors.  It also makes the return value a
      // lot easier to see.
      static_assert (Kokkos::is_view<OutputScalarViewType>::value,
                     "Template parameter OutputScalarViewType must be "
                     "a Kokkos::View.");
      static_assert (Kokkos::is_view<LocalIndicesViewType>::value,
                     "Template parameter LocalIndicesViewType must be "
                     "a Kokkos::View.");
      static_assert (Kokkos::is_view<InputScalarViewType>::value,
                     "Template parameter InputScalarViewType must be a "
                     "Kokkos::View.");
      static_assert (static_cast<int> (OutputScalarViewType::rank) == 1,
                     "Template parameter OutputScalarViewType must "
                     "have rank 1.");
      static_assert (static_cast<int> (LocalIndicesViewType::rank) == 1,
                     "Template parameter LocalIndicesViewType must "
                     "have rank 1.");
      static_assert (static_cast<int> (InputScalarViewType::rank) == 1,
                     "Template parameter InputScalarViewType must have "
                     "rank 1.");
      static_assert (std::is_same<
                       typename OutputScalarViewType::non_const_value_type,
                       typename InputScalarViewType::non_const_value_type>::value,
                     "Template parameters OutputScalarViewType and "
                     "InputScalarViewType must contain values of the same "
                     "type.");
      static_assert (std::is_same<
                       typename LocalIndicesViewType::non_const_value_type,
                       local_ordinal_type>::value,
                     "Template parameter LocalIndicesViewType must "
                     "contain values of type local_ordinal_type.");

      typedef LocalOrdinal LO;
      typedef GlobalOrdinal GO;

      // Don't call this->hasColMap(), because that changes RCP's
      // reference count, which is not thread safe.  Just
      // dereferencing an RCP or calling RCP::is_null() does not
      // change its reference count.
      if (colMap_.is_null ()) {
        // No such thing as local column indices without a column Map.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }
      else if (newVals.dimension_0 () != inds.dimension_0 ()) {
        // The dimensions of the input arrays must match.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }

      const bool sorted = this->isSorted ();

      size_t hint = 0; // Guess for the current index k into rowVals
      LO numValid = 0; // number of valid local column indices

      // NOTE (mfh 11 Oct 2015) This method assumes UVM.  More
      // accurately, it assumes that the host execution space can
      // access data in both InputMemorySpace and ValsMemorySpace.

      if (isLocallyIndexed ()) {
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getLocalKokkosRowView (rowInfo);

        const LO numElts = static_cast<LO> (inds.dimension_0 ());
        for (LO j = 0; j < numElts; ++j) {
          const LO lclColInd = inds(j);
          const size_t offset =
            KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                         lclColInd, hint, sorted);
          if (offset != rowInfo.numEntries) {
            if (atomic) {
              Kokkos::atomic_add (&rowVals(offset), newVals(j));
            }
            else {
              rowVals(offset) += newVals(j);
            }
            hint = offset + 1;
            ++numValid;
          }
        }
      }
      else if (isGloballyIndexed ()) {
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getGlobalKokkosRowView (rowInfo);

        const LO numElts = static_cast<LO> (inds.dimension_0 ());
        for (LO j = 0; j < numElts; ++j) {
          const GO gblColInd = this->colMap_->getGlobalElement (inds(j));
          if (gblColInd != Teuchos::OrdinalTraits<GO>::invalid ()) {
            const size_t offset =
              KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                           gblColInd, hint, sorted);
            if (offset != rowInfo.numEntries) {
              if (atomic) {
                Kokkos::atomic_add (&rowVals(offset), newVals(j));
              }
              else {
                rowVals(offset) += newVals(j);
              }
              hint = offset + 1;
              ++numValid;
            }
          }
        }
      }
      // NOTE (mfh 26 Jun 2014, 26 Nov 2015) In the current version of
      // CrsGraph and CrsMatrix, it's possible for a matrix (or graph)
      // to be neither locally nor globally indexed on a process.
      // This means that the graph or matrix has no entries on that
      // process.  Epetra also works like this.  It's related to lazy
      // allocation (on first insertion, not at graph / matrix
      // construction).  Lazy allocation will go away because it is
      // not thread scalable.

      return numValid;
    }

    template<class OutputScalarViewType,
             class LocalIndicesViewType,
             class InputScalarViewType>
    LocalOrdinal
    replaceLocalValues (const RowInfo& rowInfo,
                        const typename UnmanagedView<OutputScalarViewType>::type& rowVals,
                        const typename UnmanagedView<LocalIndicesViewType>::type& inds,
                        const typename UnmanagedView<InputScalarViewType>::type& newVals) const
    {
      // We use static_assert here to check the template parameters,
      // rather than std::enable_if (e.g., on the return value, to
      // enable compilation only if the template parameters match the
      // desired attributes).  This turns obscure link errors into
      // clear compilation errors.  It also makes the return value a
      // lot easier to see.
      static_assert (Kokkos::is_view<OutputScalarViewType>::value,
                     "Template parameter OutputScalarViewType must be "
                     "a Kokkos::View.");
      static_assert (Kokkos::is_view<LocalIndicesViewType>::value,
                     "Template parameter LocalIndicesViewType must be "
                     "a Kokkos::View.");
      static_assert (Kokkos::is_view<InputScalarViewType>::value,
                     "Template parameter InputScalarViewType must be a "
                     "Kokkos::View.");
      static_assert (static_cast<int> (OutputScalarViewType::rank) == 1,
                     "Template parameter OutputScalarViewType must "
                     "have rank 1.");
      static_assert (static_cast<int> (LocalIndicesViewType::rank) == 1,
                     "Template parameter LocalIndicesViewType must "
                     "have rank 1.");
      static_assert (static_cast<int> (InputScalarViewType::rank) == 1,
                     "Template parameter InputScalarViewType must have "
                     "rank 1.");
      static_assert (std::is_same<
                       typename OutputScalarViewType::non_const_value_type,
                       typename InputScalarViewType::non_const_value_type>::value,
                     "Template parameters OutputScalarViewType and "
                     "InputScalarViewType must contain values of the same "
                     "type.");
      static_assert (std::is_same<
                       typename LocalIndicesViewType::non_const_value_type,
                       local_ordinal_type>::value,
                     "Template parameter LocalIndicesViewType must "
                     "contain values of type local_ordinal_type.");

      typedef LocalOrdinal LO;
      typedef GlobalOrdinal GO;

      // Don't call this->hasColMap(), because that changes RCP's
      // reference count, which is not thread safe.  Just
      // dereferencing an RCP or calling RCP::is_null() does not
      // change its reference count.
      if (colMap_.is_null ()) {
        // No such thing as local column indices without a column Map.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }
      else if (newVals.dimension_0 () != inds.dimension_0 ()) {
        // The dimensions of the input arrays must match.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }

      const bool sorted = this->isSorted ();

      size_t hint = 0; // Guess for the current index k into rowVals
      LO numValid = 0; // number of valid local column indices

      // NOTE (mfh 11 Oct 2015) This method assumes UVM.  More
      // accurately, it assumes that the host execution space can
      // access data in all the Views.

      if (isLocallyIndexed ()) {
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getLocalKokkosRowView (rowInfo);

        const LO numElts = static_cast<LO> (inds.dimension_0 ());
        for (LO j = 0; j < numElts; ++j) {
          const LO lclColInd = inds(j);
          const size_t offset =
            KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                         lclColInd, hint, sorted);
          if (offset != rowInfo.numEntries) {
            rowVals(offset) = newVals(j);
            hint = offset + 1;
            ++numValid;
          }
        }
      }
      else if (isGloballyIndexed ()) {
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getGlobalKokkosRowView (rowInfo);

        const LO numElts = static_cast<LO> (inds.dimension_0 ());
        for (LO j = 0; j < numElts; ++j) {
          const GO gblColInd = this->colMap_->getGlobalElement (inds(j));
          if (gblColInd != Teuchos::OrdinalTraits<GO>::invalid ()) {
            const size_t offset =
              KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                           gblColInd, hint, sorted);
            if (offset != rowInfo.numEntries) {
              rowVals(offset) = newVals(j);
              hint = offset + 1;
              ++numValid;
            }
          }
        }
      }
      // NOTE (mfh 26 Jun 2014, 26 Nov 2015) In the current version of
      // CrsGraph and CrsMatrix, it's possible for a matrix (or graph)
      // to be neither locally nor globally indexed on a process.
      // This means that the graph or matrix has no entries on that
      // process.  Epetra also works like this.  It's related to lazy
      // allocation (on first insertion, not at graph / matrix
      // construction).  Lazy allocation will go away because it is
      // not thread scalable.

      return numValid;
    }

    template<class Scalar, class InputMemorySpace, class ValsMemorySpace>
    LocalOrdinal
    sumIntoGlobalValues (const RowInfo& rowInfo,
                         const Kokkos::View<Scalar*, ValsMemorySpace,
                           Kokkos::MemoryUnmanaged>& rowVals,
                         const Kokkos::View<const GlobalOrdinal*, InputMemorySpace,
                           Kokkos::MemoryUnmanaged>& inds,
                         const Kokkos::View<const Scalar*, InputMemorySpace,
                           Kokkos::MemoryUnmanaged>& newVals,
                         const bool atomic = useAtomicUpdatesByDefault) const
    {
      typedef LocalOrdinal LO;
      typedef GlobalOrdinal GO;

      if (newVals.dimension_0 () != inds.dimension_0 ()) {
        // The dimensions of the input arrays must match.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }

      const bool sorted = this->isSorted ();

      size_t hint = 0; // guess at the index's relative offset in the row
      LO numValid = 0; // number of valid input column indices

      // NOTE (mfh 11 Oct 2015) This method assumes UVM.  More
      // accurately, it assumes that the host execution space can
      // access data in both InputMemorySpace and ValsMemorySpace.

      if (isLocallyIndexed ()) {
        // NOTE (mfh 04 Nov 2015) Dereferencing an RCP or reading its
        // pointer does NOT change its reference count.  Thus, this
        // code is still thread safe.
        if (colMap_.is_null ()) {
          // NO input column indices are valid in this case, since if
          // the column Map is null on the calling process, then the
          // calling process owns no graph entries.
          return numValid;
        }
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getLocalKokkosRowView (rowInfo);
        const LO LINV = Teuchos::OrdinalTraits<LO>::invalid ();

        const LO numElts = static_cast<LO> (inds.dimension_0 ());
        for (LO j = 0; j < numElts; ++j) {
          const LO lclColInd = this->colMap_->getLocalElement (inds(j));
          if (lclColInd != LINV) {
            const size_t offset =
              KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                           lclColInd, hint, sorted);
            if (offset != rowInfo.numEntries) {
              if (atomic) {
                Kokkos::atomic_add (&rowVals(offset), newVals(j));
              }
              else {
                rowVals(offset) += newVals(j);
              }
              hint = offset + 1;
              numValid++;
            }
          }
        }
      }
      else if (isGloballyIndexed ()) {
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getGlobalKokkosRowView (rowInfo);

        const LO numElts = static_cast<LO> (inds.dimension_0 ());
        for (LO j = 0; j < numElts; ++j) {
          const GO gblColInd = inds(j);
          const size_t offset =
            KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                         gblColInd, hint, sorted);
          if (offset != rowInfo.numEntries) {
            if (atomic) {
              Kokkos::atomic_add (&rowVals(offset), newVals(j));
            }
            else {
              rowVals(offset) += newVals(j);
            }
            hint = offset + 1;
            numValid++;
          }
        }
      }
      // If the graph is neither locally nor globally indexed on the
      // calling process, that means the calling process has no graph
      // entries.  Thus, none of the input column indices are valid.

      return numValid;
    }

    template<class OutputScalarViewType,
             class GlobalIndicesViewType,
             class InputScalarViewType>
    LocalOrdinal
    replaceGlobalValues (const RowInfo& rowInfo,
                         const typename UnmanagedView<OutputScalarViewType>::type& rowVals,
                         const typename UnmanagedView<GlobalIndicesViewType>::type& inds,
                         const typename UnmanagedView<InputScalarViewType>::type& newVals) const
    {
      // We use static_assert here to check the template parameters,
      // rather than std::enable_if (e.g., on the return value, to
      // enable compilation only if the template parameters match the
      // desired attributes).  This turns obscure link errors into
      // clear compilation errors.  It also makes the return value a
      // lot easier to see.
      static_assert (Kokkos::is_view<OutputScalarViewType>::value,
                     "Template parameter OutputScalarViewType must be "
                     "a Kokkos::View.");
      static_assert (Kokkos::is_view<GlobalIndicesViewType>::value,
                     "Template parameter GlobalIndicesViewType must be "
                     "a Kokkos::View.");
      static_assert (Kokkos::is_view<InputScalarViewType>::value,
                     "Template parameter InputScalarViewType must be a "
                     "Kokkos::View.");
      static_assert (static_cast<int> (OutputScalarViewType::rank) == 1,
                     "Template parameter OutputScalarViewType must "
                     "have rank 1.");
      static_assert (static_cast<int> (GlobalIndicesViewType::rank) == 1,
                     "Template parameter GlobalIndicesViewType must "
                     "have rank 1.");
      static_assert (static_cast<int> (InputScalarViewType::rank) == 1,
                     "Template parameter InputScalarViewType must have "
                     "rank 1.");
      static_assert (std::is_same<
                       typename OutputScalarViewType::non_const_value_type,
                       typename InputScalarViewType::non_const_value_type>::value,
                     "Template parameters OutputScalarViewType and "
                     "InputScalarViewType must contain values of the same "
                     "type.");
      static_assert (std::is_same<
                       typename GlobalIndicesViewType::non_const_value_type,
                       global_ordinal_type>::value,
                     "Template parameter GlobalIndicesViewType must "
                     "contain values of type global_ordinal_type.");

      typedef LocalOrdinal LO;
      typedef GlobalOrdinal GO;

      if (newVals.dimension_0 () != inds.dimension_0 ()) {
        // The dimensions of the input arrays must match.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }

      const bool sorted = this->isSorted ();

      size_t hint = 0; // guess at the index's relative offset in the row
      LO numValid = 0; // number of valid input column indices

      // NOTE (mfh 11 Oct 2015) This method assumes UVM.  More
      // accurately, it assumes that the host execution space can
      // access data in all the Views.

      if (isLocallyIndexed ()) {
        // NOTE (mfh 04 Nov 2015) Dereferencing an RCP or reading its
        // pointer does NOT change its reference count.  Thus, this
        // code is still thread safe.
        if (colMap_.is_null ()) {
          // NO input column indices are valid in this case, since if
          // the column Map is null on the calling process, then the
          // calling process owns no graph entries.
          return numValid;
        }
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getLocalKokkosRowView (rowInfo);
        const LO LINV = Teuchos::OrdinalTraits<LO>::invalid ();

        const LO numElts = static_cast<LO> (inds.dimension_0 ());
        for (LO j = 0; j < numElts; ++j) {
          const LO lclColInd = this->colMap_->getLocalElement (inds(j));
          if (lclColInd != LINV) {
            const size_t offset =
              KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                           lclColInd, hint, sorted);
            if (offset != rowInfo.numEntries) {
              rowVals(offset) = newVals(j);
              hint = offset + 1;
              numValid++;
            }
          }
        }
      }
      else if (isGloballyIndexed ()) {
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getGlobalKokkosRowView (rowInfo);

        const LO numElts = static_cast<LO> (inds.dimension_0 ());
        for (LO j = 0; j < numElts; ++j) {
          const GO gblColInd = inds(j);
          const size_t offset =
            KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                         gblColInd, hint, sorted);
          if (offset != rowInfo.numEntries) {
            rowVals(offset) = newVals(j);
            hint = offset + 1;
            numValid++;
          }
        }
      }
      // If the graph is neither locally nor globally indexed on the
      // calling process, that means the calling process has no graph
      // entries.  Thus, none of the input column indices are valid.

      return numValid;
    }

    template<class Scalar,
             class BinaryFunction,
             class InputMemorySpace,
             class ValsMemorySpace>
    LocalOrdinal
    transformGlobalValues (const RowInfo& rowInfo,
                           const Kokkos::View<Scalar*, ValsMemorySpace,
                             Kokkos::MemoryUnmanaged>& rowVals,
                           const Kokkos::View<const GlobalOrdinal*,
                             InputMemorySpace,
                             Kokkos::MemoryUnmanaged>& inds,
                           const Kokkos::View<const Scalar*,
                             InputMemorySpace,
                             Kokkos::MemoryUnmanaged>& newVals,
                           BinaryFunction f,
                           const bool atomic = useAtomicUpdatesByDefault) const
    {
      typedef LocalOrdinal LO;
      typedef GlobalOrdinal GO;

      if (newVals.dimension_0 () != inds.dimension_0 ()) {
        // The sizes of the input arrays must match.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }

      const LO numElts = static_cast<LO> (inds.dimension_0 ());
      const bool sorted = this->isSorted ();

      LO numValid = 0; // number of valid input column indices
      size_t hint = 0; // guess at the index's relative offset in the row

      if (isLocallyIndexed ()) {
        // NOTE (mfh 04 Nov 2015) Dereferencing an RCP or reading its
        // pointer does NOT change its reference count.  Thus, this
        // code is still thread safe.
        if (colMap_.is_null ()) {
          // NO input column indices are valid in this case, since if
          // the column Map is null on the calling process, then the
          // calling process owns no graph entries.
          return numValid;
        }
        const map_type& colMap = *colMap_;
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getLocalKokkosRowView (rowInfo);

        const LO LINV = Teuchos::OrdinalTraits<LO>::invalid ();
        for (LO j = 0; j < numElts; ++j) {
          const LO lclColInd = colMap.getLocalElement (inds(j));
          if (lclColInd != LINV) {
            const size_t offset =
              KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                           lclColInd, hint, sorted);
            if (offset != rowInfo.numEntries) {
              if (atomic) {
                // NOTE (mfh 30 Nov 2015) The commented-out code is
                // wrong because another thread may have changed
                // rowVals(offset) between those two lines of code.
                //
                //const Scalar newVal = f (rowVals(offset), newVals(j));
                //Kokkos::atomic_assign (&rowVals(offset), newVal);

                volatile Scalar* const dest = &rowVals(offset);
                (void) atomic_binary_function_update (dest, newVals(j), f);
              }
              else {
                // use binary function f
                rowVals(offset) = f (rowVals(offset), newVals(j));
              }
              hint = offset + 1;
              numValid++;
            }
          }
        }
      }
      else if (isGloballyIndexed ()) {
        // Get a view of the column indices in the row.  This amortizes
        // the cost of getting the view over all the entries of inds.
        auto colInds = this->getGlobalKokkosRowView (rowInfo);

        for (LO j = 0; j < numElts; ++j) {
          const GO gblColInd = inds(j);
          const size_t offset =
            KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
                                         gblColInd, hint, sorted);
          if (offset != rowInfo.numEntries) {
            if (atomic) {
              // NOTE (mfh 30 Nov 2015) The commented-out code is
              // wrong because another thread may have changed
              // rowVals(offset) between those two lines of code.
              //
              //const Scalar newVal = f (rowVals(offset), newVals(j));
              //Kokkos::atomic_assign (&rowVals(offset), newVal);

              volatile Scalar* const dest = &rowVals(offset);
              (void) atomic_binary_function_update (dest, newVals(j), f);
            }
            else {
              // use binary function f
              rowVals(offset) = f (rowVals(offset), newVals(j));
            }
            hint = offset + 1;
            numValid++;
          }
        }
      }
      // If the graph is neither locally nor globally indexed on the
      // calling process, that means the calling process has no graph
      // entries.  Thus, none of the input column indices are valid.

      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)
    {
      if (rowinfo.numEntries > 0) {
        Teuchos::ArrayView<LocalOrdinal> inds_view =
          this->getLocalViewNonConst (rowinfo);
        sort2 (inds_view.begin (), inds_view.begin () + rowinfo.numEntries,
               values.begin ());
      }
    }

    void mergeAllIndices ();

    void mergeRowIndices (RowInfo rowinfo);

    template<class Scalar>
    void
    mergeRowIndicesAndValues (RowInfo rowinfo,
                              const Teuchos::ArrayView<Scalar>& rowValues)
    {
      using Teuchos::ArrayView;
      const char tfecfFuncName[] = "mergeRowIndicesAndValues: ";
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (isStorageOptimized(), std::logic_error, "It is invalid to call this "
         "method if the graph's storage has already been optimized.  Please "
         "report this bug to the Tpetra developers.");

      typedef typename ArrayView<Scalar>::iterator Iter;
      Iter rowValueIter = rowValues.begin ();
      ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst (rowinfo);
      typename ArrayView<LocalOrdinal>::iterator beg, end, newend;

      // beg,end define a half-exclusive interval over which to iterate.
      beg = inds_view.begin();
      end = inds_view.begin() + rowinfo.numEntries;
      newend = beg;
      if (beg != end) {
        typename ArrayView<LocalOrdinal>::iterator cur = beg + 1;
        Iter vcur = rowValueIter + 1;
        Iter vend = rowValueIter;
        cur = beg+1;
        while (cur != end) {
          if (*cur != *newend) {
            // new entry; save it
            ++newend;
            ++vend;
            (*newend) = (*cur);
            (*vend) = (*vcur);
          }
          else {
            // old entry; merge it
            //(*vend) = f (*vend, *vcur);
            (*vend) += *vcur;
          }
          ++cur;
          ++vcur;
        }
        ++newend; // one past the last entry, per typical [beg,end) semantics
      }
      const size_t mergedEntries = newend - beg;
#ifdef HAVE_TPETRA_DEBUG
      // merge should not have eliminated any entries; if so, the
      // assignment below will destroy the packed structure
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (isStorageOptimized() && mergedEntries != rowinfo.numEntries,
         std::logic_error,
         "Merge was incorrect; it eliminated entries from the graph.  "
         << "Please report this bug to the Tpetra developers.");
#endif // HAVE_TPETRA_DEBUG

      k_numRowEntries_(rowinfo.localRow) = mergedEntries;
      nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries);
    }


    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 LocalOrdinal myRow) const;

    RowInfo getRowInfoFromGlobalRowIndex (const GlobalOrdinal gblRow) const;

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

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

    LocalOrdinal
    getLocalViewRawConst (const LocalOrdinal*& lclInds,
                          LocalOrdinal& numEnt,
                          const RowInfo& rowinfo) const;

  private:

    Kokkos::View<const LocalOrdinal*, execution_space, Kokkos::MemoryUnmanaged>
    getLocalKokkosRowView (const RowInfo& rowinfo) const;

    Kokkos::View<LocalOrdinal*, execution_space, Kokkos::MemoryUnmanaged>
    getLocalKokkosRowViewNonConst (const RowInfo& rowinfo);

    Kokkos::View<const GlobalOrdinal*, execution_space, Kokkos::MemoryUnmanaged>
    getGlobalKokkosRowView (const RowInfo& rowinfo) const;

  protected:

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

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

    LocalOrdinal
    getGlobalViewRawConst (const GlobalOrdinal*& gblInds,
                           LocalOrdinal& numEnt,
                           const RowInfo& rowinfo) const;

  public:
    local_graph_type getLocalGraph () const;

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

    typename Kokkos::View<const size_t*, execution_space>::HostMirror
    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_;

    typename Kokkos::View<size_t*, Kokkos::LayoutLeft, device_type>::HostMirror
      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::Array;
        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