Tpetra_BlockCrsMatrix_Helpers_def.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef TPETRA_BLOCKCRSMATRIX_HELPERS_DEF_HPP
#define TPETRA_BLOCKCRSMATRIX_HELPERS_DEF_HPP


#include "Tpetra_BlockCrsMatrix.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_HashTable.hpp"
#include "Tpetra_Import.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_ScalarTraits.hpp"
#include <ctime>
#include <fstream>

namespace Tpetra {

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

    const std::string dataTypeStr =
        Teuchos::ScalarTraits<Scalar>::isComplex ? "complex" : "real";

    // Explanation of first line of file:
    // - "%%MatrixMarket" is the tag for Matrix Market format.
    // - "matrix" is what we're printing.
    // - "coordinate" means sparse (triplet format), rather than dense.
    // - "real" / "complex" is the type (in an output format sense,
    //   not in a C++ sense) of each value of the matrix.  (The
    //   other two possibilities are "integer" (not really necessary
    //   here) and "pattern" (no values, just graph).)
    os << "%%MatrixMarket matrix coordinate " << dataTypeStr << " general" << std::endl;
    os << "% time stamp: " << ctime(&now);
    os << "% written from " << numProcs << " processes" << std::endl;
    os << "% point representation of Tpetra::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->getLocalNumElements();

    // 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::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) {
  using Teuchos::RCP;
  using map_type                       = Tpetra::Map<LO, GO, Node>;
  using bcrs_type                      = BlockCrsMatrix<Scalar, LO, GO, Node>;
  using bcrs_local_inds_host_view_type = typename bcrs_type::local_inds_host_view_type;
  using bcrs_values_host_view_type     = typename bcrs_type::values_host_view_type;
  using impl_scalar_type               = typename bcrs_type::impl_scalar_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::writeMatrixStrip: "
                             "mesh row index base != mesh column index base");

  if (myRank != 0) {
    TEUCHOS_TEST_FOR_EXCEPTION(A.getLocalNumRows() != 0,
                               std::runtime_error, "Tpetra::writeMatrixStrip: pid " << myRank << " should have 0 rows but has " << A.getLocalNumRows());
    TEUCHOS_TEST_FOR_EXCEPTION(A.getLocalNumCols() != 0,
                               std::runtime_error, "Tpetra::writeMatrixStrip: pid " << myRank << " should have 0 columns but has " << A.getLocalNumCols());

  } else {
    TEUCHOS_TEST_FOR_EXCEPTION(numRows != A.getLocalNumRows(),
                               std::runtime_error,
                               "Tpetra::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.getLocalNumRows();
    bool precisionChanged     = false;
    int oldPrecision          = 0;  // avoid "unused variable" warning
    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.
      bcrs_local_inds_host_view_type localColInds;
      bcrs_values_host_view_type vals;
      LO numEntries;
      A.getLocalRowView(localRowInd, localColInds, vals);
      numEntries         = localColInds.extent(0);
      GO globalMeshRowID = rowMap->getGlobalElement(localRowInd) - meshRowOffset;

      for (LO k = 0; k < numEntries; ++k) {
        GO globalMeshColID               = colMap->getGlobalElement(localColInds[k]) - meshColOffset;
        const impl_scalar_type *curBlock = vals.data() + 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 impl_scalar_type curVal = curBlock[i + j * blockSize];

            os << globalPointRowID << " " << globalPointColID << " ";
            if (Teuchos::ScalarTraits<impl_scalar_type>::isComplex) {
              // Matrix Market format wants complex values to be
              // written as space-delimited pairs.  See Bug 6469.
              os << Teuchos::ScalarTraits<impl_scalar_type>::real(curVal) << " "
                 << Teuchos::ScalarTraits<impl_scalar_type>::imag(curVal);
            } else {
              os << curVal;
            }
            os << std::endl;
          }
        }
      }
    }
    if (precisionChanged)
      os.precision(oldPrecision);
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0,
                               std::runtime_error,
                               "Tpetra::writeMatrixStrip: "
                               "error getting view of local row "
                                   << localRowInd);
  }
}

