Tpetra_CrsMatrix_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_CRSMATRIX_DECL_HPP
#define TPETRA_CRSMATRIX_DECL_HPP


#include "Tpetra_ConfigDefs.hpp"
#include "Tpetra_RowMatrix_decl.hpp"
#include "Tpetra_Exceptions.hpp"
#include "Tpetra_DistObject.hpp"
#include "Tpetra_CrsGraph.hpp"
#include "Tpetra_Vector.hpp"

#include "KokkosCompat_ClassicNodeAPI_Wrapper.hpp"
// localMultiply is templated on DomainScalar and RangeScalar, so we
// have to include this header file here, rather than in the _def
// header file, so that we can get KokkosSparse::spmv.
#include "Kokkos_Sparse.hpp"
// localGaussSeidel and reorderedLocalGaussSeidel are templated on
// DomainScalar and RangeScalar, so we have to include this header
// file here, rather than in the _def header file, so that we can get
// the interfaces to the corresponding local computational kernels.
#include "Kokkos_Sparse_impl_sor.hpp"

namespace Tpetra {
  template <class Scalar = Details::DefaultTypes::scalar_type,
            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 CrsMatrix :
    public RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>,
    public DistObject<char, LocalOrdinal, GlobalOrdinal, Node, classic>
  {
  public:


    typedef Scalar scalar_type;
    typedef typename Kokkos::Details::ArithTraits<Scalar>::val_type impl_scalar_type;
    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 typename Kokkos::Details::ArithTraits<impl_scalar_type>::mag_type mag_type;

    typedef Map<LocalOrdinal, GlobalOrdinal, Node> map_type;

    typedef Import<LocalOrdinal, GlobalOrdinal, Node> import_type;

    typedef Export<LocalOrdinal, GlobalOrdinal, Node> export_type;

    typedef CrsGraph<LocalOrdinal, GlobalOrdinal, Node, classic> crs_graph_type;

    typedef typename crs_graph_type::local_graph_type local_graph_type;

    typedef Kokkos::CrsMatrix<impl_scalar_type, LocalOrdinal, execution_space, void,
                              typename local_graph_type::size_type> local_matrix_type;

    typedef typename local_matrix_type::row_map_type t_RowPtrs TPETRA_DEPRECATED;
    typedef typename local_matrix_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 typename local_matrix_type::values_type t_ValuesType TPETRA_DEPRECATED;

    typedef local_matrix_type k_local_matrix_type TPETRA_DEPRECATED;



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

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

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

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

    explicit CrsMatrix (const Teuchos::RCP<const crs_graph_type>& graph,
                        const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);

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

    CrsMatrix (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::ArrayRCP<Scalar>& values,
               const Teuchos::RCP<ParameterList>& params = null);

    CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
               const Teuchos::RCP<const map_type>& colMap,
               const local_matrix_type& lclMatrix,
               const Teuchos::RCP<Teuchos::ParameterList>& params = null);

    // This friend declaration makes the clone() method work.
    template <class S2, class LO2, class GO2, class N2, const bool isClassic>
    friend class CrsMatrix;

    template <class Node2>
    Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2, Node2::classic> >
    clone (const Teuchos::RCP<Node2>& node2,
           const Teuchos::RCP<Teuchos::ParameterList>& params = null) const
    {
      using Teuchos::ArrayRCP;
      using Teuchos::null;
      using Teuchos::ParameterList;
      using Teuchos::RCP;
      using Teuchos::rcp;
      using Teuchos::sublist;
      typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2, Node2::classic> CrsMatrix2;
      typedef Map<LocalOrdinal, GlobalOrdinal, Node2> Map2;
      const char tfecfFuncName[] = "clone";

      // Get parameter values.  Set them initially to their default values.
      bool fillCompleteClone = true;
      bool useLocalIndices = this->hasColMap ();
      ProfileType pftype = StaticProfile;
      if (! params.is_null ()) {
        fillCompleteClone = params->get ("fillComplete clone", fillCompleteClone);
        useLocalIndices = params->get ("Locally indexed clone", useLocalIndices);

        bool staticProfileClone = true;
        staticProfileClone = params->get ("Static profile clone", staticProfileClone);
        pftype = staticProfileClone ? StaticProfile : DynamicProfile;
      }

      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
        ! this->hasColMap () && useLocalIndices, std::runtime_error,
        ": You requested that the returned clone have local indices, but the "
        "the source matrix does not have a column Map yet.");

      RCP<const Map2> clonedRowMap = this->getRowMap ()->template clone<Node2> (node2);

