Tpetra_Experimental_BlockCrsMatrix_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_EXPERIMENTAL_BLOCKCRSMATRIX_DECL_HPP
#define TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_DECL_HPP


#include <ctime>
#include <Tpetra_ConfigDefs.hpp>
#include <Tpetra_CrsGraph.hpp>
#include <Tpetra_RowMatrix.hpp>
#include <Tpetra_Experimental_BlockMultiVector.hpp>

namespace Tpetra {
namespace Experimental {

template<class Scalar = Details::DefaultTypes::scalar_type,
         class LO = Details::DefaultTypes::local_ordinal_type,
         class GO = Details::DefaultTypes::global_ordinal_type,
         class Node = Details::DefaultTypes::node_type>
class BlockCrsMatrix :
  virtual public Tpetra::RowMatrix<Scalar, LO, GO, Node>,
  virtual public Tpetra::DistObject<char, LO, GO, Node>
{
private:
  typedef Tpetra::DistObject<char, LO, GO, Node> dist_object_type;
  typedef BlockMultiVector<Scalar, LO, GO, Node> BMV;
  typedef Teuchos::ScalarTraits<Scalar> STS;

protected:
  typedef char packet_type;

public:


  typedef Scalar scalar_type;

  typedef typename BlockMultiVector<Scalar, LO, GO, Node>::impl_scalar_type impl_scalar_type;

  typedef LO local_ordinal_type;
  typedef GO global_ordinal_type;
  typedef Node node_type;

  typedef ::Tpetra::Map<LO, GO, node_type> map_type;
  typedef Tpetra::MultiVector<Scalar, LO, GO, node_type> mv_type;
  typedef Tpetra::CrsGraph<LO, GO, node_type> crs_graph_type;

  typedef LittleBlock<impl_scalar_type, LO> little_block_type;
  typedef LittleBlock<const impl_scalar_type, LO> const_little_block_type;
  typedef LittleVector<impl_scalar_type, LO> little_vec_type;
  typedef LittleVector<const impl_scalar_type, LO> const_little_vec_type;



  BlockCrsMatrix ();

  BlockCrsMatrix (const crs_graph_type& graph, const LO blockSize);

  BlockCrsMatrix (const crs_graph_type& graph,
                  const map_type& domainPointMap,
                  const map_type& rangePointMap,
                  const LO blockSize);

  virtual ~BlockCrsMatrix () {}



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

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

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

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

  global_size_t getGlobalNumRows() const;

  size_t getNodeNumRows() const;

  size_t getNodeMaxNumRowEntries() const;

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

  bool hasTransposeApply () const {
    // FIXME (mfh 04 May 2014) Transpose and conjugate transpose modes
    // are not implemented yet.  Fill in applyBlockTrans() to fix this.
    return false;
  }

  void setAllToScalar (const Scalar& alpha);



  std::string description () const;

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



  LO getBlockSize () const { return blockSize_; }

  virtual Teuchos::RCP<const Tpetra::RowGraph<LO,GO,Node> > getGraph () const;

  const crs_graph_type & getCrsGraph () const { return graph_; }

  void
  applyBlock (const BlockMultiVector<Scalar, LO, GO, Node>& X,
              BlockMultiVector<Scalar, LO, GO, Node>& Y,
              Teuchos::ETransp mode = Teuchos::NO_TRANS,
              const Scalar alpha = Teuchos::ScalarTraits<Scalar>::one (),
              const Scalar beta = Teuchos::ScalarTraits<Scalar>::zero ());

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

  void
  reorderedGaussSeidelCopy (MultiVector<Scalar,LO,GO,Node>& X,
                            const MultiVector<Scalar,LO,GO,Node>& B,
                            const MultiVector<Scalar,LO,GO,Node>& D,
                            const ArrayView<LO>& rowIndices,
                            const Scalar& dampingFactor,
                            const ESweepDirection direction,
                            const int numSweeps,
                            const bool zeroInitialGuess) const;


  void
  localGaussSeidel (const BlockMultiVector<Scalar, LO, GO, Node>& Residual,
                          BlockMultiVector<Scalar, LO, GO, Node>& Solution,
                          BlockCrsMatrix<Scalar, LO, GO, Node> & factorizedDiagonal,
                          const int * factorizationPivots,
                          const Scalar omega,
                          const ESweepDirection direction) const;

  LO
  replaceLocalValues (const LO localRowInd,
                      const LO colInds[],
                      const Scalar vals[],
                      const LO numColInds) const;

  LO
  sumIntoLocalValues (const LO localRowInd,
                      const LO colInds[],
                      const Scalar vals[],
                      const LO numColInds) const;

  LO
  getLocalRowView (const LO localRowInd,
                   const LO*& colInds,
                   Scalar*& vals,
                   LO& numInds) const;

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

  void
  getLocalRowCopy (LO LocalRow,
                   const Teuchos::ArrayView<LO> &Indices,
                   const Teuchos::ArrayView<Scalar> &Values,
                   size_t &NumEntries) const;

  little_block_type
  getLocalBlock (const LO localRowInd, const LO localColInd) const;

  LO
  getLocalRowOffsets (const LO localRowInd,
                      ptrdiff_t offsets[],
                      const LO colInds[],
                      const LO numColInds) const;

  LO
  replaceLocalValuesByOffsets (const LO localRowInd,
                               const ptrdiff_t offsets[],
                               const Scalar vals[],
                               const LO numOffsets) const;

  LO
  sumIntoLocalValuesByOffsets (const LO localRowInd,
                               const ptrdiff_t offsets[],
                               const Scalar vals[],
                               const LO numOffsets) const;

  size_t getNumEntriesInLocalRow (const LO localRowInd) const;

  bool localError () const {
    return *localError_;
  }

  std::string errorMessages () const {
    return (*errs_).is_null () ? std::string ("") : (*errs_)->str ();
  }

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

  void
  getLocalDiagCopy (BlockCrsMatrix<Scalar,LO,GO,Node>& diag,
                    const Teuchos::ArrayView<const size_t>& offsets) const;


  void computeDiagonalGraph ();

  bool isComputedDiagonalGraph() const { return computedDiagonalGraph_;}

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

protected:
  LO
  absMaxLocalValues (const LO localRowInd,
                     const LO colInds[],
                     const Scalar vals[],
                     const LO numColInds) const;

  LO
  absMaxLocalValuesByOffsets (const LO localRowInd,
                              const ptrdiff_t offsets[],
                              const Scalar vals[],
                              const LO numOffsets) const;



  virtual bool checkSizes (const Tpetra::SrcDistObject& source);

  virtual void
  copyAndPermute (const Tpetra::SrcDistObject& source,
                  size_t numSameIDs,
                  const Teuchos::ArrayView<const LO>& permuteToLIDs,
                  const Teuchos::ArrayView<const LO>& permuteFromLIDs);

  virtual void
  packAndPrepare (const Tpetra::SrcDistObject& source,
                  const Teuchos::ArrayView<const LO>& exportLIDs,
                  Teuchos::Array<packet_type>& exports,
                  const Teuchos::ArrayView<size_t>& numPacketsPerLID,
                  size_t& constantNumPackets,
                  Tpetra::Distributor& distor);

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

private:
  crs_graph_type graph_;
  Teuchos::RCP<crs_graph_type> graphRCP_;
  map_type rowMeshMap_;
  map_type domainPointMap_;
  map_type rangePointMap_;
  LO blockSize_;

  typename crs_graph_type::local_graph_type::row_map_type::HostMirror ptr_;

  const LO* ind_;
  Teuchos::ArrayRCP<impl_scalar_type> valView_;
  impl_scalar_type* val_;

  Teuchos::RCP<Teuchos::RCP<BMV> > X_colMap_;
  Teuchos::RCP<Teuchos::RCP<BMV> > Y_rowMap_;

  LO columnPadding_;

  bool rowMajor_;

  Teuchos::RCP<bool> localError_;

  Teuchos::RCP<Teuchos::RCP<std::ostringstream> > errs_;

  std::ostream& markLocalErrorAndGetStream ();

  // //! Clear the local error state and stream.
  // void clearLocalErrorStateAndStream ();

  void
  applyBlockTrans (const BlockMultiVector<Scalar, LO, GO, Node>& X,
                   BlockMultiVector<Scalar, LO, GO, Node>& Y,
                   const Teuchos::ETransp mode,
                   const Scalar alpha,
                   const Scalar beta);

  void
  applyBlockNoTrans (const BlockMultiVector<Scalar, LO, GO, Node>& X,
                     BlockMultiVector<Scalar, LO, GO, Node>& Y,
                     const Scalar alpha,
                     const Scalar beta);

  void
  localApplyBlockNoTrans (const BlockMultiVector<Scalar, LO, GO, Node>& X,
                          BlockMultiVector<Scalar, LO, GO, Node>& Y,
                          const Scalar alpha,
                          const Scalar beta);

  Teuchos::RCP<crs_graph_type> diagonalGraph_;
  bool computedDiagonalGraph_;

  size_t
  findRelOffsetOfColumnIndex (const LO localRowIndex,
                              const LO colIndexToFind,
                              const size_t hint) const;

  LO offsetPerBlock () const;

  const_little_block_type
  getConstLocalBlockFromInput (const impl_scalar_type* val, const size_t pointOffset) const;

  little_block_type
  getNonConstLocalBlockFromInput (impl_scalar_type* val, const size_t pointOffset) const;

  const_little_block_type
  getConstLocalBlockFromAbsOffset (const size_t absBlockOffset) const;

  little_block_type
  getNonConstLocalBlockFromAbsOffset (const size_t absBlockOffset) const;

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

  virtual Teuchos::RCP<Node> getNode() const;

  virtual global_size_t getGlobalNumCols() const;

  virtual size_t getNodeNumCols() const;

  virtual GO getIndexBase() const;

  virtual global_size_t getGlobalNumEntries() const;

  virtual size_t getNodeNumEntries() const;

  virtual size_t getNumEntriesInGlobalRow (GO globalRow) const;

  virtual global_size_t getGlobalNumDiags() const;

  virtual size_t getNodeNumDiags() const;

  virtual size_t getGlobalMaxNumRowEntries() const;

  virtual bool hasColMap() const;

  virtual bool isLowerTriangular() const;

  virtual bool isUpperTriangular() const;

  virtual bool isLocallyIndexed() const;

  virtual bool isGloballyIndexed() const;

  virtual bool isFillComplete() const;

  virtual bool supportsRowViews() const;



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

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

  virtual void getLocalDiagCopy (Vector<Scalar,LO,GO,Node>& diag) const;



  virtual void leftScale (const Vector<Scalar, LO, GO, Node>& x);

  virtual void rightScale (const Vector<Scalar, LO, GO, Node>& x);

  virtual typename Tpetra::RowMatrix<Scalar, LO, GO, Node>::mag_type
  getFrobeniusNorm () const;
};