template <class Scalar, class LO, class GO, class Node>
Teuchos::RCP<Tpetra::CrsGraph<LO, GO, Node>>
getBlockCrsGraph(const Tpetra::CrsMatrix<Scalar, LO, GO, Node> &pointMatrix, const LO &blockSize, bool use_LID) {
  TEUCHOS_FUNC_TIME_MONITOR_DIFF("Tpetra::getBlockCrsGraph", getBlockCrsGraph0);
  /*
    ASSUMPTIONS:

       1) In point matrix, all entries associated with a little block are present (even if they are zero).
       2) For given mesh DOF, point DOFs appear consecutively and in ascending order in row & column maps.
       3) Point column map and block column map are ordered consistently.
  */

  using Teuchos::Array;
  using Teuchos::ArrayView;
  using Teuchos::RCP;

  typedef Tpetra::Map<LO, GO, Node> map_type;
  typedef Tpetra::CrsGraph<LO, GO, Node> crs_graph_type;
  typedef Tpetra::CrsMatrix<Scalar, LO, GO, Node> crs_matrix_type;

  using local_graph_device_type = typename crs_matrix_type::local_graph_device_type;
  using row_map_type            = typename local_graph_device_type::row_map_type::non_const_type;
  using entries_type            = typename local_graph_device_type::entries_type::non_const_type;

  using offset_type = typename row_map_type::non_const_value_type;

  using execution_space = typename Node::execution_space;
  using range_type      = Kokkos::RangePolicy<execution_space, LO>;

  const map_type &pointRowMap = *(pointMatrix.getRowMap());
  RCP<const map_type> meshRowMap, meshColMap, meshDomainMap, meshRangeMap;

  const map_type &pointColMap    = *(pointMatrix.getColMap());
  const map_type &pointDomainMap = *(pointMatrix.getDomainMap());
  const map_type &pointRangeMap  = *(pointMatrix.getRangeMap());

  {
    TEUCHOS_FUNC_TIME_MONITOR_DIFF("Tpetra::getBlockCrsGraph::createMeshMaps", getBlockCrsGraph1);
    meshRowMap    = createMeshMap<LO, GO, Node>(blockSize, pointRowMap, use_LID);
    meshColMap    = createMeshMap<LO, GO, Node>(blockSize, pointColMap, use_LID);
    meshDomainMap = createMeshMap<LO, GO, Node>(blockSize, pointDomainMap, use_LID);
    meshRangeMap  = createMeshMap<LO, GO, Node>(blockSize, pointRangeMap, use_LID);
    Kokkos::DefaultExecutionSpace().fence();
  }

  if (meshColMap.is_null()) throw std::runtime_error("ERROR: Cannot create mesh colmap");

  auto localMeshColMap  = meshColMap->getLocalMap();
  auto localPointColMap = pointColMap.getLocalMap();
  auto localPointRowMap = pointRowMap.getLocalMap();

  RCP<crs_graph_type> meshCrsGraph;

  const offset_type bs2 = blockSize * blockSize;

  if (use_LID) {
    TEUCHOS_FUNC_TIME_MONITOR_DIFF("Tpetra::getBlockCrsGraph::LID", getBlockCrsGraph2);
    auto pointLocalGraph = pointMatrix.getCrsGraph()->getLocalGraphDevice();
    auto pointRowptr     = pointLocalGraph.row_map;
    auto pointColind     = pointLocalGraph.entries;

    TEUCHOS_TEST_FOR_EXCEPTION(pointColind.extent(0) % bs2 != 0,
                               std::runtime_error,
                               "Tpetra::getBlockCrsGraph: "
                               "local number of non zero entries is not a multiple of blockSize^2 ");

    LO block_rows = pointRowptr.extent(0) == 0 ? 0 : (pointRowptr.extent(0) - 1) / blockSize;
    row_map_type blockRowptr("blockRowptr", block_rows + 1);
    entries_type blockColind("blockColind", pointColind.extent(0) / (bs2));

    Kokkos::parallel_for(
        "fillMesh", range_type(0, block_rows), KOKKOS_LAMBDA(const LO i) {
          const LO offset_b     = pointRowptr(i * blockSize) / bs2;
          const LO offset_b_max = pointRowptr((i + 1) * blockSize) / bs2;

          if (i == block_rows - 1)
            blockRowptr(i + 1) = offset_b_max;
          blockRowptr(i) = offset_b;

          const LO offset_p = pointRowptr(i * blockSize);

          for (LO k = 0; k < offset_b_max - offset_b; ++k) {
            blockColind(offset_b + k) = pointColind(offset_p + k * blockSize) / blockSize;
          }
        });

    meshCrsGraph = rcp(new crs_graph_type(meshRowMap, meshColMap, blockRowptr, blockColind));
    meshCrsGraph->fillComplete(meshDomainMap, meshRangeMap);
    Kokkos::DefaultExecutionSpace().fence();
  } else {
    TEUCHOS_FUNC_TIME_MONITOR_DIFF("Tpetra::getBlockCrsGraph::GID", getBlockCrsGraph3);
    auto pointLocalGraph = pointMatrix.getCrsGraph()->getLocalGraphDevice();
    auto pointRowptr     = pointLocalGraph.row_map;
    auto pointColind     = pointLocalGraph.entries;

    TEUCHOS_TEST_FOR_EXCEPTION(pointColind.extent(0) % bs2 != 0,
                               std::runtime_error,
                               "Tpetra::getBlockCrsGraph: "
                               "local number of non zero entries is not a multiple of blockSize^2 ");

    LO block_rows = pointRowptr.extent(0) == 0 ? 0 : (pointRowptr.extent(0) - 1) / blockSize;
    row_map_type blockRowptr("blockRowptr", block_rows + 1);
    entries_type blockColind("blockColind", pointColind.extent(0) / (bs2));

    Kokkos::parallel_for(
        "fillMesh", range_type(0, block_rows), KOKKOS_LAMBDA(const LO i) {
          const LO offset_b     = pointRowptr(i * blockSize) / bs2;
          const LO offset_b_max = pointRowptr((i + 1) * blockSize) / bs2;

          if (i == block_rows - 1)
            blockRowptr(i + 1) = offset_b_max;
          blockRowptr(i) = offset_b;

          const LO offset_p     = pointRowptr(i * blockSize);
          const LO offset_p_max = pointRowptr((i + 1) * blockSize);

          LO filled_block = 0;
          for (LO p_i = 0; p_i < offset_p_max - offset_p; ++p_i) {
            auto bcol_GID = localPointColMap.getGlobalElement(pointColind(offset_p + p_i)) / blockSize;
            auto bcol_LID = localMeshColMap.getLocalElement(bcol_GID);

            bool visited = false;
            for (LO k = 0; k < filled_block; ++k) {
              if (blockColind(offset_b + k) == bcol_LID)
                visited = true;
            }
            if (!visited) {
              blockColind(offset_b + filled_block) = bcol_LID;
              ++filled_block;
            }
          }
        });

    meshCrsGraph = rcp(new crs_graph_type(meshRowMap, meshColMap, blockRowptr, blockColind));
    meshCrsGraph->fillComplete(meshDomainMap, meshRangeMap);
    Kokkos::DefaultExecutionSpace().fence();
  }

  return meshCrsGraph;
}

