Ifpack2_Container.hpp

Go to the documentation of this file.

/*@HEADER
// ***********************************************************************
//
//       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2009) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// 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 IFPACK2_CONTAINER_HPP
#define IFPACK2_CONTAINER_HPP


#include "Ifpack2_ConfigDefs.hpp"
#include "Tpetra_RowMatrix.hpp"
#include "Teuchos_Describable.hpp"
#include <Ifpack2_Partitioner.hpp>
#include <Tpetra_Map.hpp>
#include <Tpetra_Experimental_BlockCrsMatrix_decl.hpp>
#include <Teuchos_ParameterList.hpp>
#include <Teuchos_Time.hpp>
#include <iostream>
#include <type_traits>

#ifndef DOXYGEN_SHOULD_SKIP_THIS
namespace Teuchos {
  // Forward declaration to avoid include.
  class ParameterList;
} // namespace Teuchos
#endif // DOXYGEN_SHOULD_SKIP_THIS

namespace Ifpack2 {

template<class MatrixType>
class Container : public Teuchos::Describable {

protected:
  typedef typename MatrixType::scalar_type scalar_type;
  typedef typename MatrixType::local_ordinal_type local_ordinal_type;
  typedef typename MatrixType::global_ordinal_type global_ordinal_type;
  typedef typename MatrixType::node_type node_type;
  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> mv_type;
  typedef Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> vector_type;
  typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
  typedef Teuchos::ScalarTraits<scalar_type> STS;
  typedef Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
  typedef Partitioner<Tpetra::RowGraph<local_ordinal_type, global_ordinal_type, node_type> > partitioner_type;
  typedef Tpetra::Experimental::BlockCrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> block_crs_matrix_type;
  typedef Tpetra::RowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> row_matrix_type;

  static_assert(std::is_same<MatrixType, row_matrix_type>::value,
                "Ifpack2::Container: Please use MatrixType = Tpetra::RowMatrix.");

  typedef typename Kokkos::Details::ArithTraits<scalar_type>::val_type impl_scalar_type;

public:
  typedef typename mv_type::dual_view_type::t_host HostView;

  Container (const Teuchos::RCP<const row_matrix_type>& matrix,
             const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
             const Teuchos::RCP<const import_type>& importer,
             int OverlapLevel,
             scalar_type DampingFactor) :
    inputMatrix_ (matrix),
    OverlapLevel_ (OverlapLevel),
    DampingFactor_ (DampingFactor),
    Importer_ (importer)
  {
    using Teuchos::Ptr;
    using Teuchos::RCP;
    using Teuchos::Array;
    using Teuchos::ArrayView;
    using Teuchos::Comm;
    // typedef Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type> import_type; // unused
    NumLocalRows_ = inputMatrix_->getNodeNumRows();
    NumGlobalRows_ = inputMatrix_->getGlobalNumRows();
    NumGlobalNonzeros_ = inputMatrix_->getGlobalNumEntries();
    IsParallel_ = inputMatrix_->getRowMap()->getComm()->getSize() != 1;
    auto global_bcrs =
      Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(inputMatrix_);
    hasBlockCrs_ = !global_bcrs.is_null();
    if(hasBlockCrs_)
      bcrsBlockSize_ = global_bcrs->getBlockSize();
    else
      bcrsBlockSize_ = 1;
    setBlockSizes(partitions);
  }

  Container (const Teuchos::RCP<const row_matrix_type>& matrix,
             const Teuchos::Array<local_ordinal_type>& localRows) :
    inputMatrix_ (matrix),
    numBlocks_ (1),
    partitions_ (localRows.size()),
    blockRows_ (1),
    partitionIndices_ (1),
    OverlapLevel_ (0),
    DampingFactor_ (STS::one()),
    Importer_ (Teuchos::null)
  {
    NumLocalRows_ = inputMatrix_->getNodeNumRows();
    NumGlobalRows_ = inputMatrix_->getGlobalNumRows();
    NumGlobalNonzeros_ = inputMatrix_->getGlobalNumEntries();
    IsParallel_ = inputMatrix_->getRowMap()->getComm()->getSize() > 1;
    blockRows_[0] = localRows.size();
    partitionIndices_[0] = 0;
    const block_crs_matrix_type* global_bcrs =
      Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(inputMatrix_).get();
    hasBlockCrs_ = global_bcrs;
    if(hasBlockCrs_)
      bcrsBlockSize_ = global_bcrs->getBlockSize();
    else
      bcrsBlockSize_ = 1;
    for(local_ordinal_type i = 0; i < localRows.size(); i++)
      partitions_[i] = localRows[i];
  }

  virtual ~Container() {};

  Teuchos::ArrayView<const local_ordinal_type> getLocalRows(int blockIndex) const
  {
    return Teuchos::ArrayView<const local_ordinal_type>
      (&partitions_[partitionIndices_[blockIndex]], blockRows_[blockIndex]);
  }

  virtual void initialize () = 0;

  void setBlockSizes(const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions)
  {
    //First, create a grand total of all the rows in all the blocks
    //Note: If partitioner allowed overlap, this could be greater than the # of local rows
    local_ordinal_type totalBlockRows = 0;
    numBlocks_ = partitions.size();
    blockRows_.resize(numBlocks_);
    partitionIndices_.resize(numBlocks_);
    for(int i = 0; i < numBlocks_; i++)
    {
      local_ordinal_type rowsInBlock = partitions[i].size();
      blockRows_[i] = rowsInBlock;
      partitionIndices_[i] = totalBlockRows;
      totalBlockRows += rowsInBlock;
    }
    partitions_.resize(totalBlockRows);
    //set partitions_: each entry is the partition/block of the row
    local_ordinal_type iter = 0;
    for(int i = 0; i < numBlocks_; i++)
    {
      for(int j = 0; j < blockRows_[i]; j++)
      {
        partitions_[iter++] = partitions[i][j];
      }
    }
  }

  void getMatDiag() const
  {
    if(Diag_.is_null())
    {
      Diag_ = rcp(new vector_type(inputMatrix_->getDomainMap()));
      inputMatrix_->getLocalDiagCopy(*Diag_);
    }
  }

  virtual void compute () = 0;

  void DoJacobi(HostView& X, HostView& Y, int stride) const;
  void DoOverlappingJacobi(HostView& X, HostView& Y, HostView& W, int stride) const;
  void DoGaussSeidel(HostView& X, HostView& Y, HostView& Y2, int stride) const;
  void DoSGS(HostView& X, HostView& Y, HostView& Y2, int stride) const;

  virtual void setParameters (const Teuchos::ParameterList& List) = 0;

  virtual bool isInitialized () const = 0;

  virtual bool isComputed () const = 0;

  virtual void
  apply (HostView& X,
         HostView& Y,
         int blockIndex,
         int stride,
         Teuchos::ETransp mode = Teuchos::NO_TRANS,
         scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one(),
         scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero()) const = 0;

  // The underlying container implements the splitting <tt>A = D + R</tt>. Only
  // it can have efficient access to D and R, as these are constructed in the
  // symbolic and numeric phases.
  //
  // This is the first performance-portable implementation of a block
  // relaxation, and it is supported currently only by BlockTriDiContainer.
  virtual void applyInverseJacobi (const mv_type& /* X */, mv_type& /* Y */,
                                   bool /* zeroStartingSolution = false */,
                                   int /* numSweeps = 1 */) const
  { TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Not implemented."); }

  void applyMV (mv_type& X, mv_type& Y) const
  {
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE
    HostView XView = X.template getLocalView<Kokkos::HostSpace>();
    HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
#else
    HostView XView = X.getLocalViewHost();
    HostView YView = Y.getLocalViewHost();
#endif
    this->apply (XView, YView, 0, X.getStride());
  }

  virtual void
  weightedApply (HostView& X,
                 HostView& Y,
                 HostView& W,
                 int blockIndex,
                 int stride,
                 Teuchos::ETransp mode = Teuchos::NO_TRANS,
                 scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one(),
                 scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero()) const = 0;

  void weightedApplyMV (mv_type& X,
                        mv_type& Y,
                        vector_type& W)
  {
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE
    HostView XView = X.template getLocalView<Kokkos::HostSpace>();
    HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
    HostView WView = W.template getLocalView<Kokkos::HostSpace>();
#else
    HostView XView = X.getLocalViewHost();
    HostView YView = Y.getLocalViewHost();
    HostView WView = W.getLocalViewHost();
#endif
    weightedApply (XView, YView, WView, 0, X.getStride());
  }

  virtual void clearBlocks();

  virtual std::ostream& print (std::ostream& os) const = 0;

  static std::string getName()
  {
    return "Generic";
  }

protected:
  Teuchos::RCP<const row_matrix_type> inputMatrix_;

  int numBlocks_;
  Teuchos::Array<local_ordinal_type> partitions_;      //size: total # of local rows (in all local blocks)
  Teuchos::Array<local_ordinal_type> blockRows_;     //size: # of blocks
  Teuchos::Array<local_ordinal_type> partitionIndices_;
  mutable Teuchos::RCP<vector_type> Diag_;
  bool IsParallel_;
  int OverlapLevel_;
  scalar_type DampingFactor_;
  Teuchos::RCP<const Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type>> Importer_;
  local_ordinal_type NumLocalRows_;
  global_ordinal_type NumGlobalRows_;
  global_ordinal_type NumGlobalNonzeros_;
  bool hasBlockCrs_;
  int bcrsBlockSize_;
};