  template<class Scalar, class LO, class GO, class Node>
  void blockCrsMatrixWriter(BlockCrsMatrix<Scalar,LO,GO,Node> const &A, std::string const &fileName) {
    Teuchos::ParameterList pl;
    std::ofstream out;
    out.open(fileName.c_str());
    blockCrsMatrixWriter(A, out, pl);
  }

  template<class Scalar, class LO, class GO, class Node>
  void blockCrsMatrixWriter(BlockCrsMatrix<Scalar,LO,GO,Node> const &A, std::string const &fileName, Teuchos::ParameterList const &params) {
    std::ofstream out;
    out.open(fileName.c_str());
    blockCrsMatrixWriter(A, out, params);
  }

  template<class Scalar, class LO, class GO, class Node>
  void blockCrsMatrixWriter(BlockCrsMatrix<Scalar,LO,GO,Node> const &A, std::ostream &os) {
    Teuchos::ParameterList pl;
    blockCrsMatrixWriter(A, os, pl);
  }

  template<class Scalar, class LO, class GO, class Node>
  void blockCrsMatrixWriter(BlockCrsMatrix<Scalar,LO,GO,Node> const &A, std::ostream &os, Teuchos::ParameterList const &params) {

    using Teuchos::RCP;
    using Teuchos::rcp;

    typedef Teuchos::OrdinalTraits<GO>                     TOT;
    typedef BlockCrsMatrix<Scalar, LO, GO, Node>           block_crs_matrix_type;
    typedef Tpetra::Import<LO, GO, Node>                   import_type;
    typedef Tpetra::Map<LO, GO, Node>                      map_type;
    typedef Tpetra::MultiVector<GO, LO, GO, Node>          mv_type;
    typedef Tpetra::CrsGraph<LO, GO, Node>                 crs_graph_type;

    RCP<const map_type> const rowMap = A.getRowMap(); //"mesh" map
    RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
    const int myRank = comm->getRank();
    const size_t numProcs = comm->getSize();

    // If true, force use of the import strip-mining infrastructure.  This is useful for debugging on one process.
    bool alwaysUseParallelAlgorithm = false;
    if (params.isParameter("always use parallel algorithm"))
      alwaysUseParallelAlgorithm = params.get<bool>("always use parallel algorithm");
    bool printMatrixMarketHeader = true;
    if (params.isParameter("print MatrixMarket header"))
      printMatrixMarketHeader = params.get<bool>("print MatrixMarket header");

    if (printMatrixMarketHeader && myRank==0) {
      std::time_t now = std::time(NULL);
      os << "%%MatrixMarket matrix coordinate real general" << std::endl;
      os << "% time stamp: " << ctime(&now);
      os << "% written from " << numProcs << " processes" << std::endl;
      os << "% point representation of Tpetra::Experimental::BlockCrsMatrix" << std::endl;
      size_t numRows = A.getGlobalNumRows();
      size_t numCols = A.getGlobalNumCols();
      os << "% " << numRows << " block rows, " << numCols << " block columns" << std::endl;
      const LO blockSize = A.getBlockSize();
      os << "% block size " << blockSize << std::endl;
      os << numRows*blockSize << " " << numCols*blockSize << " " << A.getGlobalNumEntries()*blockSize*blockSize << std::endl;
    }

    if (numProcs==1 && !alwaysUseParallelAlgorithm) {
      writeMatrixStrip(A,os,params);
    } else {
      size_t numRows = rowMap->getNodeNumElements();

      //Create source map
      RCP<const map_type> allMeshGidsMap = rcp(new map_type(TOT::invalid(), numRows, A.getIndexBase(), comm));
      //Create and populate vector of mesh GIDs corresponding to this pid's rows.
      //This vector will be imported one pid's worth of information at a time to pid 0.
      mv_type allMeshGids(allMeshGidsMap,1);
      Teuchos::ArrayRCP<GO> allMeshGidsData = allMeshGids.getDataNonConst(0);

      for (size_t i=0; i<numRows; i++)
        allMeshGidsData[i] = rowMap->getGlobalElement(i);
      allMeshGidsData = Teuchos::null;

      // Now construct a RowMatrix on PE 0 by strip-mining the rows of the input matrix A.
      size_t stripSize = allMeshGids.getGlobalLength() / numProcs;
      size_t remainder = allMeshGids.getGlobalLength() % numProcs;
      size_t curStart = 0;
      size_t curStripSize = 0;
      Teuchos::Array<GO> importMeshGidList;
      for (size_t i=0; i<numProcs; i++) {
        if (myRank==0) { // Only PE 0 does this part
          curStripSize = stripSize;
          if (i<remainder) curStripSize++; // handle leftovers
          importMeshGidList.resize(curStripSize); // Set size of vector to max needed
          for (size_t j=0; j<curStripSize; j++) importMeshGidList[j] = j + curStart + A.getIndexBase();
          curStart += curStripSize;
        }
        // The following import map should be non-trivial only on PE 0.
        TEUCHOS_TEST_FOR_EXCEPTION(myRank>0 && curStripSize!=0,
          std::runtime_error, "Tpetra::Experimental::blockCrsMatrixWriter: (pid "
          << myRank << ") map size should be zero, but is " << curStripSize);
        RCP<map_type> importMeshGidMap = rcp(new map_type(TOT::invalid(), importMeshGidList(), A.getIndexBase(), comm));
        import_type gidImporter(allMeshGidsMap, importMeshGidMap);
        mv_type importMeshGids(importMeshGidMap, 1);
        importMeshGids.doImport(allMeshGids, gidImporter, INSERT);

        // importMeshGids now has a list of GIDs for the current strip of matrix rows.
        // Use these values to build another importer that will get rows of the matrix.

        // The following import map will be non-trivial only on PE 0.
        Teuchos::ArrayRCP<const GO> importMeshGidsData = importMeshGids.getData(0);
        Teuchos::Array<GO> importMeshGidsGO;
        importMeshGidsGO.reserve(importMeshGidsData.size());
        for (typename Teuchos::ArrayRCP<const GO>::size_type j=0; j<importMeshGidsData.size(); ++j)
          importMeshGidsGO.push_back(importMeshGidsData[j]);
        RCP<const map_type> importMap = rcp(new map_type(TOT::invalid(), importMeshGidsGO(), rowMap->getIndexBase(), comm) );

        import_type importer(rowMap,importMap );
        size_t numEntriesPerRow = A.getCrsGraph().getGlobalMaxNumRowEntries();
        RCP<crs_graph_type> graph = createCrsGraph(importMap,numEntriesPerRow);
        RCP<const map_type> domainMap = A.getCrsGraph().getDomainMap();
        graph->doImport(A.getCrsGraph(), importer, INSERT);
        graph->fillComplete(domainMap, importMap);

        block_crs_matrix_type importA(*graph, A.getBlockSize());
        importA.doImport(A, importer, INSERT);

        // Finally we are ready to write this strip of the matrix
        writeMatrixStrip(importA, os, params);
      }
    }
  }