template <class Scalar, class LO, class GO, class Node>
Teuchos::RCP<BlockCrsMatrix<Scalar, LO, GO, Node>>
convertToBlockCrsMatrix(const Tpetra::CrsMatrix<Scalar, LO, GO, Node> &pointMatrix, const LO &blockSize, bool use_LID) {
  TEUCHOS_FUNC_TIME_MONITOR_DIFF("Tpetra::convertToBlockCrsMatrix", convertToBlockCrsMatrix0);
  /*
    ASSUMPTIONS:

       1) In point matrix, all entries associated with a little block are present (even if they are zero).
       2) For given mesh DOF, point DOFs appear consecutively and in ascending order in row & column maps.
       3) Point column map and block column map are ordered consistently.
  */

  using Teuchos::Array;
  using Teuchos::ArrayView;
  using Teuchos::RCP;

  typedef Tpetra::BlockCrsMatrix<Scalar, LO, GO, Node> block_crs_matrix_type;
  typedef Tpetra::CrsMatrix<Scalar, LO, GO, Node> crs_matrix_type;

  using local_graph_device_type  = typename crs_matrix_type::local_graph_device_type;
  using local_matrix_device_type = typename crs_matrix_type::local_matrix_device_type;
  using row_map_type             = typename local_graph_device_type::row_map_type::non_const_type;
  using values_type              = typename local_matrix_device_type::values_type::non_const_type;

  using offset_type = typename row_map_type::non_const_value_type;

  using execution_space = typename Node::execution_space;
  using range_type      = Kokkos::RangePolicy<execution_space, LO>;

  RCP<block_crs_matrix_type> blockMatrix;

  const offset_type bs2 = blockSize * blockSize;

  auto meshCrsGraph = getBlockCrsGraph(pointMatrix, blockSize, use_LID);

  if (use_LID) {
    TEUCHOS_FUNC_TIME_MONITOR_DIFF("Tpetra::convertToBlockCrsMatrix::LID", convertToBlockCrsMatrix1);
    auto pointLocalGraph = pointMatrix.getCrsGraph()->getLocalGraphDevice();
    auto pointRowptr     = pointLocalGraph.row_map;
    auto pointColind     = pointLocalGraph.entries;

    offset_type block_rows = pointRowptr.extent(0) == 0 ? 0 : (pointRowptr.extent(0) - 1) / blockSize;
    values_type blockValues("values", meshCrsGraph->getLocalNumEntries() * bs2);
    auto pointValues = pointMatrix.getLocalValuesDevice(Access::ReadOnly);
    auto blockRowptr = meshCrsGraph->getLocalGraphDevice().row_map;

    Kokkos::parallel_for(
        "copyblockValues", range_type(0, block_rows), KOKKOS_LAMBDA(const LO i) {
          const offset_type blkBeg    = blockRowptr[i];
          const offset_type numBlocks = blockRowptr[i + 1] - blkBeg;

          // For each block in the row...
          for (offset_type block = 0; block < numBlocks; block++) {
            // For each entry in the block...
            for (LO little_row = 0; little_row < blockSize; little_row++) {
              offset_type point_row_offset = pointRowptr[i * blockSize + little_row];
              for (LO little_col = 0; little_col < blockSize; little_col++) {
                blockValues((blkBeg + block) * bs2 + little_row * blockSize + little_col) =
                    pointValues[point_row_offset + block * blockSize + little_col];
              }
            }
          }
        });
    blockMatrix = rcp(new block_crs_matrix_type(*meshCrsGraph, blockValues, blockSize));
    Kokkos::DefaultExecutionSpace().fence();
  } else {
    TEUCHOS_FUNC_TIME_MONITOR_DIFF("Tpetra::convertToBlockCrsMatrix::GID", convertToBlockCrsMatrix2);
    auto localMeshColMap  = meshCrsGraph->getColMap()->getLocalMap();
    auto localPointColMap = pointMatrix.getColMap()->getLocalMap();

    values_type blockValues("values", meshCrsGraph->getLocalNumEntries() * bs2);
    {
      auto pointLocalGraph = pointMatrix.getCrsGraph()->getLocalGraphDevice();
      auto pointRowptr     = pointLocalGraph.row_map;
      auto pointColind     = pointLocalGraph.entries;

      offset_type block_rows = pointRowptr.extent(0) == 0 ? 0 : (pointRowptr.extent(0) - 1) / blockSize;
      auto pointValues       = pointMatrix.getLocalValuesDevice(Access::ReadOnly);
      auto blockRowptr       = meshCrsGraph->getLocalGraphDevice().row_map;
      auto blockColind       = meshCrsGraph->getLocalGraphDevice().entries;

      row_map_type pointGColind("pointGColind", pointColind.extent(0));

      Kokkos::parallel_for(
          "computePointGColind", range_type(0, pointColind.extent(0)), KOKKOS_LAMBDA(const LO i) {
            pointGColind(i) = localPointColMap.getGlobalElement(pointColind(i));
          });

      row_map_type blockGColind("blockGColind", blockColind.extent(0));

      Kokkos::parallel_for(
          "computeBlockGColind", range_type(0, blockGColind.extent(0)), KOKKOS_LAMBDA(const LO i) {
            blockGColind(i) = localMeshColMap.getGlobalElement(blockColind(i));
          });

      Kokkos::parallel_for(
          "copyblockValues", range_type(0, block_rows), KOKKOS_LAMBDA(const LO i) {
            const offset_type blkBeg    = blockRowptr[i];
            const offset_type numBlocks = blockRowptr[i + 1] - blkBeg;

            for (offset_type point_i = 0; point_i < pointRowptr[i * blockSize + 1] - pointRowptr[i * blockSize]; point_i++) {
              offset_type block_inv      = static_cast<offset_type>(-1);
              offset_type little_col_inv = static_cast<offset_type>(-1);
              for (offset_type block_2 = 0; block_2 < numBlocks; block_2++) {
                for (LO little_col_2 = 0; little_col_2 < blockSize; little_col_2++) {
                  if (blockGColind(blkBeg + block_2) * blockSize + little_col_2 == pointGColind(pointRowptr[i * blockSize] + point_i)) {
                    block_inv      = block_2;
                    little_col_inv = little_col_2;
                    break;
                  }
                }
                if (block_inv != static_cast<offset_type>(-1))
                  break;
              }

              for (LO little_row = 0; little_row < blockSize; little_row++) {
                blockValues((blkBeg + block_inv) * bs2 + little_row * blockSize + little_col_inv) = pointValues[pointRowptr[i * blockSize + little_row] + point_i];
              }
            }
          });
    }
    blockMatrix = rcp(new block_crs_matrix_type(*meshCrsGraph, blockValues, blockSize));
    Kokkos::DefaultExecutionSpace().fence();
  }

  return blockMatrix;
}