template <class MatrixType>
inline std::ostream&
operator<< (std::ostream& os, const Ifpack2::Container<MatrixType>& obj)
{
  return obj.print (os);
}

template <class MatrixType>
void Container<MatrixType>::DoJacobi(HostView& X, HostView& Y, int stride) const
{
  const scalar_type one = STS::one();
  // Note: Flop counts copied naively from Ifpack.
  // use partitions_ and blockRows_
  size_t numVecs = X.extent(1);
  // Non-overlapping Jacobi
  for (local_ordinal_type i = 0; i < numBlocks_; i++)
  {
    // may happen that a partition is empty
    if(blockRows_[i] != 1 || hasBlockCrs_)
    {
      if(blockRows_[i] == 0 )
        continue;
      apply(X, Y, i, stride, Teuchos::NO_TRANS, DampingFactor_, one);
    }
    else    // singleton, can't access Containers_[i] as it was never filled and may be null.
    {
      local_ordinal_type LRID = partitions_[partitionIndices_[i]];
      getMatDiag();
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE
      HostView diagView = Diag_->template getLocalView<Kokkos::HostSpace>();
#else
      HostView diagView = Diag_->getLocalViewHost();
#endif
      impl_scalar_type d = (impl_scalar_type) one / diagView(LRID, 0);
      for(size_t nv = 0; nv < numVecs; nv++)
      {
        impl_scalar_type x = X(LRID, nv);
        Y(LRID, nv) = x * d;
      }
    }
  }
}