  template<class Scalar, class LO, class GO, class Node>
  void writeMatrixStrip(BlockCrsMatrix<Scalar,LO,GO,Node> const &A, std::ostream &os, Teuchos::ParameterList const &params) {

    typedef Tpetra::Map<LO, GO, Node>                      map_type;

    size_t numRows = A.getGlobalNumRows();
    RCP<const map_type> rowMap = A.getRowMap();
    RCP<const map_type> colMap = A.getColMap();
    RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
    const int myRank = comm->getRank();

    const size_t meshRowOffset = rowMap->getIndexBase();
    const size_t meshColOffset = colMap->getIndexBase();
    TEUCHOS_TEST_FOR_EXCEPTION(meshRowOffset != meshColOffset,
      std::runtime_error, "Tpetra::Experimental::writeMatrixStrip: "
      "mesh row index base != mesh column index base");

    if (myRank !=0) {

      TEUCHOS_TEST_FOR_EXCEPTION(A.getNodeNumRows() != 0,
        std::runtime_error, "Tpetra::Experimental::writeMatrixStrip: pid "
        << myRank << " should have 0 rows but has " << A.getNodeNumRows());
      TEUCHOS_TEST_FOR_EXCEPTION(A.getNodeNumCols() != 0,
        std::runtime_error, "Tpetra::Experimental::writeMatrixStrip: pid "
        << myRank << " should have 0 columns but has " << A.getNodeNumCols());

    } else {

      TEUCHOS_TEST_FOR_EXCEPTION(numRows != A.getNodeNumRows(),
        std::runtime_error, "Tpetra::Experimental::writeMatrixStrip: "
        "number of rows on pid 0 does not match global number of rows");


      int err = 0;
      const LO blockSize = A.getBlockSize();
      const size_t numLocalRows = A.getNodeNumRows();
      bool precisionChanged=false;
      int oldPrecision;
      if (params.isParameter("precision")) {
        oldPrecision = os.precision(params.get<int>("precision"));
        precisionChanged=true;
      }
      int pointOffset = 1;
      if (params.isParameter("zero-based indexing")) {
        if (params.get<bool>("zero-based indexing") == true)
          pointOffset = 0;
      }

      size_t localRowInd;
      for (localRowInd = 0; localRowInd < numLocalRows; ++localRowInd) {

        // Get a view of the current row.
        const LO*     localColInds;
        Scalar* vals;
        LO numEntries;
        err = A.getLocalRowView (localRowInd, localColInds, vals, numEntries);
        if (err != 0)
          break;
        GO globalMeshRowID = rowMap->getGlobalElement(localRowInd) - meshRowOffset;

        for (LO k = 0; k < numEntries; ++k) {
          GO globalMeshColID = colMap->getGlobalElement(localColInds[k]) - meshColOffset;
          Scalar* const curBlock = vals + blockSize * blockSize * k;
          // Blocks are stored in row-major format.
          for (LO j = 0; j < blockSize; ++j) {
            GO globalPointRowID = globalMeshRowID * blockSize + j + pointOffset;
            for (LO i = 0; i < blockSize; ++i) {
              GO globalPointColID = globalMeshColID * blockSize + i + pointOffset;
              const Scalar curVal = curBlock[i + j * blockSize];
              os << globalPointRowID << " " << globalPointColID << " " << curVal << std::endl;
            }
          }
        }
      }
      if (precisionChanged)
        os.precision(oldPrecision);
      TEUCHOS_TEST_FOR_EXCEPTION(err != 0,
        std::runtime_error, "Tpetra::Experimental::writeMatrixStrip: "
        "error getting view of local row " << localRowInd);

    }

  }

} // namespace Experimental
} // namespace Tpetra

#endif // TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_DECL_HPP