      // Get an upper bound on the number of entries per row.
      RCP<CrsMatrix2> clonedMatrix;
      ArrayRCP<const size_t> numEntriesPerRow;
      size_t numEntriesForAll = 0;
      bool boundSameForAllLocalRows = false;
      staticGraph_->getNumEntriesPerLocalRowUpperBound (numEntriesPerRow,
                                                        numEntriesForAll,
                                                        boundSameForAllLocalRows);
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
        numEntriesForAll != 0 &&
        static_cast<size_t> (numEntriesPerRow.size ()) != 0,
        std::logic_error, ": getNumEntriesPerLocalRowUpperBound returned a "
        "nonzero numEntriesForAll = " << numEntriesForAll << " , as well as a "
        "numEntriesPerRow array of nonzero length " << numEntriesPerRow.size ()
        << ".  This should never happen.  Please report this bug to the Tpetra "
        "developers.");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
        numEntriesForAll != 0 && ! boundSameForAllLocalRows,
        std::logic_error, ": getNumEntriesPerLocalRowUpperBound returned a "
        "nonzero numEntriesForAll = " << numEntriesForAll << " , but claims "
        "(via its third output value) that the upper bound is not the same for "
        "all rows.  This should never happen.  Please report this bug to the "
        "Tpetra developers.");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
        numEntriesPerRow.size () != 0 && boundSameForAllLocalRows,
        std::logic_error, ": getNumEntriesPerLocalRowUpperBound returned a "
        "numEntriesPerRow array of nonzero length " << numEntriesPerRow.size ()
        << ", but claims (via its third output value) that the upper bound is "
        "not the same for all rows.  This should never happen.  Please report "
        "this bug to the Tpetra developers.");

      RCP<ParameterList> matParams =
        params.is_null () ? null : sublist (params,"CrsMatrix");
      if (useLocalIndices) {
        RCP<const Map2> clonedColMap =
          this->getColMap ()->template clone<Node2> (node2);
        if (numEntriesPerRow.is_null ()) {
          clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap,
                                              numEntriesForAll, pftype,
                                              matParams));
        }
        else {
          clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap,
                                              numEntriesPerRow, pftype,
                                              matParams));
        }
      }
      else {
        if (numEntriesPerRow.is_null ()) {
          clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntriesForAll,
                                              pftype, matParams));
        }
        else {
          clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntriesPerRow,
                                              pftype, matParams));
        }
      }
      // done with these
      numEntriesPerRow = Teuchos::null;
      numEntriesForAll = 0;

      if (useLocalIndices) {
        clonedMatrix->allocateValues (LocalIndices,
                                      CrsMatrix2::GraphNotYetAllocated);
        if (this->isLocallyIndexed ()) {
          ArrayView<const LocalOrdinal> linds;
          ArrayView<const Scalar> vals;
          for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex ();
               lrow <= clonedRowMap->getMaxLocalIndex ();
               ++lrow) {
            this->getLocalRowView (lrow, linds, vals);
            if (linds.size ()) {
              clonedMatrix->insertLocalValues (lrow, linds, vals);
            }
          }
        }
        else { // this->isGloballyIndexed()
          Array<LocalOrdinal> linds;
          Array<Scalar> vals;
          for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex ();
               lrow <= clonedRowMap->getMaxLocalIndex ();
               ++lrow) {
            size_t theNumEntries = this->getNumEntriesInLocalRow (lrow);
            if (theNumEntries > static_cast<size_t> (linds.size ())) {
              linds.resize (theNumEntries);
            }
            if (theNumEntries > static_cast<size_t> (vals.size ())) {
              vals.resize (theNumEntries);
            }
            this->getLocalRowCopy (clonedRowMap->getGlobalElement (lrow),
                                   linds (), vals (), theNumEntries);
            if (theNumEntries != 0) {
              clonedMatrix->insertLocalValues (lrow, linds (0, theNumEntries),
                                               vals (0, theNumEntries));
            }
          }
        }
      }
      else { // useGlobalIndices
        clonedMatrix->allocateValues (GlobalIndices,
                                      CrsMatrix2::GraphNotYetAllocated);
        if (this->isGloballyIndexed ()) {
          ArrayView<const GlobalOrdinal> ginds;
          ArrayView<const Scalar> vals;
          for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex ();
               grow <= clonedRowMap->getMaxGlobalIndex ();
               ++grow) {
            this->getGlobalRowView (grow, ginds, vals);
            if (ginds.size () > 0) {
              clonedMatrix->insertGlobalValues (grow, ginds, vals);
            }
          }
        }
        else { // this->isLocallyIndexed()
          Array<GlobalOrdinal> ginds;
          Array<Scalar> vals;
          for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex ();
               grow <= clonedRowMap->getMaxGlobalIndex ();
               ++grow) {
            size_t theNumEntries = this->getNumEntriesInGlobalRow (grow);
            if (theNumEntries > static_cast<size_t> (ginds.size ())) {
              ginds.resize (theNumEntries);
            }
            if (theNumEntries > static_cast<size_t> (vals.size ())) {
              vals.resize (theNumEntries);
            }
            this->getGlobalRowCopy (grow, ginds (), vals (), theNumEntries);
            if (theNumEntries != 0) {
              clonedMatrix->insertGlobalValues (grow, ginds (0, theNumEntries),
                                                vals (0, theNumEntries));
            }
          }
        }
      }

      if (fillCompleteClone) {
        RCP<const Map2> clonedRangeMap;
        RCP<const Map2> clonedDomainMap;
        try {
          if (! this->getRangeMap ().is_null () &&
              this->getRangeMap () != clonedRowMap) {
            clonedRangeMap  = this->getRangeMap ()->template clone<Node2> (node2);
          }
          else {
            clonedRangeMap = clonedRowMap;
          }
          if (! this->getDomainMap ().is_null () &&
              this->getDomainMap () != clonedRowMap) {
            clonedDomainMap = this->getDomainMap ()->template clone<Node2> (node2);
          }
          else {
            clonedDomainMap = clonedRowMap;
          }
        }
        catch (std::exception &e) {
          const bool caughtExceptionOnClone = true;
          TEUCHOS_TEST_FOR_EXCEPTION
            (caughtExceptionOnClone, std::runtime_error,
             Teuchos::typeName (*this) << "::clone: Caught the following "
             "exception while cloning range and domain Maps on a clone of "
             "type " << Teuchos::typeName (*clonedMatrix) << ": " << e.what ());
        }

        RCP<ParameterList> fillparams =
          params.is_null () ? Teuchos::null : sublist (params, "fillComplete");
        try {
          clonedMatrix->fillComplete (clonedDomainMap, clonedRangeMap,
                                      fillparams);
        }
        catch (std::exception &e) {
          const bool caughtExceptionOnClone = true;
          TEUCHOS_TEST_FOR_EXCEPTION(
            caughtExceptionOnClone, std::runtime_error,
            Teuchos::typeName (*this) << "::clone: Caught the following "
            "exception while calling fillComplete() on a clone of type "
            << Teuchos::typeName (*clonedMatrix) << ": " << e.what ());
        }
      }
      return clonedMatrix;
    }

    virtual ~CrsMatrix ();



    //
    void
    insertGlobalValues (const GlobalOrdinal globalRow,
                        const Teuchos::ArrayView<const GlobalOrdinal>& cols,
                        const Teuchos::ArrayView<const Scalar>& vals);

    void
    insertLocalValues (const LocalOrdinal localRow,
                       const ArrayView<const LocalOrdinal> &cols,
                       const ArrayView<const Scalar> &vals);

    LocalOrdinal
    replaceGlobalValues (GlobalOrdinal globalRow,
                         const ArrayView<const GlobalOrdinal>& cols,
                         const ArrayView<const Scalar>& vals);

    LocalOrdinal
    replaceLocalValues (const LocalOrdinal localRow,
                        const ArrayView<const LocalOrdinal>& cols,
                        const ArrayView<const Scalar>& vals);

    LocalOrdinal
    sumIntoGlobalValues (const GlobalOrdinal globalRow,
                         const ArrayView<const GlobalOrdinal> &cols,
                         const ArrayView<const Scalar>        &vals);

    LocalOrdinal
    sumIntoLocalValues (const LocalOrdinal localRow,
                        const ArrayView<const LocalOrdinal>& cols,
                        const ArrayView<const Scalar>& vals);

    void setAllToScalar (const Scalar &alpha);

    void scale (const Scalar &alpha);


    void
    setAllValues (const typename local_matrix_type::row_map_type& rowPointers,
                  const typename local_graph_type::entries_type::non_const_type& columnIndices,
                  const typename local_matrix_type::values_type& values);


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

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



    void globalAssemble();

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

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

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

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

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

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

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

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



    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 RowGraph<LocalOrdinal, GlobalOrdinal, Node> > getGraph () const;

    Teuchos::RCP<const crs_graph_type> getCrsGraph () const;

    local_matrix_type getLocalMatrix () const {return lclMatrix_; }

    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;


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

    ProfileType getProfileType () const;

    bool isStaticGraph () const;

    mag_type getFrobeniusNorm () const;

    virtual bool supportsRowViews () const;

    void
    getGlobalRowCopy (GlobalOrdinal GlobalRow,
                      const Teuchos::ArrayView<GlobalOrdinal>& Indices,
                      const Teuchos::ArrayView<Scalar>& Values,
                      size_t& NumEntries) const;

    void
    getLocalRowCopy (LocalOrdinal localRow,
                     const Teuchos::ArrayView<LocalOrdinal>& colInds,
                     const Teuchos::ArrayView<Scalar>& vals,
                     size_t& numEntries) const;

    void
    getGlobalRowView (GlobalOrdinal GlobalRow,
                      Teuchos::ArrayView<const GlobalOrdinal>& indices,
                      Teuchos::ArrayView<const Scalar>& values) const;

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

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

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

    void
    getLocalDiagCopy (Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& diag,
                      const Teuchos::ArrayView<const size_t>& offsets) const;



    void
    leftScale (const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& x);

    void
    rightScale (const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& x);



    template <class DomainScalar, class RangeScalar>
    void
    localMultiply (const MultiVector<DomainScalar, LocalOrdinal, GlobalOrdinal, Node, classic>& X,
                   MultiVector<RangeScalar, LocalOrdinal, GlobalOrdinal, Node, classic>& Y,
                   Teuchos::ETransp mode,
                   RangeScalar alpha,
                   RangeScalar beta) const
    {
      using Teuchos::NO_TRANS;
      // Just like Scalar and impl_scalar_type may differ in CrsMatrix,
      // RangeScalar and its corresponding impl_scalar_type may differ in
      // MultiVector.
      typedef typename MultiVector<RangeScalar, LocalOrdinal, GlobalOrdinal,
        Node, classic>::impl_scalar_type range_impl_scalar_type;
#ifdef HAVE_TPETRA_DEBUG
      const char tfecfFuncName[] = "localMultiply: ";
#endif // HAVE_TPETRA_DEBUG

      const range_impl_scalar_type theAlpha = static_cast<range_impl_scalar_type> (alpha);
      const range_impl_scalar_type theBeta = static_cast<range_impl_scalar_type> (beta);
      const bool conjugate = (mode == Teuchos::CONJ_TRANS);
      const bool transpose = (mode != Teuchos::NO_TRANS);

#ifdef HAVE_TPETRA_DEBUG
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (X.getNumVectors () != Y.getNumVectors (), std::runtime_error,
         "X.getNumVectors() = " << X.getNumVectors () << " != Y.getNumVectors() = "
         << Y.getNumVectors () << ".");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (! transpose && X.getLocalLength () != getColMap ()->getNodeNumElements (),
         std::runtime_error, "NO_TRANS case: X has the wrong number of local rows.  "
         "X.getLocalLength() = " << X.getLocalLength () << " != getColMap()->"
         "getNodeNumElements() = " << getColMap ()->getNodeNumElements () << ".");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (! transpose && Y.getLocalLength () != getRowMap ()->getNodeNumElements (),
         std::runtime_error, "NO_TRANS case: Y has the wrong number of local rows.  "
         "Y.getLocalLength() = " << Y.getLocalLength () << " != getRowMap()->"
         "getNodeNumElements() = " << getRowMap ()->getNodeNumElements () << ".");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (transpose && X.getLocalLength () != getRowMap ()->getNodeNumElements (),
         std::runtime_error, "TRANS or CONJ_TRANS case: X has the wrong number of "
         "local rows.  X.getLocalLength() = " << X.getLocalLength () << " != "
         "getRowMap()->getNodeNumElements() = "
         << getRowMap ()->getNodeNumElements () << ".");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (transpose && Y.getLocalLength () != getColMap ()->getNodeNumElements (),
         std::runtime_error, "TRANS or CONJ_TRANS case: X has the wrong number of "
         "local rows.  Y.getLocalLength() = " << Y.getLocalLength () << " != "
         "getColMap()->getNodeNumElements() = "
         << getColMap ()->getNodeNumElements () << ".");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (! isFillComplete (), std::runtime_error, "The matrix is not fill "
         "complete.  You must call fillComplete() (possibly with domain and range "
         "Map arguments) without an intervening resumeFill() call before you may "
         "call this method.");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (! X.isConstantStride () || ! Y.isConstantStride (), std::runtime_error,
         "X and Y must be constant stride.");
      // If the two pointers are NULL, then they don't alias one
      // another, even though they are equal.
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
        X.getDualView ().d_view.ptr_on_device () == Y.getDualView ().d_view.ptr_on_device () &&
        X.getDualView ().d_view.ptr_on_device () != NULL,
        std::runtime_error, "X and Y may not alias one another.");