template <class Scalar, class LO, class GO, class Node>
Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LO, GO, Node>>
fillLogicalBlocks(const Tpetra::CrsMatrix<Scalar, LO, GO, Node> &pointMatrix, const LO &blockSize) {
  using crs_t           = Tpetra::CrsMatrix<Scalar, LO, GO, Node>;
  using execution_space = typename Node::execution_space;
  using dev_row_view_t  = typename crs_t::local_graph_device_type::row_map_type::non_const_type;
  using dev_col_view_t  = typename crs_t::local_graph_device_type::entries_type::non_const_type;
  using dev_val_view_t  = typename crs_t::local_matrix_device_type::values_type::non_const_type;
  using range_type      = Kokkos::RangePolicy<execution_space, size_t>;
  using team_policy     = Kokkos::TeamPolicy<execution_space>;
  using member_type     = typename team_policy::member_type;
  using scratch_view    = Kokkos::View<bool *, typename execution_space::scratch_memory_space>;
  using Ordinal         = typename dev_row_view_t::non_const_value_type;

  // Get structure / values
  auto local       = pointMatrix.getLocalMatrixDevice();
  auto row_ptrs    = local.graph.row_map;
  auto col_inds    = local.graph.entries;
  auto values      = local.values;
  const auto nrows = pointMatrix.getLocalNumRows();

  //
  // Populate all active blocks, they must be fully populated with entries so fill with zeroes
  //

  // Make row workspace views
  dev_row_view_t new_rowmap("new_rowmap", nrows + 1);
  const auto blocks_per_row = nrows / blockSize;  // assumes square matrix
  dev_row_view_t active_block_row_map("active_block_row_map", blocks_per_row + 1);
  const int max_threads = execution_space().concurrency();
  assert(blockSize > 1);
  assert(nrows % blockSize == 0);
  const int mem_level = 1;
  const int bytes     = scratch_view::shmem_size(blocks_per_row);

  if (max_threads >= blockSize) {
    // Prefer 1 team per block since this will require a lot less scratch memory
    team_policy tp(blocks_per_row, blockSize);

    // Count active blocks
    Kokkos::parallel_for(
        "countActiveBlocks", tp.set_scratch_size(mem_level, Kokkos::PerTeam(bytes)), KOKKOS_LAMBDA(const member_type &team) {
          Ordinal block_row = team.league_rank();

          scratch_view row_block_active(team.team_scratch(mem_level), blocks_per_row);
          Kokkos::single(
              Kokkos::PerTeam(team), [&]() {
                for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
                  row_block_active(row_block_idx) = false;
                }
              });
          team.team_barrier();

          // All threads in a team scan a blocks-worth of rows to see which
          // blocks are active
          Kokkos::parallel_for(
              Kokkos::TeamThreadRange(team, blockSize), [&](Ordinal block_offset) {
                Ordinal row     = block_row * blockSize + block_offset;
                Ordinal row_itr = row_ptrs(row);
                Ordinal row_end = row_ptrs(row + 1);

                for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
                  const Ordinal first_possible_col_in_block      = row_block_idx * blockSize;
                  const Ordinal first_possible_col_in_next_block = (row_block_idx + 1) * blockSize;
                  Ordinal curr_nnz_col                           = col_inds(row_itr);
                  while (curr_nnz_col >= first_possible_col_in_block && curr_nnz_col < first_possible_col_in_next_block && row_itr < row_end) {
                    // This block has at least one nnz in this row
                    row_block_active(row_block_idx) = true;
                    ++row_itr;
                    if (row_itr == row_end) break;
                    curr_nnz_col = col_inds(row_itr);
                  }
                }
              });

          team.team_barrier();

          Kokkos::single(
              Kokkos::PerTeam(team), [&]() {
                Ordinal count = 0;
                for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
                  if (row_block_active(row_block_idx)) {
                    ++count;
                  }
                }
                active_block_row_map(block_row) = count;
              });
        });
  } else {
    // We don't have enough parallelism to make a thread team handle a block, so just
    // have 1 thread per block
    team_policy tp(blocks_per_row, 1);

    // Count active blocks
    Kokkos::parallel_for(
        "countActiveBlocks", tp.set_scratch_size(mem_level, Kokkos::PerTeam(bytes)), KOKKOS_LAMBDA(const member_type &team) {
          Ordinal block_row = team.league_rank();

          scratch_view row_block_active(team.team_scratch(mem_level), blocks_per_row);
          for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
            row_block_active(row_block_idx) = false;
          }

          // One thread scans a blocks-worth of rows to see which blocks are active
          for (int block_offset = 0; block_offset < blockSize; ++block_offset) {
            Ordinal row     = block_row * blockSize + block_offset;
            Ordinal row_itr = row_ptrs(row);
            Ordinal row_end = row_ptrs(row + 1);

            for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
              const Ordinal first_possible_col_in_block      = row_block_idx * blockSize;
              const Ordinal first_possible_col_in_next_block = (row_block_idx + 1) * blockSize;
              Ordinal curr_nnz_col                           = col_inds(row_itr);
              while (curr_nnz_col >= first_possible_col_in_block && curr_nnz_col < first_possible_col_in_next_block && row_itr < row_end) {
                // This block has at least one nnz in this row
                row_block_active(row_block_idx) = true;
                ++row_itr;
                if (row_itr == row_end) break;
                curr_nnz_col = col_inds(row_itr);
              }
            }
          }

          Ordinal count = 0;
          for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
            if (row_block_active(row_block_idx)) {
              ++count;
            }
          }
          active_block_row_map(block_row) = count;
        });
  }

  Ordinal nnz_block_count = 0;
  KokkosKernels::Impl::kk_exclusive_parallel_prefix_sum<
      execution_space>(active_block_row_map.extent(0), active_block_row_map, nnz_block_count);
  dev_col_view_t block_col_ids("block_col_ids", nnz_block_count);

  // Find active blocks
  if (max_threads >= blockSize) {
    // Prefer 1 team per block since this will require a lot less scratch memory
    team_policy tp(blocks_per_row, blockSize);

    Kokkos::parallel_for(
        "findActiveBlocks", tp.set_scratch_size(mem_level, Kokkos::PerTeam(bytes)), KOKKOS_LAMBDA(const member_type &team) {
          Ordinal block_row = team.league_rank();

          scratch_view row_block_active(team.team_scratch(mem_level), blocks_per_row);
          Kokkos::single(
              Kokkos::PerTeam(team), [&]() {
                for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
                  row_block_active(row_block_idx) = false;
                }
              });
          team.team_barrier();

          // All threads in a team scan a blocks-worth of rows to see which
          // blocks are active
          Kokkos::parallel_for(
              Kokkos::TeamThreadRange(team, blockSize), [&](Ordinal block_offset) {
                Ordinal row     = block_row * blockSize + block_offset;
                Ordinal row_itr = row_ptrs(row);
                Ordinal row_end = row_ptrs(row + 1);

                for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
                  const Ordinal first_possible_col_in_block      = row_block_idx * blockSize;
                  const Ordinal first_possible_col_in_next_block = (row_block_idx + 1) * blockSize;
                  Ordinal curr_nnz_col                           = col_inds(row_itr);
                  while (curr_nnz_col >= first_possible_col_in_block && curr_nnz_col < first_possible_col_in_next_block && row_itr < row_end) {
                    // This block has at least one nnz in this row
                    row_block_active(row_block_idx) = true;
                    ++row_itr;
                    if (row_itr == row_end) break;
                    curr_nnz_col = col_inds(row_itr);
                  }
                }
              });

          team.team_barrier();

          Kokkos::single(
              Kokkos::PerTeam(team), [&]() {
                Ordinal offset = active_block_row_map[block_row];
                for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
                  if (row_block_active(row_block_idx)) {
                    block_col_ids(offset) = row_block_idx;
                    ++offset;
                  }
                }
              });
        });
  } else {
    team_policy tp(blocks_per_row, 1);

    Kokkos::parallel_for(
        "findActiveBlocks", tp.set_scratch_size(mem_level, Kokkos::PerTeam(bytes)), KOKKOS_LAMBDA(const member_type &team) {
          Ordinal block_row = team.league_rank();

          scratch_view row_block_active(team.team_scratch(mem_level), blocks_per_row);
          for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
            row_block_active(row_block_idx) = false;
          }

          // One thread scans a blocks-worth of rows to see which blocks are active
          for (int block_offset = 0; block_offset < blockSize; ++block_offset) {
            Ordinal row     = block_row * blockSize + block_offset;
            Ordinal row_itr = row_ptrs(row);
            Ordinal row_end = row_ptrs(row + 1);

            for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
              const Ordinal first_possible_col_in_block      = row_block_idx * blockSize;
              const Ordinal first_possible_col_in_next_block = (row_block_idx + 1) * blockSize;
              Ordinal curr_nnz_col                           = col_inds(row_itr);
              while (curr_nnz_col >= first_possible_col_in_block && curr_nnz_col < first_possible_col_in_next_block && row_itr < row_end) {
                // This block has at least one nnz in this row
                row_block_active(row_block_idx) = true;
                ++row_itr;
                if (row_itr == row_end) break;
                curr_nnz_col = col_inds(row_itr);
              }
            }
          }

          Ordinal offset = active_block_row_map[block_row];
          for (size_t row_block_idx = 0; row_block_idx < blocks_per_row; ++row_block_idx) {
            if (row_block_active(row_block_idx)) {
              block_col_ids(offset) = row_block_idx;
              ++offset;
            }
          }
        });
  }

  // Sizing
  Kokkos::parallel_for(
      "sizing", range_type(0, nrows), KOKKOS_LAMBDA(const size_t row) {
        const auto block_row          = row / blockSize;
        const Ordinal block_row_begin = active_block_row_map(block_row);
        const Ordinal block_row_end   = active_block_row_map(block_row + 1);
        const Ordinal row_nnz         = (block_row_end - block_row_begin) * blockSize;
        new_rowmap(row)               = row_nnz;
      });

  Ordinal new_nnz_count = 0;
  KokkosKernels::Impl::kk_exclusive_parallel_prefix_sum<
      execution_space>(new_rowmap.extent(0), new_rowmap, new_nnz_count);
  // Now populate cols and vals
  dev_col_view_t new_col_ids("new_col_ids", new_nnz_count);
  dev_val_view_t new_vals("new_vals", new_nnz_count);
  Kokkos::parallel_for(
      "entries", range_type(0, nrows), KOKKOS_LAMBDA(const size_t row) {
        Ordinal row_itr     = row_ptrs(row);
        Ordinal row_end     = row_ptrs(row + 1);
        Ordinal row_itr_new = new_rowmap(row);

        Ordinal block_row       = row / blockSize;
        Ordinal block_row_begin = active_block_row_map(block_row);
        Ordinal block_row_end   = active_block_row_map(block_row + 1);

        for (Ordinal row_block_idx = block_row_begin; row_block_idx < block_row_end; ++row_block_idx) {
          const Ordinal block_col                        = block_col_ids(row_block_idx);
          const Ordinal first_possible_col_in_block      = block_col * blockSize;
          const Ordinal first_possible_col_in_next_block = (block_col + 1) * blockSize;
          for (Ordinal possible_col = first_possible_col_in_block; possible_col < first_possible_col_in_next_block; ++possible_col, ++row_itr_new) {
            new_col_ids(row_itr_new) = possible_col;
            Ordinal curr_nnz_col     = col_inds(row_itr);
            if (possible_col == curr_nnz_col && row_itr < row_end) {
              // Already a non-zero entry
              new_vals(row_itr_new) = values(row_itr);
              ++row_itr;
            }
          }
        }
      });

  // Create new, filled CRS
  auto crs_row_map           = pointMatrix.getRowMap();
  auto crs_col_map           = pointMatrix.getColMap();
  Teuchos::RCP<crs_t> result = Teuchos::rcp(new crs_t(crs_row_map, crs_col_map, new_rowmap, new_col_ids, new_vals));
  result->fillComplete();
  return result;
}