template <class MatrixType>
void Container<MatrixType>::DoOverlappingJacobi(HostView& X, HostView& Y, HostView& W, int stride) const
{
  // Overlapping Jacobi
  for(local_ordinal_type i = 0; i < numBlocks_; i++)
  {
    // may happen that a partition is empty
    if(blockRows_[i] == 0)
      continue;
    if(blockRows_[i] != 1)
      weightedApply(X, Y, W, i, stride, Teuchos::NO_TRANS, DampingFactor_, STS::one());
  }
}

template<class MatrixType>
void Container<MatrixType>::DoGaussSeidel(HostView& X, HostView& Y, HostView& Y2, int stride) const
{
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::Ptr;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  // Note: Flop counts copied naively from Ifpack.
  const scalar_type one = STS::one();
  const size_t Length = inputMatrix_->getNodeMaxNumRowEntries();
  auto numVecs = X.extent(1);
  Array<scalar_type> Values;
  Array<local_ordinal_type> Indices;
  Indices.resize(Length);

  Values.resize(bcrsBlockSize_ * bcrsBlockSize_ * Length);    //note: if A was not a BlockCRS, bcrsBlockSize_ is 1
  // I think I decided I need two extra vectors:
  // One to store the sum of the corrections (initialized to zero)
  // One to store the temporary residual (doesn't matter if it is zeroed or not)
  // My apologies for making the names clear and meaningful. (X=RHS, Y=guess?! Nice.)
  HostView Resid("", X.extent(0), X.extent(1));
  for(local_ordinal_type i = 0; i < numBlocks_; i++)
  {
    if(blockRows_[i] > 1 || hasBlockCrs_)
    {
      if (blockRows_[i] == 0)
        continue;
      // update from previous block
      ArrayView<const local_ordinal_type> localRows = getLocalRows (i);
      const size_t localNumRows = blockRows_[i];
      for(size_t j = 0; j < localNumRows; j++)
      {
        const local_ordinal_type LID = localRows[j]; // Containers_[i]->ID (j);
        size_t NumEntries;
        inputMatrix_->getLocalRowCopy(LID, Indices(), Values(), NumEntries);
        for(size_t m = 0; m < numVecs; m++)
        {
          if(hasBlockCrs_)
          {
            for (int localR = 0; localR < bcrsBlockSize_; localR++)
              Resid(LID * bcrsBlockSize_ + localR, m) = X(LID * bcrsBlockSize_ + localR, m);
            for (size_t k = 0; k < NumEntries; ++k)
            {
              const local_ordinal_type col = Indices[k];
              for (int localR = 0; localR < bcrsBlockSize_; localR++)
              {
                for(int localC = 0; localC < bcrsBlockSize_; localC++)
                {
                  Resid(LID * bcrsBlockSize_ + localR, m) -=
                    Values[k * bcrsBlockSize_ * bcrsBlockSize_ + localR + localC * bcrsBlockSize_]
                    * Y2(col * bcrsBlockSize_ + localC, m);
                }
              }
            }
          }
          else
          {
            Resid(LID, m) = X(LID, m);
            for (size_t k = 0; k < NumEntries; ++k)
            {
              const local_ordinal_type col = Indices[k];
              Resid(LID, m) -= Values[k] * Y2(col, m);
            }
          }
        }
      }
      // solve with this block
      //
      // Note: I'm abusing the ordering information, knowing that X/Y
      // and Y2 have the same ordering for on-proc unknowns.
      //
      // Note: Add flop counts for inverse apply
      apply(Resid, Y2, i, stride, Teuchos::NO_TRANS, DampingFactor_, one);
    }
    else if(blockRows_[i] == 1)
    {
      // singleton, can't access Containers_[i] as it was never filled and may be null.
      // a singleton calculation is exact, all residuals should be zero.
      local_ordinal_type LRID = partitionIndices_[i];  // by definition, a singleton 1 row in block.
      getMatDiag();
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE
      HostView diagView = Diag_->template getLocalView<Kokkos::HostSpace>();
#else
      HostView diagView = Diag_->getLocalViewHost();
#endif
      impl_scalar_type d = (impl_scalar_type) one / diagView(LRID, 0);
      for(size_t m = 0; m < numVecs; m++)
      {
        impl_scalar_type x = X(LRID, m);
        impl_scalar_type newy = x * d;
        Y2(LRID, m) = newy;
     }
    } // end else
  } // end for numBlocks_
  if(IsParallel_)
  {
    auto numMyRows = inputMatrix_->getNodeNumRows();
    for (size_t m = 0; m < numVecs; ++m)
    {
      for (size_t i = 0; i < numMyRows * bcrsBlockSize_; ++i)
      {
        Y(i, m) = Y2(i, m);
      }
    }
  }
}