#endif // HAVE_TPETRA_DEBUG

      // Y = alpha*op(M) + beta*Y
      if (transpose) {
        KokkosSparse::spmv (conjugate ? KokkosSparse::ConjugateTranspose : KokkosSparse::Transpose,
                            theAlpha,
                            lclMatrix_,
                            X.template getLocalView<device_type> (),
                            theBeta,
                            Y.template getLocalView<device_type> ());
      }
      else {
        KokkosSparse::spmv (KokkosSparse::NoTranspose,
                            theAlpha,
                            lclMatrix_,
                            X.template getLocalView<device_type> (),
                            theBeta,
                            Y.template getLocalView<device_type> ());
      }
    }

    template <class DomainScalar, class RangeScalar>
    void
    localGaussSeidel (const MultiVector<DomainScalar, LocalOrdinal, GlobalOrdinal, Node, classic> &B,
                      MultiVector<RangeScalar, LocalOrdinal, GlobalOrdinal, Node, classic> &X,
                      const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> &D,
                      const RangeScalar& dampingFactor,
                      const KokkosClassic::ESweepDirection direction) const
    {
      typedef LocalOrdinal LO;
      typedef GlobalOrdinal GO;
      typedef Tpetra::MultiVector<DomainScalar, LO, GO, Node, classic> DMV;
      typedef Tpetra::MultiVector<RangeScalar, LO, GO, Node, classic> RMV;
      typedef Tpetra::MultiVector<Scalar, LO, GO, Node, classic> MMV;
      typedef typename DMV::dual_view_type::host_mirror_space HMDT ;
      typedef typename Graph::local_graph_type k_local_graph_type;
      typedef typename k_local_graph_type::size_type offset_type;
      const char prefix[] = "Tpetra::CrsMatrix::localGaussSeidel: ";

      TEUCHOS_TEST_FOR_EXCEPTION
        (! this->isFillComplete (), std::runtime_error,
         prefix << "The matrix is not fill complete.");
      const size_t lclNumRows = this->getNodeNumRows ();
      const size_t numVecs = B.getNumVectors ();
      TEUCHOS_TEST_FOR_EXCEPTION
        (X.getNumVectors () != numVecs, std::invalid_argument,
         prefix << "B.getNumVectors() = " << numVecs << " != "
         "X.getNumVectors() = " << X.getNumVectors () << ".");
      TEUCHOS_TEST_FOR_EXCEPTION
        (B.getLocalLength () != lclNumRows, std::invalid_argument,
         prefix << "B.getLocalLength() = " << B.getLocalLength ()
         << " != this->getNodeNumRows() = " << lclNumRows << ".");

      typename DMV::dual_view_type::t_host B_lcl = B.template getLocalView<HMDT> ();
      typename RMV::dual_view_type::t_host X_lcl = X.template getLocalView<HMDT> ();
      typename MMV::dual_view_type::t_host D_lcl = D.template getLocalView<HMDT> ();

      offset_type B_stride[8], X_stride[8], D_stride[8];
      B_lcl.stride (B_stride);
      X_lcl.stride (X_stride);
      D_lcl.stride (D_stride);

      local_matrix_type lclMatrix = this->getLocalMatrix ();
      k_local_graph_type lclGraph = lclMatrix.graph;
      typename local_matrix_type::row_map_type ptr = lclGraph.row_map;
      typename local_matrix_type::index_type ind = lclGraph.entries;
      typename local_matrix_type::values_type val = lclMatrix.values;
      const offset_type* const ptrRaw = ptr.ptr_on_device ();
      const LO* const indRaw = ind.ptr_on_device ();
      const impl_scalar_type* const valRaw = val.ptr_on_device ();

      const std::string dir ((direction == KokkosClassic::Forward) ? "F" : "B");
      KokkosSparse::Impl::Sequential::gaussSeidel (static_cast<LO> (lclNumRows),
                                                   static_cast<LO> (numVecs),
                                                   ptrRaw, indRaw, valRaw,
                                                   B_lcl.ptr_on_device (), B_stride[1],
                                                   X_lcl.ptr_on_device (), X_stride[1],
                                                   D_lcl.ptr_on_device (),
                                                   static_cast<impl_scalar_type> (dampingFactor),
                                                   dir.c_str ());
    }

    template <class DomainScalar, class RangeScalar>
    void
    reorderedLocalGaussSeidel (const MultiVector<DomainScalar, LocalOrdinal, GlobalOrdinal, Node, classic>& B,
                               MultiVector<RangeScalar, LocalOrdinal, GlobalOrdinal, Node, classic>& X,
                               const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& D,
                               const Teuchos::ArrayView<LocalOrdinal>& rowIndices,
                               const RangeScalar& dampingFactor,
                               const KokkosClassic::ESweepDirection direction) const
    {
      typedef LocalOrdinal LO;
      typedef GlobalOrdinal GO;
      typedef Tpetra::MultiVector<DomainScalar, LO, GO, Node, classic> DMV;
      typedef Tpetra::MultiVector<RangeScalar, LO, GO, Node, classic> RMV;
      typedef Tpetra::MultiVector<Scalar, LO, GO, Node, classic> MMV;
      typedef typename DMV::dual_view_type::host_mirror_space HMDT ;
      typedef typename Graph::local_graph_type k_local_graph_type;
      typedef typename k_local_graph_type::size_type offset_type;
      const char prefix[] = "Tpetra::CrsMatrix::reorderedLocalGaussSeidel: ";

      TEUCHOS_TEST_FOR_EXCEPTION
        (! this->isFillComplete (), std::runtime_error,
         prefix << "The matrix is not fill complete.");
      const size_t lclNumRows = this->getNodeNumRows ();
      const size_t numVecs = B.getNumVectors ();
      TEUCHOS_TEST_FOR_EXCEPTION
        (X.getNumVectors () != numVecs, std::invalid_argument,
         prefix << "B.getNumVectors() = " << numVecs << " != "
         "X.getNumVectors() = " << X.getNumVectors () << ".");
      TEUCHOS_TEST_FOR_EXCEPTION
        (B.getLocalLength () != lclNumRows, std::invalid_argument,
         prefix << "B.getLocalLength() = " << B.getLocalLength ()
         << " != this->getNodeNumRows() = " << lclNumRows << ".");
      TEUCHOS_TEST_FOR_EXCEPTION
        (static_cast<size_t> (rowIndices.size ()) < lclNumRows,
         std::invalid_argument, prefix << "rowIndices.size() = "
         << rowIndices.size () << " < this->getNodeNumRows() = "
         << lclNumRows << ".");

      typename DMV::dual_view_type::t_host B_lcl = B.template getLocalView<HMDT> ();
      typename RMV::dual_view_type::t_host X_lcl = X.template getLocalView<HMDT> ();
      typename MMV::dual_view_type::t_host D_lcl = D.template getLocalView<HMDT> ();

      offset_type B_stride[8], X_stride[8], D_stride[8];
      B_lcl.stride (B_stride);
      X_lcl.stride (X_stride);
      D_lcl.stride (D_stride);

      local_matrix_type lclMatrix = this->getLocalMatrix ();
      typename Graph::local_graph_type lclGraph = lclMatrix.graph;
      typename local_matrix_type::index_type ind = lclGraph.entries;
      typename local_matrix_type::row_map_type ptr = lclGraph.row_map;
      typename local_matrix_type::values_type val = lclMatrix.values;
      const offset_type* const ptrRaw = ptr.ptr_on_device ();
      const LO* const indRaw = ind.ptr_on_device ();
      const impl_scalar_type* const valRaw = val.ptr_on_device ();

      const std::string dir = (direction == KokkosClassic::Forward) ? "F" : "B";
      KokkosSparse::Impl::Sequential::reorderedGaussSeidel (static_cast<LO> (lclNumRows),
                                                            static_cast<LO> (numVecs),
                                                            ptrRaw, indRaw, valRaw,
                                                            B_lcl.ptr_on_device (),
                                                            B_stride[1],
                                                            X_lcl.ptr_on_device (),
                                                            X_stride[1],
                                                            D_lcl.ptr_on_device (),
                                                            rowIndices.getRawPtr (),
                                                            static_cast<LO> (lclNumRows),
                                                            static_cast<impl_scalar_type> (dampingFactor),
                                                            dir.c_str ());
    }

    template <class DomainScalar, class RangeScalar>
    void
    localSolve (const MultiVector<RangeScalar, LocalOrdinal, GlobalOrdinal, Node, classic>& Y,
                MultiVector<DomainScalar, LocalOrdinal, GlobalOrdinal, Node, classic>& X,
                Teuchos::ETransp mode) const
    {
      using Teuchos::CONJ_TRANS;
      using Teuchos::NO_TRANS;
      using Teuchos::TRANS;
      typedef LocalOrdinal LO;
      typedef GlobalOrdinal GO;
      typedef Tpetra::MultiVector<DomainScalar, LO, GO, Node, classic> DMV;
      typedef Tpetra::MultiVector<RangeScalar, LO, GO, Node, classic> RMV;
      typedef typename DMV::dual_view_type::host_mirror_space HMDT ;

      const char tfecfFuncName[] = "localSolve: ";

      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (! isFillComplete (), std::runtime_error,
         "The matrix is not fill complete.");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (! X.isConstantStride () || ! Y.isConstantStride (), std::invalid_argument,
         "X and Y must be constant stride.");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (! isUpperTriangular () && ! isLowerTriangular (), std::runtime_error,
         "The matrix is neither upper triangular or lower triangular.  "
         "You may only call this method if the matrix is triangular.  "
         "Remember that this is a local (per MPI process) property, and that "
         "Tpetra only knows how to do a local (per process) triangular solve.");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
        (STS::isComplex && mode == TRANS, std::logic_error, "This method does "
         "not currently support non-conjugated transposed solve (mode == "
         "Teuchos::TRANS) for complex scalar types.");

      // FIXME (mfh 27 Aug 2014) Tpetra has always made the odd decision
      // that if _some_ diagonal entries are missing locally, then it
      // assumes that the matrix has an implicitly stored unit diagonal.
      // Whether the matrix has an implicit unit diagonal or not should
      // be up to the user to decide.  What if the graph has no diagonal
      // entries, and the user wants it that way?  The only reason this
      // matters, though, is for the triangular solve, and in that case,
      // missing diagonal entries will cause trouble anyway.  However,
      // it would make sense to warn the user if they ask for a
      // triangular solve with an incomplete diagonal.  Furthermore,
      // this code should only assume an implicitly stored unit diagonal
      // if the matrix has _no_ explicitly stored diagonal entries.

      const std::string uplo = isUpperTriangular () ? "U" :
        (isLowerTriangular () ? "L" : "N");
      const std::string trans = (mode == Teuchos::CONJ_TRANS) ? "C" :
        (mode == Teuchos::TRANS ? "T" : "N");
      const std::string diag =
        (getNodeNumDiags () < getNodeNumRows ()) ? "U" : "N";

      local_matrix_type A_lcl = this->getLocalMatrix ();

      // FIXME (mfh 23 Apr 2015) We currently only have a host,
      // sequential kernel for local sparse triangular solve.

      Y.getDualView ().template sync<HMDT> (); // Y is read-only
      X.getDualView ().template sync<HMDT> ();
      X.getDualView ().template modify<HMDT> (); // we will write to X

      if (X.isConstantStride () && Y.isConstantStride ()) {
        typename DMV::dual_view_type::t_host X_lcl = X.template getLocalView<HMDT> ();
        typename RMV::dual_view_type::t_host Y_lcl = Y.template getLocalView<HMDT> ();
        KokkosSparse::trsv (uplo.c_str (), trans.c_str (), diag.c_str (), A_lcl, Y_lcl, X_lcl);
      }
      else {
        const size_t numVecs = std::min (X.getNumVectors (), Y.getNumVectors ());
        for (size_t j = 0; j < numVecs; ++j) {
          auto X_j = X.getVector (j);
          auto Y_j = X.getVector (j);
          auto X_lcl = X_j->template getLocalView<HMDT> ();
          auto Y_lcl = Y_j->template getLocalView<HMDT> ();
          KokkosSparse::trsv (uplo.c_str (), trans.c_str (), diag.c_str (), A_lcl, Y_lcl, X_lcl);
        }
      }
    }

    template <class T>
    Teuchos::RCP<CrsMatrix<T, LocalOrdinal, GlobalOrdinal, Node, classic> >
    convert () const;



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

    bool hasTransposeApply () const;

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

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



    void
    gaussSeidel (const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> &B,
                 MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> &X,
                 const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> &D,
                 const Scalar& dampingFactor,
                 const ESweepDirection direction,
                 const int numSweeps) const;

    void
    reorderedGaussSeidel (const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& B,
                          MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& X,
                          const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& D,
                          const Teuchos::ArrayView<LocalOrdinal>& rowIndices,
                          const Scalar& dampingFactor,
                          const ESweepDirection direction,
                          const int numSweeps) const;

    void
    gaussSeidelCopy (MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> &X,
                     const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> &B,
                     const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> &D,
                     const Scalar& dampingFactor,
                     const ESweepDirection direction,
                     const int numSweeps,
                     const bool zeroInitialGuess) const;

    void
    reorderedGaussSeidelCopy (MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& X,
                              const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& B,
                              const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& D,
                              const Teuchos::ArrayView<LocalOrdinal>& rowIndices,
                              const Scalar& dampingFactor,
                              const ESweepDirection direction,
                              const int numSweeps,
                              const bool zeroInitialGuess) const;

    virtual Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
    add (const Scalar& alpha,
         const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
         const Scalar& beta,
         const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& domainMap,
         const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& rangeMap,
         const Teuchos::RCP<Teuchos::ParameterList>& params) 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<char>& exports,
                    const Teuchos::ArrayView<size_t>& numPacketsPerLID,
                    size_t& constantNumPackets,
                    Distributor& distor);

  private:
    void
    unpackAndCombineImpl (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
                          const Teuchos::ArrayView<const char> &imports,
                          const Teuchos::ArrayView<size_t> &numPacketsPerLID,
                          size_t constantNumPackets,
                          Distributor& distor,
                          CombineMode combineMode);

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


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

  private:
    bool
    packRow (char* const numEntOut,
             char* const valOut,
             char* const indOut,
             const size_t numEnt,
             const LocalOrdinal lclRow) const;

    bool
    unpackRow (Scalar* const valInTmp,
               GlobalOrdinal* const indInTmp,
               const size_t tmpNumEnt,
               const char* const valIn,
               const char* const indIn,
               const size_t numEnt,
               const LocalOrdinal lclRow,
               const Tpetra::CombineMode combineMode);

    void
    allocatePackSpace (Teuchos::Array<char>& exports,
                       size_t& totalNumEntries,
                       const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs) const;

  public:
    typename local_matrix_type::values_type getLocalValuesView () const {
      return k_values1D_;
    }

  private:
    // Friend declaration for nonmember function.
    template<class CrsMatrixType>
    friend Teuchos::RCP<CrsMatrixType>
    importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                                    const Import<typename CrsMatrixType::local_ordinal_type,
                                                 typename CrsMatrixType::global_ordinal_type,
                                                 typename CrsMatrixType::node_type>& importer,
                                    const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                 typename CrsMatrixType::global_ordinal_type,
                                                                 typename CrsMatrixType::node_type> >& domainMap,
                                    const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                 typename CrsMatrixType::global_ordinal_type,
                                                                 typename CrsMatrixType::node_type> >& rangeMap,
                                    const Teuchos::RCP<Teuchos::ParameterList>& params);

    // Friend declaration for nonmember function.
    template<class CrsMatrixType>
    friend Teuchos::RCP<CrsMatrixType>
    exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                                    const Export<typename CrsMatrixType::local_ordinal_type,
                                                 typename CrsMatrixType::global_ordinal_type,
                                                 typename CrsMatrixType::node_type>& exporter,
                                    const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                 typename CrsMatrixType::global_ordinal_type,
                                                                 typename CrsMatrixType::node_type> >& domainMap,
                                    const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                                 typename CrsMatrixType::global_ordinal_type,
                                                                 typename CrsMatrixType::node_type> >& rangeMap,
                                    const Teuchos::RCP<Teuchos::ParameterList>& params);

  public:
    void
    importAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >& destMatrix,
                           const import_type& importer,
                           const Teuchos::RCP<const map_type>& domainMap,
                           const Teuchos::RCP<const map_type>& rangeMap,
                           const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;

    void
    exportAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >& destMatrix,
                           const export_type& exporter,
                           const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
                           const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
                           const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;



  private:
    void
    transferAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >& destMatrix,
                             const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& rowTransfer,
                             const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
                             const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
                             const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
    // We forbid copy construction by declaring this method private
    // and not implementing it.
    CrsMatrix (const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>& rhs);

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

    void
    insertGlobalValuesFiltered (const GlobalOrdinal globalRow,
                                const Teuchos::ArrayView<const GlobalOrdinal>& indices,
                                const Teuchos::ArrayView<const Scalar>& values);

    void
    insertLocalValuesFiltered (const LocalOrdinal localRow,
                               const Teuchos::ArrayView<const LocalOrdinal>& indices,
                               const Teuchos::ArrayView<const Scalar>& values);

    void
    combineGlobalValues (const GlobalOrdinal globalRowIndex,
                         const Teuchos::ArrayView<const GlobalOrdinal>& columnIndices,
                         const Teuchos::ArrayView<const Scalar>& values,
                         const Tpetra::CombineMode combineMode);

    template<class BinaryFunction>
    LocalOrdinal
    transformGlobalValues (GlobalOrdinal globalRow,
                           const Teuchos::ArrayView<const GlobalOrdinal>& indices,
                           const Teuchos::ArrayView<const Scalar>& values,
                           BinaryFunction f)
    {
      using Teuchos::Array;
      using Teuchos::ArrayView;
      typedef LocalOrdinal LO;
      typedef GlobalOrdinal GO;
      typedef impl_scalar_type ST;
      typedef typename ArrayView<const GO>::size_type size_type;

      if (! isFillActive ()) {
        // Fill must be active in order to call this method.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }
      else if (values.size () != indices.size ()) {
        // The sizes of values and indices must match.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }

      const LO lrow = this->getRowMap ()->getLocalElement (globalRow);
      if (lrow == Teuchos::OrdinalTraits<LO>::invalid ()) {
        // We don't own the row, so we're not allowed to modify its values.
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }

      if (staticGraph_.is_null ()) {
        return Teuchos::OrdinalTraits<LO>::invalid ();
      }
      const crs_graph_type& graph = *staticGraph_;
      RowInfo rowInfo = graph.getRowInfo (lrow);
      if (indices.size () == 0) {
        return static_cast<LO> (0);
      }
      else {
        ArrayView<const ST> valsIn =
          Teuchos::av_reinterpret_cast<const ST> (values);
        ArrayView<ST> curVals = this->getViewNonConst (rowInfo);
        if (isLocallyIndexed ()) {
          // Convert the given global indices to local indices.
          //
          // FIXME (mfh 08 Jul 2014) Why can't we ask the graph to do
          // that?  It could do the conversions in place, so that we
          // wouldn't need temporary storage.
          const map_type& colMap = * (this->getColMap ());
          const size_type numInds = indices.size ();
          Array<LO> lclInds (numInds);
          for (size_type k = 0; k < numInds; ++k) {
            // There is no need to filter out indices not in the
            // column Map.  Those that aren't will be mapped to
            // invalid(), which the graph's transformGlobalValues()
            // will filter out (but not count in its return value).
            lclInds[k] = colMap.getLocalElement (indices[k]);
          }
          return graph.template transformLocalValues<ST, BinaryFunction> (rowInfo,
                                                                          curVals,
                                                                          lclInds (),
                                                                          valsIn, f);
        }
        else if (isGloballyIndexed ()) {
          return graph.template transformGlobalValues<ST, BinaryFunction> (rowInfo,
                                                                           curVals,
                                                                           indices,
                                                                           valsIn, f);
        }
        else {
          // If the graph is neither locally nor globally indexed on
          // the calling process, that means that the calling process
          // can't possibly have any entries in the owned row.  Thus,
          // there are no entries to transform, so we return zero.
          return static_cast<LO> (0);
        }
      }
    }

  private:
    void
    insertNonownedGlobalValues (const GlobalOrdinal globalRow,
                                const Teuchos::ArrayView<const GlobalOrdinal>& indices,
                                const Teuchos::ArrayView<const Scalar>& values);

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

  protected:
    // useful typedefs
    typedef Teuchos::OrdinalTraits<LocalOrdinal> OTL;
    typedef Kokkos::Details::ArithTraits<impl_scalar_type> STS;
    typedef Kokkos::Details::ArithTraits<mag_type> STM;
    typedef MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> MV;
    typedef Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>      V;
    typedef crs_graph_type Graph;

    // Enums
    enum GraphAllocationStatus {
      GraphAlreadyAllocated,
      GraphNotYetAllocated
    };

    void allocateValues (ELocalGlobal lg, GraphAllocationStatus gas);

    void sortEntries();

    void mergeRedundantEntries();

    void clearGlobalConstants();

    void computeGlobalConstants();

    mutable Teuchos::RCP<MV> importMV_;

    mutable Teuchos::RCP<MV> exportMV_;

    Teuchos::RCP<MV>
    getColumnMapMultiVector (const MV& X_domainMap,
                             const bool force = false) const;

    Teuchos::RCP<MV>
    getRowMapMultiVector (const MV& Y_rangeMap,
                          const bool force = false) const;

    void
    applyNonTranspose (const MV& X_in,
                       MV& Y_in,
                       Scalar alpha,
                       Scalar beta) const;

    void
    applyTranspose (const MV& X_in,
                    MV& Y_in,
                    const Teuchos::ETransp mode,
                    Scalar alpha,
                    Scalar beta) const;

    // matrix data accessors

    Teuchos::ArrayView<const impl_scalar_type> getView (RowInfo rowinfo) const;

    Teuchos::ArrayView<impl_scalar_type> getViewNonConst (RowInfo rowinfo);

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

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

    void checkInternalState () const;


    Teuchos::RCP<const Graph> staticGraph_;
    Teuchos::RCP<      Graph>     myGraph_;

    local_matrix_type lclMatrix_;


    typename local_matrix_type::values_type k_values1D_;
    Teuchos::ArrayRCP<Teuchos::Array<impl_scalar_type> > values2D_;

    Details::EStorageStatus storageStatus_;

    bool fillComplete_;

    std::map<GlobalOrdinal, std::pair<Teuchos::Array<GlobalOrdinal>,
                                      Teuchos::Array<Scalar> > > nonlocals_;

    mutable mag_type frobNorm_;

  public:
    // FIXME (mfh 24 Feb 2014) Is it _really_ necessary to make this a
    // public inner class of CrsMatrix?  It looks like it doesn't
    // depend on any implementation details of CrsMatrix at all.  It
    // should really be declared and defined outside of CrsMatrix.
    template<class ViewType, class OffsetViewType>
    struct pack_functor {
      typedef typename ViewType::execution_space execution_space;
      ViewType src_;
      ViewType dst_;
      OffsetViewType src_offset_;
      OffsetViewType dst_offset_;
      typedef typename OffsetViewType::non_const_value_type scalar_index_type;

      pack_functor (ViewType dst, ViewType src,
                    OffsetViewType dst_offset, OffsetViewType src_offset) :
        src_ (src),
        dst_ (dst),
        src_offset_ (src_offset),
        dst_offset_ (dst_offset)
      {}

      KOKKOS_INLINE_FUNCTION
      void operator () (const LocalOrdinal row) const {
        scalar_index_type srcPos = src_offset_(row);
        const scalar_index_type dstEnd = dst_offset_(row+1);
        scalar_index_type dstPos = dst_offset_(row);
        for ( ; dstPos < dstEnd; ++dstPos, ++srcPos) {
          dst_(dstPos) = src_(srcPos);
        }
      }
    };
  }; // class CrsMatrix


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

  template<class CrsMatrixType>
  Teuchos::RCP<CrsMatrixType>
  importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                                  const Import<typename CrsMatrixType::local_ordinal_type,
                                               typename CrsMatrixType::global_ordinal_type,
                                               typename CrsMatrixType::node_type>& importer,
                                  const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                               typename CrsMatrixType::global_ordinal_type,
                                                               typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
                                  const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                               typename CrsMatrixType::global_ordinal_type,
                                                               typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
                                  const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
  {
    Teuchos::RCP<CrsMatrixType> destMatrix;
    sourceMatrix->importAndFillComplete (destMatrix,importer, domainMap, rangeMap, params);
    return destMatrix;
  }

  template<class CrsMatrixType>
  Teuchos::RCP<CrsMatrixType>
  exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
                                  const Export<typename CrsMatrixType::local_ordinal_type,
                                               typename CrsMatrixType::global_ordinal_type,
                                               typename CrsMatrixType::node_type>& exporter,
                                  const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                               typename CrsMatrixType::global_ordinal_type,
                                                               typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
                                  const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
                                                               typename CrsMatrixType::global_ordinal_type,
                                                               typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
                                  const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
  {
    Teuchos::RCP<CrsMatrixType> destMatrix;
    sourceMatrix->exportAndFillComplete (destMatrix,exporter, domainMap, rangeMap, params);
    return destMatrix;
  }
} // namespace Tpetra

#endif // TPETRA_CRSMATRIX_DECL_HPP