template <class Scalar, class LO, class GO, class Node>
Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LO, GO, Node>>
unfillFormerBlockCrs(const Tpetra::CrsMatrix<Scalar, LO, GO, Node> &pointMatrix) {
  using crs_t          = Tpetra::CrsMatrix<Scalar, LO, GO, Node>;
  using dev_row_view_t = typename crs_t::local_graph_device_type::row_map_type::non_const_type;
  using dev_col_view_t = typename crs_t::local_graph_device_type::entries_type::non_const_type;
  using dev_val_view_t = typename crs_t::local_matrix_device_type::values_type::non_const_type;
#if KOKKOS_VERSION >= 40799
  using impl_scalar_t = typename KokkosKernels::ArithTraits<Scalar>::val_type;
#else
  using impl_scalar_t = typename Kokkos::ArithTraits<Scalar>::val_type;
#endif
#if KOKKOS_VERSION >= 40799
  using STS = KokkosKernels::ArithTraits<impl_scalar_t>;
#else
  using STS           = Kokkos::ArithTraits<impl_scalar_t>;
#endif
  using Ordinal         = typename dev_row_view_t::non_const_value_type;
  using execution_space = typename Node::execution_space;
  using range_type      = Kokkos::RangePolicy<execution_space, size_t>;

  // Get structure / values
  auto local       = pointMatrix.getLocalMatrixHost();
  auto row_ptrs    = local.graph.row_map;
  auto col_inds    = local.graph.entries;
  auto values      = local.values;
  const auto nrows = pointMatrix.getLocalNumRows();

  // Sizing and rows
  dev_row_view_t new_rowmap("new_rowmap", nrows + 1);
  const impl_scalar_t zero = STS::zero();
  Kokkos::parallel_for(
      "sizing", range_type(0, nrows), KOKKOS_LAMBDA(const size_t row) {
        const Ordinal row_nnz_begin = row_ptrs(row);
        Ordinal row_nnzs            = 0;
        for (Ordinal row_nnz = row_nnz_begin; row_nnz < row_ptrs(row + 1); ++row_nnz) {
          const impl_scalar_t value = values(row_nnz);
          if (value != zero) {
            ++row_nnzs;
          }
        }
        new_rowmap(row) = row_nnzs;
      });

  Ordinal real_nnzs = 0;
  KokkosKernels::Impl::kk_exclusive_parallel_prefix_sum<
      execution_space>(new_rowmap.extent(0), new_rowmap, real_nnzs);
  // Now populate cols and vals
  dev_col_view_t new_col_ids("new_col_ids", real_nnzs);
  dev_val_view_t new_vals("new_vals", real_nnzs);
  Kokkos::parallel_for(
      "entries", range_type(0, nrows), KOKKOS_LAMBDA(const size_t row) {
        Ordinal row_nnzs                = 0;
        const Ordinal old_row_nnz_begin = row_ptrs(row);
        const Ordinal new_row_nnz_begin = new_rowmap(row);
        for (Ordinal old_row_nnz = old_row_nnz_begin; old_row_nnz < row_ptrs(row + 1); ++old_row_nnz) {
          const impl_scalar_t value = values(old_row_nnz);
          if (value != zero) {
            new_col_ids(new_row_nnz_begin + row_nnzs) = col_inds(old_row_nnz);
            new_vals(new_row_nnz_begin + row_nnzs)    = value;
            ++row_nnzs;
          }
        }
      });

  // Create new, unfilled CRS
  auto crs_row_map           = pointMatrix.getRowMap();
  auto crs_col_map           = pointMatrix.getColMap();
  Teuchos::RCP<crs_t> result = Teuchos::rcp(new crs_t(crs_row_map, crs_col_map, new_rowmap, new_col_ids, new_vals));
  result->fillComplete();
  return result;
}