template<class MatrixType>
void Container<MatrixType>::DoSGS(HostView& X, HostView& Y, HostView& Y2, int stride) const
{
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::Ptr;
  using Teuchos::ptr;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  const scalar_type one = STS::one();
  auto numVecs = X.extent(1);
  const size_t Length = inputMatrix_->getNodeMaxNumRowEntries();
  Array<scalar_type> Values;
  Array<local_ordinal_type> Indices(Length);
  Values.resize(bcrsBlockSize_ * bcrsBlockSize_ * Length);
  // I think I decided I need two extra vectors:
  // One to store the sum of the corrections (initialized to zero)
  // One to store the temporary residual (doesn't matter if it is zeroed or not)
  // My apologies for making the names clear and meaningful. (X=RHS, Y=guess?! Nice.)
  HostView Resid("", X.extent(0), X.extent(1));
  // Forward Sweep
  for(local_ordinal_type i = 0; i < numBlocks_; i++)
  {
    if(blockRows_[i] != 1 || hasBlockCrs_)
    {
      if(blockRows_[i] == 0)
        continue; // Skip empty partitions
      // update from previous block
      ArrayView<const local_ordinal_type> localRows = getLocalRows(i);
      for(local_ordinal_type j = 0; j < blockRows_[i]; j++)
      {
        const local_ordinal_type LID = localRows[j]; // Containers_[i]->ID (j);
        size_t NumEntries;
        inputMatrix_->getLocalRowCopy(LID, Indices(), Values(), NumEntries);
        //set tmpresid = initresid - A*correction
        for(size_t m = 0; m < numVecs; m++)
        {
          if(hasBlockCrs_)
          {
            for(int localR = 0; localR < bcrsBlockSize_; localR++)
              Resid(LID * bcrsBlockSize_ + localR, m) = X(LID * bcrsBlockSize_ + localR, m);
            for(size_t k = 0; k < NumEntries; ++k)
            {
              const local_ordinal_type col = Indices[k];
              for (int localR = 0; localR < bcrsBlockSize_; localR++)
              {
                for(int localC = 0; localC < bcrsBlockSize_; localC++)
                  Resid(LID * bcrsBlockSize_ + localR, m) -=
                    Values[k * (bcrsBlockSize_ * bcrsBlockSize_) + (localR + localC * bcrsBlockSize_)]
                    * Y2(col * bcrsBlockSize_ + localC, m);
              }
            }
          }
          else
          {
            Resid(LID, m) = X(LID, m);
            for(size_t k = 0; k < NumEntries; k++)
            {
              local_ordinal_type col = Indices[k];
              Resid(LID, m) -= Values[k] * Y2(col, m);
            }
          }
        }
      }
      // solve with this block
      //
      // Note: I'm abusing the ordering information, knowing that X/Y
      // and Y2 have the same ordering for on-proc unknowns.
      //
      // Note: Add flop counts for inverse apply
      apply(Resid, Y2, i, stride, Teuchos::NO_TRANS, DampingFactor_, one);
      // operations for all getrow's
    }
    else // singleton, can't access Containers_[i] as it was never filled and may be null.
    {
      local_ordinal_type LRID  = partitions_[partitionIndices_[i]];
      getMatDiag();
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE
      HostView diagView = Diag_->template getLocalView<Kokkos::HostSpace>();
#else
      HostView diagView = Diag_->getLocalViewHost();
#endif
      impl_scalar_type d = (impl_scalar_type) one / diagView(LRID, 0);
      for(size_t m = 0; m < numVecs; m++)
      {
        impl_scalar_type x = X(LRID, m);
        impl_scalar_type newy = x * d;
        Y2(LRID, m) = newy;
      }
    } // end else
  } // end forward sweep over NumLocalBlocks
  // Reverse Sweep
  //
  // mfh 12 July 2013: The unusual iteration bounds, and the use of
  // i-1 rather than i in the loop body, ensure correctness even if
  // local_ordinal_type is unsigned.  "i = numBlocks_-1; i >= 0;
  // i--" will loop forever if local_ordinal_type is unsigned, because
  // unsigned integers are (trivially) always nonnegative.
  for(local_ordinal_type i = numBlocks_; i > 0; --i)
  {
    if(hasBlockCrs_ || blockRows_[i-1] != 1)
    {
      if(blockRows_[i - 1] == 0)
        continue;
      // update from previous block
      ArrayView<const local_ordinal_type> localRows = getLocalRows(i-1);
      for(local_ordinal_type j = 0; j < blockRows_[i-1]; j++)
      {
        const local_ordinal_type LID = localRows[j]; // Containers_[i-1]->ID (j);
        size_t NumEntries;
        inputMatrix_->getLocalRowCopy(LID, Indices(), Values(), NumEntries);
        //set tmpresid = initresid - A*correction
        for (size_t m = 0; m < numVecs; m++)
        {
          if(hasBlockCrs_)
          {
            for(int localR = 0; localR < bcrsBlockSize_; localR++)
              Resid(LID * bcrsBlockSize_ + localR, m) = X(LID * bcrsBlockSize_ + localR, m);
            for(size_t k = 0; k < NumEntries; ++k)
            {
              const local_ordinal_type col = Indices[k];
              for(int localR = 0; localR < bcrsBlockSize_; localR++)
              {
                for(int localC = 0; localC < bcrsBlockSize_; localC++)
                  Resid(LID*bcrsBlockSize_+localR, m) -=
                    Values[k * (bcrsBlockSize_ * bcrsBlockSize_) + (localR + localC * bcrsBlockSize_)]
                    * Y2(col * bcrsBlockSize_ + localC, m);
              }
            }
          }
          else
          {
            Resid(LID, m) = X(LID, m);
            for(size_t k = 0; k < NumEntries; ++k)
            {
              local_ordinal_type col = Indices[k];
              Resid(LID, m) -= Values[k] * Y2(col, m);
            }
          }
        }
      }
      // solve with this block
      //
      // Note: I'm abusing the ordering information, knowing that X/Y
      // and Y2 have the same ordering for on-proc unknowns.
      //
      // Note: Add flop counts for inverse apply
      apply(Resid, Y2, i - 1, stride, Teuchos::NO_TRANS, DampingFactor_, one);
      // operations for all getrow's
    } // end  Partitioner_->numRowsInPart (i) != 1 )
    // else do nothing, as by definition with a singleton, the residuals are zero.
  } //end reverse sweep
  if(IsParallel_)
  {
    auto numMyRows = inputMatrix_->getNodeNumRows();
    for (size_t m = 0; m < numVecs; ++m)
    {
      for (size_t i = 0; i < numMyRows * bcrsBlockSize_; ++i)
      {
        Y(i, m) = Y2(i, m);
      }
    }
  }
}

template<class MatrixType>
void Container<MatrixType>::clearBlocks()
{
  numBlocks_ = 0;
  partitions_.clear();
  blockRows_.clear();
  partitionIndices_.clear();
  Diag_ = Teuchos::null;      //Diag_ will be recreated if needed
}

} // namespace Ifpack2

namespace Teuchos {

template<class MatrixType>
class TEUCHOSCORE_LIB_DLL_EXPORT TypeNameTraits< ::Ifpack2::Container<MatrixType> >
{
 public:
  static std::string name () {
    return std::string ("Ifpack2::Container<") +
      TypeNameTraits<MatrixType>::name () + " >";
  }

  static std::string concreteName (const ::Ifpack2::Container<MatrixType>&) {
    return name ();
  }
};

} // namespace Teuchos

#endif // IFPACK2_CONTAINER_HPP