template <class LO, class GO, class Node>
Teuchos::RCP<const Tpetra::Map<LO, GO, Node>>
createMeshMap(const LO &blockSize, const Tpetra::Map<LO, GO, Node> &pointMap, bool use_LID) {
  typedef Teuchos::OrdinalTraits<Tpetra::global_size_t> TOT;
  typedef Tpetra::Map<LO, GO, Node> map_type;

  if (use_LID) {
    using execution_space = typename Node::execution_space;
    using range_type      = Kokkos::RangePolicy<execution_space, size_t>;

    auto pointGlobalID = pointMap.getMyGlobalIndicesDevice();
    LO block_rows      = pointGlobalID.extent(0) / blockSize;
    Kokkos::View<GO *, typename map_type::device_type> meshGlobalID("meshGlobalID", block_rows);

    Kokkos::parallel_for(
        "fillMeshMap", range_type(0, block_rows), KOKKOS_LAMBDA(const LO i) {
          meshGlobalID(i) = pointGlobalID(i * blockSize) / blockSize;
        });

    Teuchos::RCP<const map_type> meshMap = Teuchos::rcp(new map_type(TOT::invalid(), meshGlobalID, 0, pointMap.getComm()));
    return meshMap;
  } else {
    // calculate mesh GIDs
    Teuchos::ArrayView<const GO> pointGids = pointMap.getLocalElementList();
    Teuchos::Array<GO> meshGids;
    GO indexBase = pointMap.getIndexBase();

    // Use hash table to track whether we've encountered this GID previously.  This will happen
    // when striding through the point DOFs in a block.  It should not happen otherwise.
    // I don't use sort/make unique because I don't want to change the ordering.
    meshGids.reserve(pointGids.size());
    Tpetra::Details::HashTable<GO, int> hashTable(pointGids.size());
    for (int i = 0; i < pointGids.size(); ++i) {
      GO meshGid = (pointGids[i] - indexBase) / blockSize + indexBase;
      if (hashTable.get(meshGid) == -1) {
        hashTable.add(meshGid, 1);  //(key,value)
        meshGids.push_back(meshGid);
      }
    }

    Teuchos::RCP<const map_type> meshMap = Teuchos::rcp(new map_type(TOT::invalid(), meshGids(), 0, pointMap.getComm()));
    return meshMap;
  }
}

template <class LO, class GO, class Node>
Teuchos::RCP<const Tpetra::Map<LO, GO, Node>>
createPointMap(const LO &blockSize, const Tpetra::Map<LO, GO, Node> &blockMap) {
  typedef Teuchos::OrdinalTraits<Tpetra::global_size_t> TOT;
  typedef Tpetra::Map<LO, GO, Node> map_type;

  // calculate mesh GIDs
  Teuchos::ArrayView<const GO> blockGids = blockMap.getLocalElementList();
  Teuchos::Array<GO> pointGids(blockGids.size() * blockSize);
  GO indexBase = blockMap.getIndexBase();

  for (LO i = 0, ct = 0; i < (LO)blockGids.size(); i++) {
    GO base = (blockGids[i] - indexBase) * blockSize + indexBase;
    for (LO j = 0; j < blockSize; j++, ct++)
      pointGids[i * blockSize + j] = base + j;
  }

  Teuchos::RCP<const map_type> pointMap = Teuchos::rcp(new map_type(TOT::invalid(), pointGids(), 0, blockMap.getComm()));
  return pointMap;
}

template <class Scalar, class LO, class GO, class Node>
Teuchos::RCP<CrsMatrix<Scalar, LO, GO, Node>>
convertToCrsMatrix(const Tpetra::BlockCrsMatrix<Scalar, LO, GO, Node> &blockMatrix) {
  using Teuchos::Array;
  using Teuchos::ArrayView;
  using Teuchos::RCP;

  typedef Tpetra::BlockCrsMatrix<Scalar, LO, GO, Node> block_crs_matrix_type;
  typedef Tpetra::Map<LO, GO, Node> map_type;
  typedef Tpetra::CrsMatrix<Scalar, LO, GO, Node> crs_matrix_type;

  using crs_graph_type           = typename block_crs_matrix_type::crs_graph_type;
  using local_graph_device_type  = typename crs_matrix_type::local_graph_device_type;
  using local_matrix_device_type = typename crs_matrix_type::local_matrix_device_type;
  using row_map_type             = typename local_graph_device_type::row_map_type::non_const_type;
  using entries_type             = typename local_graph_device_type::entries_type::non_const_type;
  using values_type              = typename local_matrix_device_type::values_type::non_const_type;

  using row_map_type_const = typename local_graph_device_type::row_map_type;
  using entries_type_const = typename local_graph_device_type::entries_type;

  using little_block_type = typename block_crs_matrix_type::const_little_block_type;
  using offset_type       = typename row_map_type::non_const_value_type;
  using column_type       = typename entries_type::non_const_value_type;

  using execution_space = typename Node::execution_space;
  using range_type      = Kokkos::RangePolicy<execution_space, LO>;

  LO blocksize          = blockMatrix.getBlockSize();
  const offset_type bs2 = blocksize * blocksize;
  size_t block_nnz      = blockMatrix.getLocalNumEntries();
  size_t point_nnz      = block_nnz * bs2;

  // We can get these from the blockMatrix directly
  RCP<const map_type> pointDomainMap = blockMatrix.getDomainMap();
  RCP<const map_type> pointRangeMap  = blockMatrix.getRangeMap();

  // We need to generate the row/col Map ourselves.
  RCP<const map_type> blockRowMap = blockMatrix.getRowMap();
  RCP<const map_type> pointRowMap = createPointMap<LO, GO, Node>(blocksize, *blockRowMap);

  RCP<const map_type> blockColMap = blockMatrix.getColMap();
  RCP<const map_type> pointColMap = createPointMap<LO, GO, Node>(blocksize, *blockColMap);

  // Get the last few things

  const crs_graph_type &blockGraph = blockMatrix.getCrsGraph();
  LO point_rows                    = (LO)pointRowMap->getLocalNumElements();
  LO block_rows                    = (LO)blockRowMap->getLocalNumElements();
  auto blockValues                 = blockMatrix.getValuesDevice();
  auto blockLocalGraph             = blockGraph.getLocalGraphDevice();
  row_map_type_const blockRowptr   = blockLocalGraph.row_map;
  entries_type_const blockColind   = blockLocalGraph.entries;

  // Generate the point matrix rowptr / colind / values
  row_map_type rowptr("row_map", point_rows + 1);
  entries_type colind("entries", point_nnz);
  values_type values("values", point_nnz);

  // Pre-fill the rowptr
  Kokkos::parallel_for(
      "fillRowPtr", range_type(0, block_rows), KOKKOS_LAMBDA(const LO i) {
        offset_type base = blockRowptr[i];
        offset_type diff = blockRowptr[i + 1] - base;
        for (LO j = 0; j < blocksize; j++) {
          rowptr[i * blocksize + j] = base * bs2 + j * diff * blocksize;
        }

        // Fill the last guy, if we're on the final entry
        if (i == block_rows - 1) {
          rowptr[block_rows * blocksize] = blockRowptr[i + 1] * bs2;
        }
      });

  Kokkos::parallel_for(
      "copyEntriesAndValues", range_type(0, block_rows), KOKKOS_LAMBDA(const LO i) {
        const offset_type blkBeg    = blockRowptr[i];
        const offset_type numBlocks = blockRowptr[i + 1] - blkBeg;

        // For each block in the row...
        for (offset_type block = 0; block < numBlocks; block++) {
          column_type point_col_base = blockColind[blkBeg + block] * blocksize;
          little_block_type my_block(blockValues.data() + (blkBeg + block) * bs2, blocksize, blocksize);

          // For each entry in the block...
          for (LO little_row = 0; little_row < blocksize; little_row++) {
            offset_type point_row_offset = rowptr[i * blocksize + little_row];
            for (LO little_col = 0; little_col < blocksize; little_col++) {
              colind[point_row_offset + block * blocksize + little_col] = point_col_base + little_col;
              values[point_row_offset + block * blocksize + little_col] = my_block(little_row, little_col);
            }
          }
        }
      });

  // Generate the matrix
  RCP<crs_matrix_type> pointCrsMatrix = rcp(new crs_matrix_type(pointRowMap, pointColMap, 0));
  pointCrsMatrix->setAllValues(rowptr, colind, values);

  // FUTURE OPTIMIZATION: Directly compute import/export, rather than letting ESFC do it
  pointCrsMatrix->expertStaticFillComplete(pointDomainMap, pointRangeMap);

  return pointCrsMatrix;
}

}  // namespace Tpetra

//
// Explicit instantiation macro for blockCrsMatrixWriter (various
// overloads), writeMatrixStrip, and convertToBlockCrsMatrix.
//
// Must be expanded from within the Tpetra namespace!
//
#define TPETRA_BLOCKCRSMATRIX_HELPERS_INSTANT(S, LO, GO, NODE)                                                                                                      \
  template void blockCrsMatrixWriter(BlockCrsMatrix<S, LO, GO, NODE> const &A, std::string const &fileName);                                                        \
  template void blockCrsMatrixWriter(BlockCrsMatrix<S, LO, GO, NODE> const &A, std::string const &fileName, Teuchos::ParameterList const &params);                  \
  template void blockCrsMatrixWriter(BlockCrsMatrix<S, LO, GO, NODE> const &A, std::ostream &os);                                                                   \
  template void blockCrsMatrixWriter(BlockCrsMatrix<S, LO, GO, NODE> const &A, std::ostream &os, Teuchos::ParameterList const &params);                             \
  template void writeMatrixStrip(BlockCrsMatrix<S, LO, GO, NODE> const &A, std::ostream &os, Teuchos::ParameterList const &params);                                 \
  template Teuchos::RCP<BlockCrsMatrix<S, LO, GO, NODE>> convertToBlockCrsMatrix(const CrsMatrix<S, LO, GO, NODE> &pointMatrix, const LO &blockSize, bool use_LID); \
  template Teuchos::RCP<CrsMatrix<S, LO, GO, NODE>> fillLogicalBlocks(const CrsMatrix<S, LO, GO, NODE> &pointMatrix, const LO &blockSize);                          \
  template Teuchos::RCP<CrsMatrix<S, LO, GO, NODE>> unfillFormerBlockCrs(const CrsMatrix<S, LO, GO, NODE> &pointMatrix);                                            \
  template Teuchos::RCP<CrsMatrix<S, LO, GO, NODE>> convertToCrsMatrix(const BlockCrsMatrix<S, LO, GO, NODE> &blockMatrix);                                         \
  template Teuchos::RCP<CrsGraph<LO, GO, NODE>> getBlockCrsGraph(const Tpetra::CrsMatrix<S, LO, GO, NODE> &pointMatrix, const LO &blockSize, bool use_local_ID);

//
// Explicit instantiation macro for createMeshMap / createPointMap
//
#define TPETRA_CREATEMESHMAP_INSTANT(LO, GO, NODE)                                                                                         \
  template Teuchos::RCP<const Map<LO, GO, NODE>> createMeshMap(const LO &blockSize, const Map<LO, GO, NODE> &pointMap, bool use_local_ID); \
  template Teuchos::RCP<const Map<LO, GO, NODE>> createPointMap(const LO &blockSize, const Map<LO, GO, NODE> &blockMap);

#endif  // TPETRA_BLOCKCRSMATRIX_HELPERS_DEF_HPP