MatrixMarket_Tpetra_def.hpp

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef MATRIXMARKET_TPETRA_DEF_HPP
#define MATRIXMARKET_TPETRA_DEF_HPP
 
#include "MatrixMarket_Tpetra_decl.hpp"
#include "Teuchos_OrdinalTraits.hpp"
#include "Tpetra_Details_gathervPrint.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_Operator.hpp"
#include "Tpetra_Vector.hpp"
#include "Tpetra_ComputeGatherMap.hpp"
#include "Teuchos_MatrixMarket_Raw_Adder.hpp"
#include "Teuchos_MatrixMarket_Raw_Graph_Adder.hpp"
#include "Teuchos_MatrixMarket_SymmetrizingAdder.hpp"
#include "Teuchos_MatrixMarket_SymmetrizingGraphAdder.hpp"
#include "Teuchos_MatrixMarket_assignScalar.hpp"
#include "Teuchos_MatrixMarket_Banner.hpp"
#include "Teuchos_MatrixMarket_CoordDataReader.hpp"
#include "Teuchos_SetScientific.hpp"
#include "Teuchos_TimeMonitor.hpp"
 
extern "C" {
#include "mmio_Tpetra.h"
}
#include "Tpetra_Distribution.hpp"
 
#include <algorithm>
#include <fstream>
#include <iostream>
#include <iterator>
#include <vector>
#include <stdexcept>
#include <numeric>
 
namespace Tpetra {
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::makeRangeMap(const trcp_tcomm_t& pComm,
                                                                            const global_ordinal_type numRows) {
  // Return a conventional, uniformly partitioned, contiguous map.
  return rcp(new map_type(static_cast<global_size_t>(numRows),
                          static_cast<global_ordinal_type>(0),
                          pComm, GloballyDistributed));
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::makeRowMap(const Teuchos::RCP<const map_type>& pRowMap,
                                                                          const trcp_tcomm_t& pComm,
                                                                          const global_ordinal_type numRows) {
  // If the caller didn't provide a map, return a conventional,
  // uniformly partitioned, contiguous map.
  if (pRowMap.is_null()) {
    return rcp(new map_type(static_cast<global_size_t>(numRows),
                            static_cast<global_ordinal_type>(0),
                            pComm, GloballyDistributed));
  } else {
    TEUCHOS_TEST_FOR_EXCEPTION(!pRowMap->isDistributed() && pComm->getSize() > 1,
                               std::invalid_argument,
                               "The specified row map is not distributed, "
                               "but the given communicator includes more than one process (in "
                               "fact, there are "
                                   << pComm->getSize() << " processes).");
    TEUCHOS_TEST_FOR_EXCEPTION(pRowMap->getComm() != pComm, std::invalid_argument,
                               "The specified row Map's communicator (pRowMap->getComm()) "
                               "differs from the given separately supplied communicator pComm.");
    return pRowMap;
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::makeDomainMap(const Teuchos::RCP<const map_type>& pRangeMap,
                                                                             const global_ordinal_type numRows,
                                                                             const global_ordinal_type numCols) {
  // Abbreviations so that the map creation call isn't too long.
  typedef local_ordinal_type LO;
  typedef global_ordinal_type GO;
  typedef node_type NT;
 
  if (numRows == numCols) {
    return pRangeMap;
  } else {
    return createUniformContigMapWithNode<LO, GO, NT>(numCols,
                                                      pRangeMap->getComm());
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::distribute(Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
                                                                               Teuchos::ArrayRCP<size_t>& myRowPtr,
                                                                               Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
                                                                               Teuchos::ArrayRCP<scalar_type>& myValues,
                                                                               const Teuchos::RCP<const map_type>& pRowMap,
                                                                               Teuchos::ArrayRCP<size_t>& numEntriesPerRow,
                                                                               Teuchos::ArrayRCP<size_t>& rowPtr,
                                                                               Teuchos::ArrayRCP<global_ordinal_type>& colInd,
                                                                               Teuchos::ArrayRCP<scalar_type>& values,
                                                                               const bool debug) {
  using std::cerr;
  using std::endl;
  using Teuchos::arcp;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::as;
  using Teuchos::Comm;
  using Teuchos::CommRequest;
  using Teuchos::null;
  using Teuchos::RCP;
  using Teuchos::receive;
  using Teuchos::send;
 
  const bool extraDebug = false;
  trcp_tcomm_t pComm    = pRowMap->getComm();
  const int numProcs    = pComm->getSize();
  const int myRank      = pComm->getRank();
  const int rootRank    = 0;
 
  // Type abbreviations to make the code more concise.
  typedef global_ordinal_type GO;
 
  // List of the global indices of my rows.  They may or may
  // not be contiguous, and the row map need not be one-to-one.
  ArrayView<const GO> myRows = pRowMap->getLocalElementList();
  const size_type myNumRows  = myRows.size();
  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(myNumRows) !=
                                 pRowMap->getLocalNumElements(),
                             std::logic_error,
                             "pRowMap->getLocalElementList().size() = "
                                 << myNumRows
                                 << " != pRowMap->getLocalNumElements() = "
                                 << pRowMap->getLocalNumElements() << ".  "
                                                                      "Please report this bug to the Tpetra developers.");
  TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
                             std::logic_error,
                             "On Proc 0: numEntriesPerRow.size() = "
                                 << numEntriesPerRow.size()
                                 << " != pRowMap->getLocalElementList().size() = "
                                 << myNumRows << ".  Please report this bug to the "
                                                 "Tpetra developers.");
 
  // Space for my proc's number of entries per row.  Will be
  // filled in below.
  myNumEntriesPerRow = arcp<size_t>(myNumRows);
 
  if (myRank != rootRank) {
    // Tell the root how many rows we have.  If we're sending
    // none, then we don't have anything else to send, nor does
    // the root have to receive anything else.
    send(*pComm, myNumRows, rootRank);
    if (myNumRows != 0) {
      // Now send my rows' global indices.  Hopefully the cast
      // to int doesn't overflow.  This is unlikely, since it
      // should fit in a LO, even though it is a GO.
      send(*pComm, static_cast<int>(myNumRows),
           myRows.getRawPtr(), rootRank);
 
      // I (this proc) don't care if my global row indices are
      // contiguous, though the root proc does (since otherwise
      // it needs to pack noncontiguous data into contiguous
      // storage before sending).  That's why we don't check
      // for contiguousness here.
 
      // Ask the root process for my part of the array of the
      // number of entries per row.
      receive(*pComm, rootRank,
              static_cast<int>(myNumRows),
              myNumEntriesPerRow.getRawPtr());
 
      // Use the resulting array to figure out how many column
      // indices and values I should ask from the root process.
      const size_t myNumEntries =
          std::accumulate(myNumEntriesPerRow.begin(),
                          myNumEntriesPerRow.end(), size_t(0));
      TEUCHOS_TEST_FOR_EXCEPTION(
          myNumEntries > static_cast<size_t>(std::numeric_limits<int>::max()),
          std::overflow_error,
          "MatrixMarketReader::distribute: myNumEntries exceeds int range");
 
      // Make space for my entries of the sparse matrix.  Note
      // that they don't have to be sorted by row index.
      // Iterating through all my rows requires computing a
      // running sum over myNumEntriesPerRow.
      myColInd = arcp<GO>(myNumEntries);
      myValues = arcp<scalar_type>(myNumEntries);
      if (myNumEntries > 0) {
        // Ask for that many column indices and values, if
        // there are any.
        receive(*pComm, rootRank,
                static_cast<int>(myNumEntries),
                myColInd.getRawPtr());
        receive(*pComm, rootRank,
                static_cast<int>(myNumEntries),
                myValues.getRawPtr());
      }
    }     // If I own at least one row
  }       // If I am not the root processor
  else {  // I _am_ the root processor
    if (debug) {
      cerr << "-- Proc 0: Copying my data from global arrays" << endl;
    }
    // Proc 0 still needs to (allocate, if not done already)
    // and fill its part of the matrix (my*).
    for (size_type k = 0; k < myNumRows; ++k) {
      const GO myCurRow                            = myRows[k];
      const local_ordinal_type numEntriesInThisRow = numEntriesPerRow[myCurRow];
      myNumEntriesPerRow[k]                        = numEntriesInThisRow;
    }
    if (extraDebug && debug) {
      cerr << "Proc " << pRowMap->getComm()->getRank()
           << ": myNumEntriesPerRow[0.." << (myNumRows - 1) << "] = [";
      for (size_type k = 0; k < myNumRows; ++k) {
        cerr << myNumEntriesPerRow[k];
        if (k < myNumRows - 1) {
          cerr << " ";
        }
      }
      cerr << "]" << endl;
    }
    // The total number of matrix entries that my proc owns.
    const size_t myNumEntries =
        std::accumulate(myNumEntriesPerRow.begin(),
                        myNumEntriesPerRow.end(), size_t(0));
    TEUCHOS_TEST_FOR_EXCEPTION(
        myNumEntries > static_cast<size_t>(std::numeric_limits<int>::max()),
        std::overflow_error,
        "MatrixMarketReader::distribute: myNumEntries exceeds int range");
    if (debug) {
      cerr << "-- Proc 0: I own " << myNumRows << " rows and "
           << myNumEntries << " entries" << endl;
    }
    myColInd = arcp<GO>(myNumEntries);
    myValues = arcp<scalar_type>(myNumEntries);
 
    // Copy Proc 0's part of the matrix into the my* arrays.
    // It's important that myCurPos be updated _before_ k,
    // otherwise myCurPos will get the wrong number of entries
    // per row (it should be for the row in the just-completed
    // iteration, not for the next iteration's row).
    local_ordinal_type myCurPos = 0;
    for (size_type k = 0; k < myNumRows;
         myCurPos += myNumEntriesPerRow[k], ++k) {
      const local_ordinal_type curNumEntries = myNumEntriesPerRow[k];
      const GO myRow                         = myRows[k];
      const size_t curPos                    = rowPtr[myRow];
      // Only copy if there are entries to copy, in order not
      // to construct empty ranges for the ArrayRCP views.
      if (curNumEntries > 0) {
        ArrayView<GO> colIndView   = colInd(curPos, curNumEntries);
        ArrayView<GO> myColIndView = myColInd(myCurPos, curNumEntries);
        std::copy(colIndView.begin(), colIndView.end(),
                  myColIndView.begin());
 
        ArrayView<scalar_type> valuesView =
            values(curPos, curNumEntries);
        ArrayView<scalar_type> myValuesView =
            myValues(myCurPos, curNumEntries);
        std::copy(valuesView.begin(), valuesView.end(),
                  myValuesView.begin());
      }
    }
 
    // Proc 0 processes each other proc p in turn.
    for (int p = 1; p < numProcs; ++p) {
      if (debug) {
        cerr << "-- Proc 0: Processing proc " << p << endl;
      }
 
      size_type theirNumRows = 0;
      // Ask Proc p how many rows it has.  If it doesn't
      // have any, we can move on to the next proc.  This
      // has to be a standard receive so that we can avoid
      // the degenerate case of sending zero data.
      receive(*pComm, p, &theirNumRows);
      if (debug) {
        cerr << "-- Proc 0: Proc " << p << " owns "
             << theirNumRows << " rows" << endl;
      }
      if (theirNumRows != 0) {
        // Ask Proc p which rows it owns.  The resulting global
        // row indices are not guaranteed to be contiguous or
        // sorted.  Global row indices are themselves indices
        // into the numEntriesPerRow array.
        ArrayRCP<GO> theirRows = arcp<GO>(theirNumRows);
        receive(*pComm, p, as<int>(theirNumRows),
                theirRows.getRawPtr());
        // Extra test to make sure that the rows we received
        // are all sensible.  This is a good idea since we are
        // going to use the global row indices we've received
        // to index into the numEntriesPerRow array.  Better to
        // catch any bugs here and print a sensible error
        // message, rather than segfault and print a cryptic
        // error message.
        {
          const global_size_t numRows = pRowMap->getGlobalNumElements();
          const GO indexBase          = pRowMap->getIndexBase();
          bool theirRowsValid         = true;
          for (size_type k = 0; k < theirNumRows; ++k) {
            if (theirRows[k] < indexBase ||
                as<global_size_t>(theirRows[k] - indexBase) >= numRows) {
              theirRowsValid = false;
            }
          }
          if (!theirRowsValid) {
            TEUCHOS_TEST_FOR_EXCEPTION(
                !theirRowsValid, std::logic_error,
                "Proc " << p << " has at least one invalid row index.  "
                                "Here are all of them: "
                        << Teuchos::toString(theirRows()) << ".  Valid row index "
                                                             "range (zero-based): [0, "
                        << (numRows - 1) << "].");
          }
        }
 
        // Perhaps we could save a little work if we check
        // whether Proc p's row indices are contiguous.  That
        // would make lookups in the global data arrays
        // faster.  For now, we just implement the general
        // case and don't prematurely optimize.  (Remember
        // that you're making Proc 0 read the whole file, so
        // you've already lost scalability.)
 
        // Compute the number of entries in each of Proc p's
        // rows.  (Proc p will compute its row pointer array
        // on its own, after it gets the data from Proc 0.)
        ArrayRCP<size_t> theirNumEntriesPerRow;
        theirNumEntriesPerRow = arcp<size_t>(theirNumRows);
        for (size_type k = 0; k < theirNumRows; ++k) {
          theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
        }
 
        // Tell Proc p the number of entries in each of its
        // rows.  Hopefully the cast to int doesn't overflow.
        // This is unlikely, since it should fit in a LO,
        // even though it is a GO.
        send(*pComm, static_cast<int>(theirNumRows),
             theirNumEntriesPerRow.getRawPtr(), p);
 
        // Figure out how many entries Proc p owns.
        const size_t theirNumEntries =
            std::accumulate(theirNumEntriesPerRow.begin(),
                            theirNumEntriesPerRow.end(), size_t(0));
        TEUCHOS_TEST_FOR_EXCEPTION(
            theirNumEntries > static_cast<size_t>(std::numeric_limits<int>::max()),
            std::overflow_error,
            "MatrixMarketReader::distribute: theirNumEntries exceeds int range");
 
        if (debug) {
          cerr << "-- Proc 0: Proc " << p << " owns "
               << theirNumEntries << " entries" << endl;
        }
 
        // If there are no entries to send, then we're done
        // with Proc p.
        if (theirNumEntries == 0) {
          continue;
        }
 
        // Construct (views of) proc p's column indices and
        // values.  Later, we might like to optimize for the
        // (common) contiguous case, for which we don't need to
        // copy data into separate "their*" arrays (we can just
        // use contiguous views of the global arrays).
        ArrayRCP<GO> theirColInd(theirNumEntries);
        ArrayRCP<scalar_type> theirValues(theirNumEntries);
        // Copy Proc p's part of the matrix into the their*
        // arrays.  It's important that theirCurPos be updated
        // _before_ k, otherwise theirCurPos will get the wrong
        // number of entries per row (it should be for the row
        // in the just-completed iteration, not for the next
        // iteration's row).
        local_ordinal_type theirCurPos = 0;
        for (size_type k = 0; k < theirNumRows;
             theirCurPos += theirNumEntriesPerRow[k], k++) {
          const local_ordinal_type curNumEntries = theirNumEntriesPerRow[k];
          const GO theirRow                      = theirRows[k];
          const local_ordinal_type curPos        = rowPtr[theirRow];
 
          // Only copy if there are entries to copy, in order
          // not to construct empty ranges for the ArrayRCP
          // views.
          if (curNumEntries > 0) {
            ArrayView<GO> colIndView =
                colInd(curPos, curNumEntries);
            ArrayView<GO> theirColIndView =
                theirColInd(theirCurPos, curNumEntries);
            std::copy(colIndView.begin(), colIndView.end(),
                      theirColIndView.begin());
 
            ArrayView<scalar_type> valuesView =
                values(curPos, curNumEntries);
            ArrayView<scalar_type> theirValuesView =
                theirValues(theirCurPos, curNumEntries);
            std::copy(valuesView.begin(), valuesView.end(),
                      theirValuesView.begin());
          }
        }
        // Send Proc p its column indices and values.
        // Hopefully the cast to int doesn't overflow.  This
        // is unlikely, since it should fit in a LO, even
        // though it is a GO.
        send(*pComm, static_cast<int>(theirNumEntries),
             theirColInd.getRawPtr(), p);
        send(*pComm, static_cast<int>(theirNumEntries),
             theirValues.getRawPtr(), p);
 
        if (debug) {
          cerr << "-- Proc 0: Finished with proc " << p << endl;
        }
      }  // If proc p owns at least one row
    }    // For each proc p not the root proc 0
  }      // If I'm (not) the root proc 0
 
  // Invalidate the input data to save space, since we don't
  // need it anymore.
  numEntriesPerRow = null;
  rowPtr           = null;
  colInd           = null;
  values           = null;
 
  if (debug && myRank == 0) {
    cerr << "-- Proc 0: About to fill in myRowPtr" << endl;
  }
 
  // Allocate and fill in myRowPtr (the row pointer array for
  // my rank's rows).  We delay this until the end because we
  // don't need it to compute anything else in distribute().
  // Each proc can do this work for itself, since it only needs
  // myNumEntriesPerRow to do so.
  myRowPtr    = arcp<size_t>(myNumRows + 1);
  myRowPtr[0] = 0;
  for (size_type k = 1; k < myNumRows + 1; ++k) {
    myRowPtr[k] = myRowPtr[k - 1] + myNumEntriesPerRow[k - 1];
  }
  if (extraDebug && debug) {
    cerr << "Proc " << Teuchos::rank(*(pRowMap->getComm()))
         << ": myRowPtr[0.." << myNumRows << "] = [";
    for (size_type k = 0; k < myNumRows + 1; ++k) {
      cerr << myRowPtr[k];
      if (k < myNumRows) {
        cerr << " ";
      }
    }
    cerr << "]" << endl
         << endl;
  }
 
  if (debug && myRank == 0) {
    cerr << "-- Proc 0: Done with distribute" << endl;
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::makeMatrix(Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
                                                                          Teuchos::ArrayRCP<size_t>& myRowPtr,
                                                                          Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
                                                                          Teuchos::ArrayRCP<scalar_type>& myValues,
                                                                          const Teuchos::RCP<const map_type>& pRowMap,
                                                                          const Teuchos::RCP<const map_type>& pRangeMap,
                                                                          const Teuchos::RCP<const map_type>& pDomainMap,
                                                                          const bool callFillComplete) {
  using std::cerr;
  using std::endl;
  using Teuchos::ArrayView;
  using Teuchos::null;
  using Teuchos::RCP;
  using Teuchos::rcp;
  // Typedef to make certain type declarations shorter.
  typedef global_ordinal_type GO;
 
  // The row pointer array always has at least one entry, even
  // if the matrix has zero rows.  myNumEntriesPerRow, myColInd,
  // and myValues would all be empty arrays in that degenerate
  // case, but the row and domain maps would still be nonnull
  // (though they would be trivial maps).
  TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
                             "makeMatrix: myRowPtr array is null.  "
                             "Please report this bug to the Tpetra developers.");
  TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
                             "makeMatrix: domain map is null.  "
                             "Please report this bug to the Tpetra developers.");
  TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
                             "makeMatrix: range map is null.  "
                             "Please report this bug to the Tpetra developers.");
  TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
                             "makeMatrix: row map is null.  "
                             "Please report this bug to the Tpetra developers.");
 
  // Construct the CrsMatrix, using the row map, with the
  // constructor specifying the number of nonzeros for each row.
  RCP<sparse_matrix_type> A =
      rcp(new sparse_matrix_type(pRowMap, myNumEntriesPerRow()));
 
  // List of the global indices of my rows.
  // They may or may not be contiguous.
  ArrayView<const GO> myRows = pRowMap->getLocalElementList();
  const size_type myNumRows  = myRows.size();
 
  // Add this processor's matrix entries to the CrsMatrix.
  const GO indexBase = pRowMap->getIndexBase();
  for (size_type i = 0; i < myNumRows; ++i) {
    const size_type myCurPos               = myRowPtr[i];
    const local_ordinal_type curNumEntries = myNumEntriesPerRow[i];
    ArrayView<GO> curColInd                = myColInd.view(myCurPos, curNumEntries);
    ArrayView<scalar_type> curValues       = myValues.view(myCurPos, curNumEntries);
 
    // Modify the column indices in place to have the right index base.
    for (size_type k = 0; k < curNumEntries; ++k) {
      curColInd[k] += indexBase;
    }
    // Avoid constructing empty views of ArrayRCP objects.
    if (curNumEntries > 0) {
      A->insertGlobalValues(myRows[i], curColInd, curValues);
    }
  }
  // We've entered in all our matrix entries, so we can delete
  // the original data.  This will save memory when we call
  // fillComplete(), so that we never keep more than two copies
  // of the matrix's data in memory at once.
  myNumEntriesPerRow = null;
  myRowPtr           = null;
  myColInd           = null;
  myValues           = null;
 
  if (callFillComplete) {
    A->fillComplete(pDomainMap, pRangeMap);
  }
  return A;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::makeMatrix(Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
                                                                          Teuchos::ArrayRCP<size_t>& myRowPtr,
                                                                          Teuchos::ArrayRCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_ordinal_type>& myColInd,
                                                                          Teuchos::ArrayRCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::scalar_type>& myValues,
                                                                          const Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& pRowMap,
                                                                          const Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& pRangeMap,
                                                                          const Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& pDomainMap,
                                                                          const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
                                                                          const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams) {
  using std::cerr;
  using std::endl;
  using Teuchos::ArrayView;
  using Teuchos::null;
  using Teuchos::RCP;
  using Teuchos::rcp;
  // Typedef to make certain type declarations shorter.
  typedef typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_ordinal_type GO;
  using size_type = typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::size_type;
 
  // The row pointer array always has at least one entry, even
  // if the matrix has zero rows.  myNumEntriesPerRow, myColInd,
  // and myValues would all be empty arrays in that degenerate
  // case, but the row and domain maps would still be nonnull
  // (though they would be trivial maps).
  TEUCHOS_TEST_FOR_EXCEPTION(
      myRowPtr.is_null(), std::logic_error,
      "makeMatrix: myRowPtr array is null.  "
      "Please report this bug to the Tpetra developers.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      pDomainMap.is_null(), std::logic_error,
      "makeMatrix: domain map is null.  "
      "Please report this bug to the Tpetra developers.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      pRangeMap.is_null(), std::logic_error,
      "makeMatrix: range map is null.  "
      "Please report this bug to the Tpetra developers.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      pRowMap.is_null(), std::logic_error,
      "makeMatrix: row map is null.  "
      "Please report this bug to the Tpetra developers.");
 
  // Construct the CrsMatrix, using the row map, with the
  // constructor specifying the number of nonzeros for each row.
  RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type> A =
      rcp(new typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type(pRowMap, myNumEntriesPerRow(),
                                                                                                         constructorParams));
 
  // List of the global indices of my rows.
  // They may or may not be contiguous.
  ArrayView<const GO> myRows = pRowMap->getLocalElementList();
  const size_type myNumRows  = myRows.size();
 
  // Add this processor's matrix entries to the CrsMatrix.
  const GO indexBase = pRowMap->getIndexBase();
  for (size_type i = 0; i < myNumRows; ++i) {
    const size_type myCurPos                                                                                       = myRowPtr[i];
    const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_ordinal_type curNumEntries = myNumEntriesPerRow[i];
    ArrayView<GO> curColInd                                                                                        = myColInd.view(myCurPos, curNumEntries);
    ArrayView<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::scalar_type> curValues       = myValues.view(myCurPos, curNumEntries);
 
    // Modify the column indices in place to have the right index base.
    for (size_type k = 0; k < curNumEntries; ++k) {
      curColInd[k] += indexBase;
    }
    if (curNumEntries > 0) {
      A->insertGlobalValues(myRows[i], curColInd, curValues);
    }
  }
  // We've entered in all our matrix entries, so we can delete
  // the original data.  This will save memory when we call
  // fillComplete(), so that we never keep more than two copies
  // of the matrix's data in memory at once.
  myNumEntriesPerRow = null;
  myRowPtr           = null;
  myColInd           = null;
  myValues           = null;
 
  A->fillComplete(pDomainMap, pRangeMap, fillCompleteParams);
  return A;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::makeMatrix(Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
                                                                          Teuchos::ArrayRCP<size_t>& myRowPtr,
                                                                          Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
                                                                          Teuchos::ArrayRCP<scalar_type>& myValues,
                                                                          const Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& rowMap,
                                                                          Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& colMap,
                                                                          const Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& domainMap,
                                                                          const Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& rangeMap,
                                                                          const bool callFillComplete) {
  using Teuchos::ArrayView;
  using Teuchos::as;
  using Teuchos::null;
  using Teuchos::RCP;
  using Teuchos::rcp;
  typedef global_ordinal_type GO;
 
  // Construct the CrsMatrix.
 
  RCP<sparse_matrix_type> A;  // the matrix to return.
  if (colMap.is_null()) {     // the user didn't provide a column Map
    A = rcp(new sparse_matrix_type(rowMap, myNumEntriesPerRow()));
  } else {  // the user provided a column Map
    A = rcp(new sparse_matrix_type(rowMap, colMap, myNumEntriesPerRow()));
  }
 
  // List of the global indices of my rows.
  // They may or may not be contiguous.
  ArrayView<const GO> myRows = rowMap->getLocalElementList();
  const size_type myNumRows  = myRows.size();
 
  // Add this process' matrix entries to the CrsMatrix.
  const GO indexBase = rowMap->getIndexBase();
  for (size_type i = 0; i < myNumRows; ++i) {
    const size_type myCurPos         = myRowPtr[i];
    const size_type curNumEntries    = as<size_type>(myNumEntriesPerRow[i]);
    ArrayView<GO> curColInd          = myColInd.view(myCurPos, curNumEntries);
    ArrayView<scalar_type> curValues = myValues.view(myCurPos, curNumEntries);
 
    // Modify the column indices in place to have the right index base.
    for (size_type k = 0; k < curNumEntries; ++k) {
      curColInd[k] += indexBase;
    }
    if (curNumEntries > 0) {
      A->insertGlobalValues(myRows[i], curColInd, curValues);
    }
  }
  // We've entered in all our matrix entries, so we can delete
  // the original data.  This will save memory when we call
  // fillComplete(), so that we never keep more than two copies
  // of the matrix's data in memory at once.
  myNumEntriesPerRow = null;
  myRowPtr           = null;
  myColInd           = null;
  myValues           = null;
 
  if (callFillComplete) {
    A->fillComplete(domainMap, rangeMap);
    if (colMap.is_null()) {
      colMap = A->getColMap();
    }
  }
  return A;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Teuchos::MatrixMarket::Banner>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readBanner(std::istream& in,
                                                                          size_t& lineNumber,
                                                                          const bool tolerant,
                                                                          const bool /* debug */,
                                                                          const bool isGraph) {
  using std::cerr;
  using std::endl;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::MatrixMarket::Banner;
  typedef Teuchos::ScalarTraits<scalar_type> STS;
 
  RCP<Banner> pBanner;  // On output, if successful: the read-in Banner.
  std::string line;     // If read from stream successful: the Banner line
 
  // Try to read a line from the input stream.
  const bool readFailed = !getline(in, line);
  TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
                             "Failed to get Matrix Market banner line from input.");
 
  // We read a line from the input stream.
  lineNumber++;
 
  // Assume that the line we found is the Banner line.
  try {
    pBanner = rcp(new Banner(line, tolerant));
  } catch (std::exception& e) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
                               "Matrix Market banner line contains syntax error(s): "
                                   << e.what());
  }
 
  TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() != "matrix",
                             std::invalid_argument,
                             "The Matrix Market file does not contain "
                             "matrix data.  Its Banner (first) line says that its object type is \""
                                 << pBanner->matrixType() << "\", rather than the required \"matrix\".");
 
  // Validate the data type of the matrix, with respect to the
  // Scalar type of the CrsMatrix entries.
  TEUCHOS_TEST_FOR_EXCEPTION(
      !STS::isComplex && pBanner->dataType() == "complex",
      std::invalid_argument,
      "The Matrix Market file contains complex-valued data, but you are "
      "trying to read it into a matrix containing entries of the real-"
      "valued Scalar type \""
          << Teuchos::TypeNameTraits<scalar_type>::name() << "\".");
  TEUCHOS_TEST_FOR_EXCEPTION(
      !isGraph &&
          pBanner->dataType() != "real" &&
          pBanner->dataType() != "complex" &&
          pBanner->dataType() != "integer",
      std::invalid_argument,
      "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
      "Matrix Market file may not contain a \"pattern\" matrix.  A "
      "pattern matrix is really just a graph with no weights.  It "
      "should be stored in a CrsGraph, not a CrsMatrix.");
 
  TEUCHOS_TEST_FOR_EXCEPTION(
      isGraph &&
          pBanner->dataType() != "pattern",
      std::invalid_argument,
      "When reading Matrix Market data into a Tpetra::CrsGraph, the "
      "Matrix Market file must contain a \"pattern\" matrix.");
 
  return pBanner;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::Tuple<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_ordinal_type, 3>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readCoordDims(std::istream& in,
                                                                             size_t& lineNumber,
                                                                             const Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
                                                                             const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t& pComm,
                                                                             const bool tolerant,
                                                                             const bool /* debug */) {
  using Teuchos::Tuple;
  using Teuchos::MatrixMarket::readCoordinateDimensions;
 
  // Packed coordinate matrix dimensions (numRows, numCols,
  // numNonzeros); computed on Rank 0 and broadcasted to all MPI
  // ranks.
  Tuple<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_ordinal_type, 3> dims;
 
  // Read in the coordinate matrix dimensions from the input
  // stream.  "success" tells us whether reading in the
  // coordinate matrix dimensions succeeded ("Guilty unless
  // proven innocent").
  bool success = false;
  if (pComm->getRank() == 0) {
    TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() != "coordinate",
                               std::invalid_argument,
                               "The Tpetra::CrsMatrix Matrix Market reader "
                               "only accepts \"coordinate\" (sparse) matrix data.");
    // Unpacked coordinate matrix dimensions
    typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_ordinal_type numRows, numCols, numNonzeros;
    // Only MPI Rank 0 reads from the input stream
    success = readCoordinateDimensions(in, numRows, numCols,
                                       numNonzeros, lineNumber,
                                       tolerant);
    // Pack up the data into a Tuple so we can send them with
    // one broadcast instead of three.
    dims[0] = numRows;
    dims[1] = numCols;
    dims[2] = numNonzeros;
  }
  // Only Rank 0 did the reading, so it decides success.
  //
  // FIXME (mfh 02 Feb 2011) Teuchos::broadcast doesn't know how
  // to send bools.  For now, we convert to/from int instead,
  // using the usual "true is 1, false is 0" encoding.
  {
    int the_success = success ? 1 : 0;  // only matters on MPI Rank 0
    Teuchos::broadcast(*pComm, 0, &the_success);
    success = (the_success == 1);
  }
  if (success) {
    // Broadcast (numRows, numCols, numNonzeros) from Rank 0
    // to all the other MPI ranks.
    Teuchos::broadcast(*pComm, 0, dims);
  } else {
    // Perhaps in tolerant mode, we could set all the
    // dimensions to zero for now, and deduce correct
    // dimensions by reading all of the file's entries and
    // computing the max(row index) and max(column index).
    // However, for now we just error out in that case.
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
                               "Error reading Matrix Market sparse matrix: failed to read "
                               "coordinate matrix dimensions.");
  }
  return dims;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::adder_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::makeAdder(const Teuchos::RCP<const Teuchos::Comm<int>>& pComm,
                                                                         Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
                                                                         const Teuchos::Tuple<global_ordinal_type, 3>& dims,
                                                                         const bool tolerant,
                                                                         const bool debug) {
  if (pComm->getRank() == 0) {
    typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type,
                                              global_ordinal_type>
        raw_adder_type;
    Teuchos::RCP<raw_adder_type> pRaw =
        Teuchos::rcp(new raw_adder_type(dims[0], dims[1], dims[2],
                                        tolerant, debug));
    return Teuchos::rcp(new adder_type(pRaw, pBanner->symmType()));
  } else {
    return Teuchos::null;
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::graph_adder_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::makeGraphAdder(const Teuchos::RCP<const Teuchos::Comm<int>>& pComm,
                                                                              Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
                                                                              const Teuchos::Tuple<global_ordinal_type, 3>& dims,
                                                                              const bool tolerant,
                                                                              const bool debug) {
  if (pComm->getRank() == 0) {
    typedef Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> raw_adder_type;
    Teuchos::RCP<raw_adder_type> pRaw =
        Teuchos::rcp(new raw_adder_type(dims[0], dims[1], dims[2],
                                        tolerant, debug));
    return Teuchos::rcp(new graph_adder_type(pRaw, pBanner->symmType()));
  } else {
    return Teuchos::null;
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_graph_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparseGraphHelper(std::istream& in,
                                                                                     const Teuchos::RCP<const Teuchos::Comm<int>>& pComm,
                                                                                     const Teuchos::RCP<const map_type>& rowMap,
                                                                                     Teuchos::RCP<const map_type>& colMap,
                                                                                     const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
                                                                                     const bool tolerant,
                                                                                     const bool debug) {
  using std::cerr;
  using std::endl;
  using Teuchos::ptr;
  using Teuchos::RCP;
  using Teuchos::Tuple;
  using Teuchos::MatrixMarket::Banner;
 
  const int myRank   = pComm->getRank();
  const int rootRank = 0;
 
  // Current line number in the input stream.  Various calls
  // will modify this depending on the number of lines that are
  // read from the input stream.  Only Rank 0 modifies this.
  size_t lineNumber = 1;
 
  if (debug && myRank == rootRank) {
    cerr << "Matrix Market reader: readGraph:" << endl
         << "-- Reading banner line" << endl;
  }
 
  // The "Banner" tells you whether the input stream represents
  // a sparse matrix, the symmetry type of the matrix, and the
  // type of the data it contains.
  //
  // pBanner will only be nonnull on MPI Rank 0.  It will be
  // null on all other MPI processes.
  RCP<const Banner> pBanner;
  {
    // We read and validate the Banner on Proc 0, but broadcast
    // the validation result to all processes.
    // Teuchos::broadcast doesn't currently work with bool, so
    // we use int (true -> 1, false -> 0).
    int bannerIsCorrect = 1;
    std::ostringstream errMsg;
 
    if (myRank == rootRank) {
      // Read the Banner line from the input stream.
      try {
        pBanner = readBanner(in, lineNumber, tolerant, debug, true);
      } catch (std::exception& e) {
        errMsg << "Attempt to read the Matrix Market file's Banner line "
                  "threw an exception: "
               << e.what();
        bannerIsCorrect = 0;
      }
 
      if (bannerIsCorrect) {
        // Validate the Banner for the case of a sparse graph.
        // We validate on Proc 0, since it reads the Banner.
 
        // In intolerant mode, the matrix type must be "coordinate".
        if (!tolerant && pBanner->matrixType() != "coordinate") {
          bannerIsCorrect = 0;
          errMsg << "The Matrix Market input file must contain a "
                    "\"coordinate\"-format sparse graph in order to create a "
                    "Tpetra::CrsGraph object from it, but the file's matrix "
                    "type is \""
                 << pBanner->matrixType() << "\" instead.";
        }
        // In tolerant mode, we allow the matrix type to be
        // anything other than "array" (which would mean that
        // the file contains a dense matrix).
        if (tolerant && pBanner->matrixType() == "array") {
          bannerIsCorrect = 0;
          errMsg << "Matrix Market file must contain a \"coordinate\"-"
                    "format sparse graph in order to create a Tpetra::CrsGraph "
                    "object from it, but the file's matrix type is \"array\" "
                    "instead.  That probably means the file contains dense matrix "
                    "data.";
        }
      }
    }  // Proc 0: Done reading the Banner, hopefully successfully.
 
    // Broadcast from Proc 0 whether the Banner was read correctly.
    broadcast(*pComm, rootRank, ptr(&bannerIsCorrect));
 
    // If the Banner is invalid, all processes throw an
    // exception.  Only Proc 0 gets the exception message, but
    // that's OK, since the main point is to "stop the world"
    // (rather than throw an exception on one process and leave
    // the others hanging).
    TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
                               std::invalid_argument, errMsg.str());
  }  // Done reading the Banner line and broadcasting success.
  if (debug && myRank == rootRank) {
    cerr << "-- Reading dimensions line" << endl;
  }
 
  // Read the graph dimensions from the Matrix Market metadata.
  // dims = (numRows, numCols, numEntries).  Proc 0 does the
  // reading, but it broadcasts the results to all MPI
  // processes.  Thus, readCoordDims() is a collective
  // operation.  It does a collective check for correctness too.
  Tuple<global_ordinal_type, 3> dims =
      readCoordDims(in, lineNumber, pBanner, pComm, tolerant, debug);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Making Adder for collecting graph data" << endl;
  }
 
  // "Adder" object for collecting all the sparse graph entries
  // from the input stream.  This is only nonnull on Proc 0.
  // The Adder internally converts the one-based indices (native
  // Matrix Market format) into zero-based indices.
  RCP<graph_adder_type> pAdder =
      makeGraphAdder(pComm, pBanner, dims, tolerant, debug);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Reading graph data" << endl;
  }
  //
  // Read the graph entries from the input stream on Proc 0.
  //
  {
    // We use readSuccess to broadcast the results of the read
    // (succeeded or not) to all MPI processes.  Since
    // Teuchos::broadcast doesn't currently know how to send
    // bools, we convert to int (true -> 1, false -> 0).
    int readSuccess = 1;
    std::ostringstream errMsg;  // Exception message (only valid on Proc 0)
    if (myRank == rootRank) {
      try {
        // Reader for "coordinate" format sparse graph data.
        typedef Teuchos::MatrixMarket::CoordPatternReader<graph_adder_type,
                                                          global_ordinal_type>
            reader_type;
        reader_type reader(pAdder);
 
        // Read the sparse graph entries.
        std::pair<bool, std::vector<size_t>> results =
            reader.read(in, lineNumber, tolerant, debug);
        readSuccess = results.first ? 1 : 0;
      } catch (std::exception& e) {
        readSuccess = 0;
        errMsg << e.what();
      }
    }
    broadcast(*pComm, rootRank, ptr(&readSuccess));
 
    // It would be nice to add a "verbose" flag, so that in
    // tolerant mode, we could log any bad line number(s) on
    // Proc 0.  For now, we just throw if the read fails to
    // succeed.
    //
    // Question: If we're in tolerant mode, and if the read did
    // not succeed, should we attempt to call fillComplete()?
    TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
                               "Failed to read the Matrix Market sparse graph file: "
                                   << errMsg.str());
  }  // Done reading the graph entries (stored on Proc 0 for now)
 
  if (debug && myRank == rootRank) {
    cerr << "-- Successfully read the Matrix Market data" << endl;
  }
 
  // In tolerant mode, we need to rebroadcast the graph
  // dimensions, since they may be different after reading the
  // actual graph data.  We only need to broadcast the number
  // of rows and columns.  Only Rank 0 needs to know the actual
  // global number of entries, since (a) we need to merge
  // duplicates on Rank 0 first anyway, and (b) when we
  // distribute the entries, each rank other than Rank 0 will
  // only need to know how many entries it owns, not the total
  // number of entries.
  if (tolerant) {
    if (debug && myRank == rootRank) {
      cerr << "-- Tolerant mode: rebroadcasting graph dimensions"
           << endl
           << "----- Dimensions before: "
           << dims[0] << " x " << dims[1]
           << endl;
    }
    // Packed coordinate graph dimensions (numRows, numCols).
    Tuple<global_ordinal_type, 2> updatedDims;
    if (myRank == rootRank) {
      // If one or more bottom rows of the graph contain no
      // entries, then the Adder will report that the number
      // of rows is less than that specified in the
      // metadata.  We allow this case, and favor the
      // metadata so that the zero row(s) will be included.
      updatedDims[0] =
          std::max(dims[0], pAdder->getAdder()->numRows());
      updatedDims[1] = pAdder->getAdder()->numCols();
    }
    broadcast(*pComm, rootRank, updatedDims);
    dims[0] = updatedDims[0];
    dims[1] = updatedDims[1];
    if (debug && myRank == rootRank) {
      cerr << "----- Dimensions after: " << dims[0] << " x "
           << dims[1] << endl;
    }
  } else {
    // In strict mode, we require that the graph's metadata and
    // its actual data agree, at least somewhat.  In particular,
    // the number of rows must agree, since otherwise we cannot
    // distribute the graph correctly.
 
    // Teuchos::broadcast() doesn't know how to broadcast bools,
    // so we use an int with the standard 1 == true, 0 == false
    // encoding.
    int dimsMatch = 1;
    if (myRank == rootRank) {
      // If one or more bottom rows of the graph contain no
      // entries, then the Adder will report that the number of
      // rows is less than that specified in the metadata.  We
      // allow this case, and favor the metadata, but do not
      // allow the Adder to think there are more rows in the
      // graph than the metadata says.
      if (dims[0] < pAdder->getAdder()->numRows()) {
        dimsMatch = 0;
      }
    }
    broadcast(*pComm, 0, ptr(&dimsMatch));
    if (dimsMatch == 0) {
      // We're in an error state anyway, so we might as well
      // work a little harder to print an informative error
      // message.
      //
      // Broadcast the Adder's idea of the graph dimensions
      // from Proc 0 to all processes.
      Tuple<global_ordinal_type, 2> addersDims;
      if (myRank == rootRank) {
        addersDims[0] = pAdder->getAdder()->numRows();
        addersDims[1] = pAdder->getAdder()->numCols();
      }
      broadcast(*pComm, 0, addersDims);
      TEUCHOS_TEST_FOR_EXCEPTION(
          dimsMatch == 0, std::runtime_error,
          "The graph metadata says that the graph is " << dims[0] << " x "
                                                       << dims[1] << ", but the actual data says that the graph is "
                                                       << addersDims[0] << " x " << addersDims[1] << ".  That means the "
                                                                                                     "data includes more rows than reported in the metadata.  This "
                                                                                                     "is not allowed when parsing in strict mode.  Parse the graph in "
                                                                                                     "tolerant mode to ignore the metadata when it disagrees with the "
                                                                                                     "data.");
    }
  }  // Matrix dimensions (# rows, # cols, # entries) agree.
 
  // Create a map describing a distribution where the root owns EVERYTHING
  RCP<map_type> proc0Map;
  global_ordinal_type indexBase;
  if (Teuchos::is_null(rowMap)) {
    indexBase = 0;
  } else {
    indexBase = rowMap->getIndexBase();
  }
  if (myRank == rootRank) {
    proc0Map = rcp(new map_type(dims[0], dims[0], indexBase, pComm));
  } else {
    proc0Map = rcp(new map_type(dims[0], 0, indexBase, pComm));
  }
 
  // Create the graph where the root owns EVERYTHING
  std::map<global_ordinal_type, size_t> numEntriesPerRow_map;
  if (myRank == rootRank) {
    const auto& entries = pAdder()->getAdder()->getEntries();
    // This will count duplicates, but it's better than dense.
    // An even better approach would use a classic algorithm,
    // likely in Saad's old textbook, for converting COO (entries)
    // to CSR (the local part of the sparse matrix data structure).
    for (const auto& entry : entries) {
      const global_ordinal_type gblRow = entry.rowIndex() + indexBase;
      ++numEntriesPerRow_map[gblRow];
    }
  }
 
  Teuchos::Array<size_t> numEntriesPerRow(proc0Map->getLocalNumElements());
  for (const auto& ent : numEntriesPerRow_map) {
    const local_ordinal_type lclRow = proc0Map->getLocalElement(ent.first);
    TEUCHOS_TEST_FOR_EXCEPTION(
        lclRow == Teuchos::OrdinalTraits<local_ordinal_type>::invalid(),
        std::runtime_error,
        "MatrixMarketReader::readSparseGraphHelper: global row "
            << ent.first << " is not owned by proc0Map.  This should not happen.");
    numEntriesPerRow[lclRow] = ent.second;
  }
  // Free anything we don't need before allocating the graph.
  // Swapping with an empty data structure is the standard idiom
  // for freeing memory used by Standard Library containers.
  // (Just resizing to 0 doesn't promise to free memory.)
  {
    std::map<global_ordinal_type, size_t> empty_map;
    std::swap(numEntriesPerRow_map, empty_map);
  }
 
  RCP<sparse_graph_type> proc0Graph =
      rcp(new sparse_graph_type(proc0Map, numEntriesPerRow(),
                                constructorParams));
  if (myRank == rootRank) {
    typedef Teuchos::MatrixMarket::Raw::GraphElement<global_ordinal_type> element_type;
 
    // Get the entries
    const std::vector<element_type>& entries =
        pAdder->getAdder()->getEntries();
 
    // Insert them one at a time
    for (size_t curPos = 0; curPos < entries.size(); curPos++) {
      const element_type& curEntry     = entries[curPos];
      const global_ordinal_type curRow = curEntry.rowIndex() + indexBase;
      const global_ordinal_type curCol = curEntry.colIndex() + indexBase;
      Teuchos::ArrayView<const global_ordinal_type> colView(&curCol, 1);
      proc0Graph->insertGlobalIndices(curRow, colView);
    }
  }
  proc0Graph->fillComplete();
 
  RCP<sparse_graph_type> distGraph;
  if (Teuchos::is_null(rowMap)) {
    // Create a map describing the distribution we actually want
    RCP<map_type> distMap =
        rcp(new map_type(dims[0], 0, pComm, GloballyDistributed));
 
    // Create the graph with that distribution too
    distGraph = rcp(new sparse_graph_type(distMap, colMap, 0, constructorParams));
 
    // Create an importer/exporter/vandelay to redistribute the graph
    typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
    import_type importer(proc0Map, distMap);
 
    // Import the data
    distGraph->doImport(*proc0Graph, importer, INSERT);
  } else {
    distGraph = rcp(new sparse_graph_type(rowMap, colMap, 0, constructorParams));
 
    // Create an importer/exporter/vandelay to redistribute the graph
    typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
    import_type importer(proc0Map, rowMap);
 
    // Import the data
    distGraph->doImport(*proc0Graph, importer, INSERT);
  }
 
  return distGraph;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_graph_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparseGraphFile(const std::string& filename,
                                                                                   const trcp_tcomm_t& comm,
                                                                                   const bool callFillComplete,
                                                                                   const bool tolerant,
                                                                                   const bool debug) {
  std::ifstream in = MatrixMarketReader::openInFileOnRankZero(comm, filename, true);
 
  return readSparseGraph(in, comm,
                         callFillComplete,
                         tolerant, debug);
  // We can rely on the destructor of the input stream to close
  // the file on scope exit, even if readSparseGraph() throws an
  // exception.
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
                 sparse_graph_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparseGraphFile(const std::string& filename,
                                                                                   const Teuchos::RCP<const Teuchos::Comm<int>>& pComm,
                                                                                   const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
                                                                                   const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
                                                                                   const bool tolerant,
                                                                                   const bool debug) {
  std::ifstream in = MatrixMarketReader::openInFileOnRankZero(pComm, filename, true);
 
  return readSparseGraph(in, pComm,
                         constructorParams,
                         fillCompleteParams, tolerant, debug);
  // We can rely on the destructor of the input stream to close
  // the file on scope exit, even if readSparseGraph() throws an
  // exception.
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_graph_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparseGraphFile(const std::string& filename,
                                                                                   const Teuchos::RCP<const map_type>& rowMap,
                                                                                   Teuchos::RCP<const map_type>& colMap,
                                                                                   const Teuchos::RCP<const map_type>& domainMap,
                                                                                   const Teuchos::RCP<const map_type>& rangeMap,
                                                                                   const bool callFillComplete,
                                                                                   const bool tolerant,
                                                                                   const bool debug) {
  TEUCHOS_TEST_FOR_EXCEPTION(rowMap.is_null(), std::invalid_argument,
                             "Input rowMap must be nonnull.");
  trcp_tcomm_t comm = rowMap->getComm();
  if (comm.is_null()) {
    // If the input communicator is null on some process, then
    // that process does not participate in the collective.
    return Teuchos::null;
  }
 
  std::ifstream in = MatrixMarketReader::openInFileOnRankZero(comm, filename, true);
 
  return readSparseGraph(in, rowMap, colMap, domainMap, rangeMap,
                         callFillComplete, tolerant, debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_graph_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparseGraph(std::istream& in,
                                                                               const Teuchos::RCP<const Teuchos::Comm<int>>& pComm,
                                                                               const bool callFillComplete,
                                                                               const bool tolerant,
                                                                               const bool debug) {
  Teuchos::RCP<const map_type> fakeRowMap;
  Teuchos::RCP<const map_type> fakeColMap;
  Teuchos::RCP<Teuchos::ParameterList> fakeCtorParams;
 
  Teuchos::RCP<sparse_graph_type> graph =
      readSparseGraphHelper(in, pComm,
                            fakeRowMap, fakeColMap,
                            fakeCtorParams, tolerant, debug);
  if (callFillComplete) {
    graph->fillComplete();
  }
  return graph;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_graph_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparseGraph(std::istream& in,
                                                                               const Teuchos::RCP<const Teuchos::Comm<int>>& pComm,
                                                                               const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
                                                                               const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
                                                                               const bool tolerant,
                                                                               const bool debug) {
  Teuchos::RCP<const map_type> fakeRowMap;
  Teuchos::RCP<const map_type> fakeColMap;
  Teuchos::RCP<sparse_graph_type> graph =
      readSparseGraphHelper(in, pComm,
                            fakeRowMap, fakeColMap,
                            constructorParams, tolerant, debug);
  graph->fillComplete(fillCompleteParams);
  return graph;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_graph_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparseGraph(std::istream& in,
                                                                               const Teuchos::RCP<const map_type>& rowMap,
                                                                               Teuchos::RCP<const map_type>& colMap,
                                                                               const Teuchos::RCP<const map_type>& domainMap,
                                                                               const Teuchos::RCP<const map_type>& rangeMap,
                                                                               const bool callFillComplete,
                                                                               const bool tolerant,
                                                                               const bool debug) {
  Teuchos::RCP<sparse_graph_type> graph =
      readSparseGraphHelper(in, rowMap->getComm(),
                            rowMap, colMap, Teuchos::null, tolerant,
                            debug);
  if (callFillComplete) {
    graph->fillComplete(domainMap, rangeMap);
  }
  return graph;
}
 
#include "MatrixMarket_TpetraNew.hpp"
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparseFile(const std::string& filename,
                                                                              const trcp_tcomm_t& comm,
                                                                              const bool callFillComplete,
                                                                              const bool tolerant,
                                                                              const bool debug) {
  std::ifstream in = MatrixMarketReader::openInFileOnRankZero(comm, filename, true);
 
  return readSparse(in, comm, callFillComplete, tolerant, debug);
  // We can rely on the destructor of the input stream to close
  // the file on scope exit, even if readSparse() throws an
  // exception.
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparseFile(const std::string& filename,
                                                                              const trcp_tcomm_t& comm,
                                                                              const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
                                                                              const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
                                                                              const bool tolerant,
                                                                              const bool debug) {
  std::ifstream in = MatrixMarketReader::openInFileOnRankZero(comm, filename, true);
 
  return readSparse(in, comm, constructorParams,
                    fillCompleteParams, tolerant, debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparseFile(const std::string& filename,
                                                                              const Teuchos::RCP<const map_type>& rowMap,
                                                                              Teuchos::RCP<const map_type>& colMap,
                                                                              const Teuchos::RCP<const map_type>& domainMap,
                                                                              const Teuchos::RCP<const map_type>& rangeMap,
                                                                              const bool callFillComplete,
                                                                              const bool tolerant,
                                                                              const bool debug) {
  TEUCHOS_TEST_FOR_EXCEPTION(
      rowMap.is_null(), std::invalid_argument,
      "Row Map must be nonnull.");
 
  trcp_tcomm_t comm = rowMap->getComm();
 
  std::ifstream in = MatrixMarketReader::openInFileOnRankZero(comm, filename, true);
 
  return readSparse(in, rowMap, colMap, domainMap, rangeMap,
                    callFillComplete, tolerant, debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparse(std::istream& in,
                                                                          const Teuchos::RCP<const Teuchos::Comm<int>>& pComm,
                                                                          const bool callFillComplete,
                                                                          const bool tolerant,
                                                                          const bool debug) {
  using std::cerr;
  using std::endl;
  using Teuchos::arcp;
  using Teuchos::ArrayRCP;
  using Teuchos::broadcast;
  using Teuchos::null;
  using Teuchos::ptr;
  using Teuchos::RCP;
  using Teuchos::REDUCE_MAX;
  using Teuchos::reduceAll;
  using Teuchos::Tuple;
  using Teuchos::MatrixMarket::Banner;
  typedef Teuchos::ScalarTraits<scalar_type> STS;
 
  const bool extraDebug = false;
  const int myRank      = pComm->getRank();
  const int rootRank    = 0;
 
  // Current line number in the input stream.  Various calls
  // will modify this depending on the number of lines that are
  // read from the input stream.  Only Rank 0 modifies this.
  size_t lineNumber = 1;
 
  if (debug && myRank == rootRank) {
    cerr << "Matrix Market reader: readSparse:" << endl
         << "-- Reading banner line" << endl;
  }
 
  // The "Banner" tells you whether the input stream represents
  // a sparse matrix, the symmetry type of the matrix, and the
  // type of the data it contains.
  //
  // pBanner will only be nonnull on MPI Rank 0.  It will be
  // null on all other MPI processes.
  RCP<const Banner> pBanner;
  {
    // We read and validate the Banner on Proc 0, but broadcast
    // the validation result to all processes.
    // Teuchos::broadcast doesn't currently work with bool, so
    // we use int (true -> 1, false -> 0).
    int bannerIsCorrect = 1;
    std::ostringstream errMsg;
 
    if (myRank == rootRank) {
      // Read the Banner line from the input stream.
      try {
        pBanner = readBanner(in, lineNumber, tolerant, debug);
      } catch (std::exception& e) {
        errMsg << "Attempt to read the Matrix Market file's Banner line "
                  "threw an exception: "
               << e.what();
        bannerIsCorrect = 0;
      }
 
      if (bannerIsCorrect) {
        // Validate the Banner for the case of a sparse matrix.
        // We validate on Proc 0, since it reads the Banner.
 
        // In intolerant mode, the matrix type must be "coordinate".
        if (!tolerant && pBanner->matrixType() != "coordinate") {
          bannerIsCorrect = 0;
          errMsg << "The Matrix Market input file must contain a "
                    "\"coordinate\"-format sparse matrix in order to create a "
                    "Tpetra::CrsMatrix object from it, but the file's matrix "
                    "type is \""
                 << pBanner->matrixType() << "\" instead.";
        }
        // In tolerant mode, we allow the matrix type to be
        // anything other than "array" (which would mean that
        // the file contains a dense matrix).
        if (tolerant && pBanner->matrixType() == "array") {
          bannerIsCorrect = 0;
          errMsg << "Matrix Market file must contain a \"coordinate\"-"
                    "format sparse matrix in order to create a Tpetra::CrsMatrix "
                    "object from it, but the file's matrix type is \"array\" "
                    "instead.  That probably means the file contains dense matrix "
                    "data.";
        }
      }
    }  // Proc 0: Done reading the Banner, hopefully successfully.
 
    // Broadcast from Proc 0 whether the Banner was read correctly.
    broadcast(*pComm, rootRank, ptr(&bannerIsCorrect));
 
    // If the Banner is invalid, all processes throw an
    // exception.  Only Proc 0 gets the exception message, but
    // that's OK, since the main point is to "stop the world"
    // (rather than throw an exception on one process and leave
    // the others hanging).
    TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
                               std::invalid_argument, errMsg.str());
  }  // Done reading the Banner line and broadcasting success.
  if (debug && myRank == rootRank) {
    cerr << "-- Reading dimensions line" << endl;
  }
 
  // Read the matrix dimensions from the Matrix Market metadata.
  // dims = (numRows, numCols, numEntries).  Proc 0 does the
  // reading, but it broadcasts the results to all MPI
  // processes.  Thus, readCoordDims() is a collective
  // operation.  It does a collective check for correctness too.
  Tuple<global_ordinal_type, 3> dims =
      readCoordDims(in, lineNumber, pBanner, pComm, tolerant, debug);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Making Adder for collecting matrix data" << endl;
  }
 
  // "Adder" object for collecting all the sparse matrix entries
  // from the input stream.  This is only nonnull on Proc 0.
  RCP<adder_type> pAdder =
      makeAdder(pComm, pBanner, dims, tolerant, debug);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Reading matrix data" << endl;
  }
  //
  // Read the matrix entries from the input stream on Proc 0.
  //
  {
    // We use readSuccess to broadcast the results of the read
    // (succeeded or not) to all MPI processes.  Since
    // Teuchos::broadcast doesn't currently know how to send
    // bools, we convert to int (true -> 1, false -> 0).
    int readSuccess = 1;
    std::ostringstream errMsg;  // Exception message (only valid on Proc 0)
    if (myRank == rootRank) {
      try {
        // Reader for "coordinate" format sparse matrix data.
        typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
                                                       global_ordinal_type, scalar_type, STS::isComplex>
            reader_type;
        reader_type reader(pAdder);
 
        // Read the sparse matrix entries.
        std::pair<bool, std::vector<size_t>> results =
            reader.read(in, lineNumber, tolerant, debug);
        readSuccess = results.first ? 1 : 0;
      } catch (std::exception& e) {
        readSuccess = 0;
        errMsg << e.what();
      }
    }
    broadcast(*pComm, rootRank, ptr(&readSuccess));
 
    // It would be nice to add a "verbose" flag, so that in
    // tolerant mode, we could log any bad line number(s) on
    // Proc 0.  For now, we just throw if the read fails to
    // succeed.
    //
    // Question: If we're in tolerant mode, and if the read did
    // not succeed, should we attempt to call fillComplete()?
    TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
                               "Failed to read the Matrix Market sparse matrix file: "
                                   << errMsg.str());
  }  // Done reading the matrix entries (stored on Proc 0 for now)
 
  if (debug && myRank == rootRank) {
    cerr << "-- Successfully read the Matrix Market data" << endl;
  }
 
  // In tolerant mode, we need to rebroadcast the matrix
  // dimensions, since they may be different after reading the
  // actual matrix data.  We only need to broadcast the number
  // of rows and columns.  Only Rank 0 needs to know the actual
  // global number of entries, since (a) we need to merge
  // duplicates on Rank 0 first anyway, and (b) when we
  // distribute the entries, each rank other than Rank 0 will
  // only need to know how many entries it owns, not the total
  // number of entries.
  if (tolerant) {
    if (debug && myRank == rootRank) {
      cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
           << endl
           << "----- Dimensions before: "
           << dims[0] << " x " << dims[1]
           << endl;
    }
    // Packed coordinate matrix dimensions (numRows, numCols).
    Tuple<global_ordinal_type, 2> updatedDims;
    if (myRank == rootRank) {
      // If one or more bottom rows of the matrix contain no
      // entries, then the Adder will report that the number
      // of rows is less than that specified in the
      // metadata.  We allow this case, and favor the
      // metadata so that the zero row(s) will be included.
      updatedDims[0] =
          std::max(dims[0], pAdder->getAdder()->numRows());
      updatedDims[1] = pAdder->getAdder()->numCols();
    }
    broadcast(*pComm, rootRank, updatedDims);
    dims[0] = updatedDims[0];
    dims[1] = updatedDims[1];
    if (debug && myRank == rootRank) {
      cerr << "----- Dimensions after: " << dims[0] << " x "
           << dims[1] << endl;
    }
  } else {
    // In strict mode, we require that the matrix's metadata and
    // its actual data agree, at least somewhat.  In particular,
    // the number of rows must agree, since otherwise we cannot
    // distribute the matrix correctly.
 
    // Teuchos::broadcast() doesn't know how to broadcast bools,
    // so we use an int with the standard 1 == true, 0 == false
    // encoding.
    int dimsMatch = 1;
    if (myRank == rootRank) {
      // If one or more bottom rows of the matrix contain no
      // entries, then the Adder will report that the number of
      // rows is less than that specified in the metadata.  We
      // allow this case, and favor the metadata, but do not
      // allow the Adder to think there are more rows in the
      // matrix than the metadata says.
      if (dims[0] < pAdder->getAdder()->numRows()) {
        dimsMatch = 0;
      }
    }
    broadcast(*pComm, 0, ptr(&dimsMatch));
    if (dimsMatch == 0) {
      // We're in an error state anyway, so we might as well
      // work a little harder to print an informative error
      // message.
      //
      // Broadcast the Adder's idea of the matrix dimensions
      // from Proc 0 to all processes.
      Tuple<global_ordinal_type, 2> addersDims;
      if (myRank == rootRank) {
        addersDims[0] = pAdder->getAdder()->numRows();
        addersDims[1] = pAdder->getAdder()->numCols();
      }
      broadcast(*pComm, 0, addersDims);
      TEUCHOS_TEST_FOR_EXCEPTION(
          dimsMatch == 0, std::runtime_error,
          "The matrix metadata says that the matrix is " << dims[0] << " x "
                                                         << dims[1] << ", but the actual data says that the matrix is "
                                                         << addersDims[0] << " x " << addersDims[1] << ".  That means the "
                                                                                                       "data includes more rows than reported in the metadata.  This "
                                                                                                       "is not allowed when parsing in strict mode.  Parse the matrix in "
                                                                                                       "tolerant mode to ignore the metadata when it disagrees with the "
                                                                                                       "data.");
    }
  }  // Matrix dimensions (# rows, # cols, # entries) agree.
 
  if (debug && myRank == rootRank) {
    cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
  }
 
  // Now that we've read in all the matrix entries from the
  // input stream into the adder on Proc 0, post-process them
  // into CSR format (still on Proc 0).  This will facilitate
  // distributing them to all the processors.
  //
  // These arrays represent the global matrix data as a CSR
  // matrix (with numEntriesPerRow as redundant but convenient
  // metadata, since it's computable from rowPtr and vice
  // versa).  They are valid only on Proc 0.
  ArrayRCP<size_t> numEntriesPerRow;
  ArrayRCP<size_t> rowPtr;
  ArrayRCP<global_ordinal_type> colInd;
  ArrayRCP<scalar_type> values;
 
  // Proc 0 first merges duplicate entries, and then converts
  // the coordinate-format matrix data to CSR.
  {
    int mergeAndConvertSucceeded = 1;
    std::ostringstream errMsg;
 
    if (myRank == rootRank) {
      try {
        typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
                                                    global_ordinal_type>
            element_type;
 
        // Number of rows in the matrix.  If we are in tolerant
        // mode, we've already synchronized dims with the actual
        // matrix data.  If in strict mode, we should use dims
        // (as read from the file's metadata) rather than the
        // matrix data to determine the dimensions.  (The matrix
        // data will claim fewer rows than the metadata, if one
        // or more rows have no entries stored in the file.)
        const size_type numRows = dims[0];
 
        // Additively merge duplicate matrix entries.
        pAdder->getAdder()->merge();
 
        // Get a temporary const view of the merged matrix entries.
        const std::vector<element_type>& entries =
            pAdder->getAdder()->getEntries();
 
        // Number of matrix entries (after merging).
        const size_t numEntries = (size_t)entries.size();
 
        if (debug) {
          cerr << "----- Proc 0: Matrix has numRows=" << numRows
               << " rows and numEntries=" << numEntries
               << " entries." << endl;
        }
 
        // Make space for the CSR matrix data.  Converting to
        // CSR is easier if we fill numEntriesPerRow with zeros
        // at first.
        numEntriesPerRow = arcp<size_t>(numRows);
        std::fill(numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
        rowPtr = arcp<size_t>(numRows + 1);
        std::fill(rowPtr.begin(), rowPtr.end(), 0);
        colInd = arcp<global_ordinal_type>(numEntries);
        values = arcp<scalar_type>(numEntries);
 
        // Convert from array-of-structs coordinate format to CSR
        // (compressed sparse row) format.
        global_ordinal_type prvRow = 0;
        size_t curPos              = 0;
        rowPtr[0]                  = 0;
        for (curPos = 0; curPos < numEntries; ++curPos) {
          const element_type& curEntry     = entries[curPos];
          const global_ordinal_type curRow = curEntry.rowIndex();
          TEUCHOS_TEST_FOR_EXCEPTION(
              curRow < prvRow, std::logic_error,
              "Row indices are out of order, even though they are supposed "
              "to be sorted.  curRow = "
                  << curRow << ", prvRow = "
                  << prvRow << ", at curPos = " << curPos << ".  Please report "
                                                             "this bug to the Tpetra developers.");
          if (curRow > prvRow) {
            for (global_ordinal_type r = prvRow + 1; r <= curRow; ++r) {
              rowPtr[r] = curPos;
            }
            prvRow = curRow;
          }
          numEntriesPerRow[curRow]++;
          colInd[curPos] = curEntry.colIndex();
          values[curPos] = curEntry.value();
        }
        // rowPtr has one more entry than numEntriesPerRow.  The
        // last entry of rowPtr is the number of entries in
        // colInd and values.
        rowPtr[numRows] = numEntries;
      }  // Finished conversion to CSR format
      catch (std::exception& e) {
        mergeAndConvertSucceeded = 0;
        errMsg << "Failed to merge sparse matrix entries and convert to "
                  "CSR format: "
               << e.what();
      }
 
      if (debug && mergeAndConvertSucceeded) {
        // Number of rows in the matrix.
        const size_type numRows      = dims[0];
        const size_type maxToDisplay = 100;
 
        cerr << "----- Proc 0: numEntriesPerRow[0.."
             << (numEntriesPerRow.size() - 1) << "] ";
        if (numRows > maxToDisplay) {
          cerr << "(only showing first and last few entries) ";
        }
        cerr << "= [";
        if (numRows > 0) {
          if (numRows > maxToDisplay) {
            for (size_type k = 0; k < 2; ++k) {
              cerr << numEntriesPerRow[k] << " ";
            }
            cerr << "... ";
            for (size_type k = numRows - 2; k < numRows - 1; ++k) {
              cerr << numEntriesPerRow[k] << " ";
            }
          } else {
            for (size_type k = 0; k < numRows - 1; ++k) {
              cerr << numEntriesPerRow[k] << " ";
            }
          }
          cerr << numEntriesPerRow[numRows - 1];
        }  // numRows > 0
        cerr << "]" << endl;
 
        cerr << "----- Proc 0: rowPtr ";
        if (numRows > maxToDisplay) {
          cerr << "(only showing first and last few entries) ";
        }
        cerr << "= [";
        if (numRows > maxToDisplay) {
          for (size_type k = 0; k < 2; ++k) {
            cerr << rowPtr[k] << " ";
          }
          cerr << "... ";
          for (size_type k = numRows - 2; k < numRows; ++k) {
            cerr << rowPtr[k] << " ";
          }
        } else {
          for (size_type k = 0; k < numRows; ++k) {
            cerr << rowPtr[k] << " ";
          }
        }
        cerr << rowPtr[numRows] << "]" << endl;
      }
 
    }  // if myRank == rootRank
 
    // Broadcast success flag to all ranks
    int mergeAndConvertSucceededAll = mergeAndConvertSucceeded;
    broadcast(*pComm, rootRank, ptr(&mergeAndConvertSucceededAll));
    if (mergeAndConvertSucceededAll == 0) {
      TEUCHOS_TEST_FOR_EXCEPTION(
          true, std::runtime_error,
          "Failed to merge sparse matrix entries and convert to CSR format: "
              << errMsg.str());
    }
  }  // Done converting sparse matrix data to CSR format
 
  // Now we're done with the Adder, so we can release the
  // reference ("free" it) to save space.  This only actually
  // does anything on Rank 0, since pAdder is null on all the
  // other MPI processes.
  pAdder = null;
 
  if (debug && myRank == rootRank) {
    cerr << "-- Making range, domain, and row maps" << endl;
  }
 
  // Make the maps that describe the matrix's range and domain,
  // and the distribution of its rows.  Creating a Map is a
  // collective operation, so we don't have to do a broadcast of
  // a success Boolean.
  RCP<const map_type> pRangeMap = makeRangeMap(pComm, dims[0]);
  RCP<const map_type> pDomainMap =
      makeDomainMap(pRangeMap, dims[0], dims[1]);
  RCP<const map_type> pRowMap = makeRowMap(null, pComm, dims[0]);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Distributing the matrix data" << endl;
  }
 
  // Distribute the matrix data.  Each processor has to add the
  // rows that it owns.  If you try to make Proc 0 call
  // insertGlobalValues() for _all_ the rows, not just those it
  // owns, then fillComplete() will compute the number of
  // columns incorrectly.  That's why Proc 0 has to distribute
  // the matrix data and why we make all the processors (not
  // just Proc 0) call insertGlobalValues() on their own data.
  //
  // These arrays represent each processor's part of the matrix
  // data, in "CSR" format (sort of, since the row indices might
  // not be contiguous).
  ArrayRCP<size_t> myNumEntriesPerRow;
  ArrayRCP<size_t> myRowPtr;
  ArrayRCP<global_ordinal_type> myColInd;
  ArrayRCP<scalar_type> myValues;
  // Distribute the matrix data.  This is a collective operation.
  distribute(myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
             numEntriesPerRow, rowPtr, colInd, values, debug);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Inserting matrix entries on each processor";
    if (callFillComplete) {
      cerr << " and calling fillComplete()";
    }
    cerr << endl;
  }
  // Each processor inserts its part of the matrix data, and
  // then they all call fillComplete().  This method invalidates
  // the my* distributed matrix data before calling
  // fillComplete(), in order to save space.  In general, we
  // never store more than two copies of the matrix's entries in
  // memory at once, which is no worse than what Tpetra
  // promises.
  RCP<sparse_matrix_type> pMatrix =
      makeMatrix(myNumEntriesPerRow, myRowPtr, myColInd, myValues,
                 pRowMap, pRangeMap, pDomainMap, callFillComplete);
  // Only use a reduce-all in debug mode to check if pMatrix is
  // null.  Otherwise, just throw an exception.  We never expect
  // a null pointer here, so we can save a communication.
  if (debug) {
    int localIsNull  = pMatrix.is_null() ? 1 : 0;
    int globalIsNull = 0;
    reduceAll(*pComm, REDUCE_MAX, localIsNull, ptr(&globalIsNull));
    TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
                               "MatrixMarketReader::makeMatrix() returned a null pointer on at least one "
                               "process.  Please report this bug to the Tpetra developers.");
  } else {
    TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
                               "MatrixMarketReader::makeMatrix() returned a null pointer.  "
                               "Please report this bug to the Tpetra developers.");
  }
 
  // We can't get the dimensions of the matrix until after
  // fillComplete() is called.  Thus, we can't do the sanity
  // check (dimensions read from the Matrix Market data,
  // vs. dimensions reported by the CrsMatrix) unless the user
  // asked makeMatrix() to call fillComplete().
  //
  // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
  // what one might think it does, so you have to ask the range
  // resp. domain map for the number of rows resp. columns.
  if (callFillComplete) {
    const int numProcs = pComm->getSize();
 
    if (extraDebug && debug) {
      const global_size_t globalNumRows =
          pRangeMap->getGlobalNumElements();
      const global_size_t globalNumCols =
          pDomainMap->getGlobalNumElements();
      if (myRank == rootRank) {
        cerr << "-- Matrix is "
             << globalNumRows << " x " << globalNumCols
             << " with " << pMatrix->getGlobalNumEntries()
             << " entries, and index base "
             << pMatrix->getIndexBase() << "." << endl;
      }
      pComm->barrier();
      for (int p = 0; p < numProcs; ++p) {
        if (myRank == p) {
          cerr << "-- Proc " << p << " owns "
               << pMatrix->getLocalNumCols() << " columns, and "
               << pMatrix->getLocalNumEntries() << " entries." << endl;
        }
        pComm->barrier();
      }
    }  // if (extraDebug && debug)
  }    // if (callFillComplete)
 
  if (debug && myRank == rootRank) {
    cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
         << endl;
  }
  return pMatrix;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparse(std::istream& in,
                                                                          const Teuchos::RCP<const Teuchos::Comm<int>>& pComm,
                                                                          const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
                                                                          const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
                                                                          const bool tolerant,
                                                                          const bool debug) {
  using std::cerr;
  using std::endl;
  using Teuchos::arcp;
  using Teuchos::ArrayRCP;
  using Teuchos::broadcast;
  using Teuchos::null;
  using Teuchos::ptr;
  using Teuchos::RCP;
  using Teuchos::reduceAll;
  using Teuchos::Tuple;
  using Teuchos::MatrixMarket::Banner;
  typedef Teuchos::ScalarTraits<scalar_type> STS;
 
  const bool extraDebug = false;
  const int myRank      = pComm->getRank();
  const int rootRank    = 0;
 
  // Current line number in the input stream.  Various calls
  // will modify this depending on the number of lines that are
  // read from the input stream.  Only Rank 0 modifies this.
  size_t lineNumber = 1;
 
  if (debug && myRank == rootRank) {
    cerr << "Matrix Market reader: readSparse:" << endl
         << "-- Reading banner line" << endl;
  }
 
  // The "Banner" tells you whether the input stream represents
  // a sparse matrix, the symmetry type of the matrix, and the
  // type of the data it contains.
  //
  // pBanner will only be nonnull on MPI Rank 0.  It will be
  // null on all other MPI processes.
  RCP<const Banner> pBanner;
  {
    // We read and validate the Banner on Proc 0, but broadcast
    // the validation result to all processes.
    // Teuchos::broadcast doesn't currently work with bool, so
    // we use int (true -> 1, false -> 0).
    int bannerIsCorrect = 1;
    std::ostringstream errMsg;
 
    if (myRank == rootRank) {
      // Read the Banner line from the input stream.
      try {
        pBanner = readBanner(in, lineNumber, tolerant, debug);
      } catch (std::exception& e) {
        errMsg << "Attempt to read the Matrix Market file's Banner line "
                  "threw an exception: "
               << e.what();
        bannerIsCorrect = 0;
      }
 
      if (bannerIsCorrect) {
        // Validate the Banner for the case of a sparse matrix.
        // We validate on Proc 0, since it reads the Banner.
 
        // In intolerant mode, the matrix type must be "coordinate".
        if (!tolerant && pBanner->matrixType() != "coordinate") {
          bannerIsCorrect = 0;
          errMsg << "The Matrix Market input file must contain a "
                    "\"coordinate\"-format sparse matrix in order to create a "
                    "Tpetra::CrsMatrix object from it, but the file's matrix "
                    "type is \""
                 << pBanner->matrixType() << "\" instead.";
        }
        // In tolerant mode, we allow the matrix type to be
        // anything other than "array" (which would mean that
        // the file contains a dense matrix).
        if (tolerant && pBanner->matrixType() == "array") {
          bannerIsCorrect = 0;
          errMsg << "Matrix Market file must contain a \"coordinate\"-"
                    "format sparse matrix in order to create a Tpetra::CrsMatrix "
                    "object from it, but the file's matrix type is \"array\" "
                    "instead.  That probably means the file contains dense matrix "
                    "data.";
        }
      }
    }  // Proc 0: Done reading the Banner, hopefully successfully.
 
    // Broadcast from Proc 0 whether the Banner was read correctly.
    broadcast(*pComm, rootRank, ptr(&bannerIsCorrect));
 
    // If the Banner is invalid, all processes throw an
    // exception.  Only Proc 0 gets the exception message, but
    // that's OK, since the main point is to "stop the world"
    // (rather than throw an exception on one process and leave
    // the others hanging).
    TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
                               std::invalid_argument, errMsg.str());
  }  // Done reading the Banner line and broadcasting success.
  if (debug && myRank == rootRank) {
    cerr << "-- Reading dimensions line" << endl;
  }
 
  // Read the matrix dimensions from the Matrix Market metadata.
  // dims = (numRows, numCols, numEntries).  Proc 0 does the
  // reading, but it broadcasts the results to all MPI
  // processes.  Thus, readCoordDims() is a collective
  // operation.  It does a collective check for correctness too.
  Tuple<global_ordinal_type, 3> dims =
      readCoordDims(in, lineNumber, pBanner, pComm, tolerant, debug);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Making Adder for collecting matrix data" << endl;
  }
 
  // "Adder" object for collecting all the sparse matrix entries
  // from the input stream.  This is only nonnull on Proc 0.
  RCP<adder_type> pAdder =
      makeAdder(pComm, pBanner, dims, tolerant, debug);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Reading matrix data" << endl;
  }
  //
  // Read the matrix entries from the input stream on Proc 0.
  //
  {
    // We use readSuccess to broadcast the results of the read
    // (succeeded or not) to all MPI processes.  Since
    // Teuchos::broadcast doesn't currently know how to send
    // bools, we convert to int (true -> 1, false -> 0).
    int readSuccess = 1;
    std::ostringstream errMsg;  // Exception message (only valid on Proc 0)
    if (myRank == rootRank) {
      try {
        // Reader for "coordinate" format sparse matrix data.
        typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
                                                       global_ordinal_type, scalar_type, STS::isComplex>
            reader_type;
        reader_type reader(pAdder);
 
        // Read the sparse matrix entries.
        std::pair<bool, std::vector<size_t>> results =
            reader.read(in, lineNumber, tolerant, debug);
        readSuccess = results.first ? 1 : 0;
      } catch (std::exception& e) {
        readSuccess = 0;
        errMsg << e.what();
      }
    }
    broadcast(*pComm, rootRank, ptr(&readSuccess));
 
    // It would be nice to add a "verbose" flag, so that in
    // tolerant mode, we could log any bad line number(s) on
    // Proc 0.  For now, we just throw if the read fails to
    // succeed.
    //
    // Question: If we're in tolerant mode, and if the read did
    // not succeed, should we attempt to call fillComplete()?
    TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
                               "Failed to read the Matrix Market sparse matrix file: "
                                   << errMsg.str());
  }  // Done reading the matrix entries (stored on Proc 0 for now)
 
  if (debug && myRank == rootRank) {
    cerr << "-- Successfully read the Matrix Market data" << endl;
  }
 
  // In tolerant mode, we need to rebroadcast the matrix
  // dimensions, since they may be different after reading the
  // actual matrix data.  We only need to broadcast the number
  // of rows and columns.  Only Rank 0 needs to know the actual
  // global number of entries, since (a) we need to merge
  // duplicates on Rank 0 first anyway, and (b) when we
  // distribute the entries, each rank other than Rank 0 will
  // only need to know how many entries it owns, not the total
  // number of entries.
  if (tolerant) {
    if (debug && myRank == rootRank) {
      cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
           << endl
           << "----- Dimensions before: "
           << dims[0] << " x " << dims[1]
           << endl;
    }
    // Packed coordinate matrix dimensions (numRows, numCols).
    Tuple<global_ordinal_type, 2> updatedDims;
    if (myRank == rootRank) {
      // If one or more bottom rows of the matrix contain no
      // entries, then the Adder will report that the number
      // of rows is less than that specified in the
      // metadata.  We allow this case, and favor the
      // metadata so that the zero row(s) will be included.
      updatedDims[0] =
          std::max(dims[0], pAdder->getAdder()->numRows());
      updatedDims[1] = pAdder->getAdder()->numCols();
    }
    broadcast(*pComm, rootRank, updatedDims);
    dims[0] = updatedDims[0];
    dims[1] = updatedDims[1];
    if (debug && myRank == rootRank) {
      cerr << "----- Dimensions after: " << dims[0] << " x "
           << dims[1] << endl;
    }
  } else {
    // In strict mode, we require that the matrix's metadata and
    // its actual data agree, at least somewhat.  In particular,
    // the number of rows must agree, since otherwise we cannot
    // distribute the matrix correctly.
 
    // Teuchos::broadcast() doesn't know how to broadcast bools,
    // so we use an int with the standard 1 == true, 0 == false
    // encoding.
    int dimsMatch = 1;
    if (myRank == rootRank) {
      // If one or more bottom rows of the matrix contain no
      // entries, then the Adder will report that the number of
      // rows is less than that specified in the metadata.  We
      // allow this case, and favor the metadata, but do not
      // allow the Adder to think there are more rows in the
      // matrix than the metadata says.
      if (dims[0] < pAdder->getAdder()->numRows()) {
        dimsMatch = 0;
      }
    }
    broadcast(*pComm, 0, ptr(&dimsMatch));
    if (dimsMatch == 0) {
      // We're in an error state anyway, so we might as well
      // work a little harder to print an informative error
      // message.
      //
      // Broadcast the Adder's idea of the matrix dimensions
      // from Proc 0 to all processes.
      Tuple<global_ordinal_type, 2> addersDims;
      if (myRank == rootRank) {
        addersDims[0] = pAdder->getAdder()->numRows();
        addersDims[1] = pAdder->getAdder()->numCols();
      }
      broadcast(*pComm, 0, addersDims);
      TEUCHOS_TEST_FOR_EXCEPTION(
          dimsMatch == 0, std::runtime_error,
          "The matrix metadata says that the matrix is " << dims[0] << " x "
                                                         << dims[1] << ", but the actual data says that the matrix is "
                                                         << addersDims[0] << " x " << addersDims[1] << ".  That means the "
                                                                                                       "data includes more rows than reported in the metadata.  This "
                                                                                                       "is not allowed when parsing in strict mode.  Parse the matrix in "
                                                                                                       "tolerant mode to ignore the metadata when it disagrees with the "
                                                                                                       "data.");
    }
  }  // Matrix dimensions (# rows, # cols, # entries) agree.
 
  if (debug && myRank == rootRank) {
    cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
  }
 
  // Now that we've read in all the matrix entries from the
  // input stream into the adder on Proc 0, post-process them
  // into CSR format (still on Proc 0).  This will facilitate
  // distributing them to all the processors.
  //
  // These arrays represent the global matrix data as a CSR
  // matrix (with numEntriesPerRow as redundant but convenient
  // metadata, since it's computable from rowPtr and vice
  // versa).  They are valid only on Proc 0.
  ArrayRCP<size_t> numEntriesPerRow;
  ArrayRCP<size_t> rowPtr;
  ArrayRCP<global_ordinal_type> colInd;
  ArrayRCP<scalar_type> values;
 
  // Proc 0 first merges duplicate entries, and then converts
  // the coordinate-format matrix data to CSR.
  {
    int mergeAndConvertSucceeded = 1;
    std::ostringstream errMsg;
 
    if (myRank == rootRank) {
      try {
        typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
                                                    global_ordinal_type>
            element_type;
 
        // Number of rows in the matrix.  If we are in tolerant
        // mode, we've already synchronized dims with the actual
        // matrix data.  If in strict mode, we should use dims
        // (as read from the file's metadata) rather than the
        // matrix data to determine the dimensions.  (The matrix
        // data will claim fewer rows than the metadata, if one
        // or more rows have no entries stored in the file.)
        const size_type numRows = dims[0];
 
        // Additively merge duplicate matrix entries.
        pAdder->getAdder()->merge();
 
        // Get a temporary const view of the merged matrix entries.
        const std::vector<element_type>& entries =
            pAdder->getAdder()->getEntries();
 
        // Number of matrix entries (after merging).
        const size_t numEntries = (size_t)entries.size();
 
        if (debug) {
          cerr << "----- Proc 0: Matrix has numRows=" << numRows
               << " rows and numEntries=" << numEntries
               << " entries." << endl;
        }
 
        // Make space for the CSR matrix data.  Converting to
        // CSR is easier if we fill numEntriesPerRow with zeros
        // at first.
        numEntriesPerRow = arcp<size_t>(numRows);
        std::fill(numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
        rowPtr = arcp<size_t>(numRows + 1);
        std::fill(rowPtr.begin(), rowPtr.end(), 0);
        colInd = arcp<global_ordinal_type>(numEntries);
        values = arcp<scalar_type>(numEntries);
 
        // Convert from array-of-structs coordinate format to CSR
        // (compressed sparse row) format.
        global_ordinal_type prvRow = 0;
        size_t curPos              = 0;
        rowPtr[0]                  = 0;
        for (curPos = 0; curPos < numEntries; ++curPos) {
          const element_type& curEntry     = entries[curPos];
          const global_ordinal_type curRow = curEntry.rowIndex();
          TEUCHOS_TEST_FOR_EXCEPTION(
              curRow < prvRow, std::logic_error,
              "Row indices are out of order, even though they are supposed "
              "to be sorted.  curRow = "
                  << curRow << ", prvRow = "
                  << prvRow << ", at curPos = " << curPos << ".  Please report "
                                                             "this bug to the Tpetra developers.");
          if (curRow > prvRow) {
            for (global_ordinal_type r = prvRow + 1; r <= curRow; ++r) {
              rowPtr[r] = curPos;
            }
            prvRow = curRow;
          }
          numEntriesPerRow[curRow]++;
          colInd[curPos] = curEntry.colIndex();
          values[curPos] = curEntry.value();
        }
        // rowPtr has one more entry than numEntriesPerRow.  The
        // last entry of rowPtr is the number of entries in
        // colInd and values.
        rowPtr[numRows] = numEntries;
      }  // Finished conversion to CSR format
      catch (std::exception& e) {
        mergeAndConvertSucceeded = 0;
        errMsg << "Failed to merge sparse matrix entries and convert to "
                  "CSR format: "
               << e.what();
      }
 
      if (debug && mergeAndConvertSucceeded) {
        // Number of rows in the matrix.
        const size_type numRows      = dims[0];
        const size_type maxToDisplay = 100;
 
        cerr << "----- Proc 0: numEntriesPerRow[0.."
             << (numEntriesPerRow.size() - 1) << "] ";
        if (numRows > maxToDisplay) {
          cerr << "(only showing first and last few entries) ";
        }
        cerr << "= [";
        if (numRows > 0) {
          if (numRows > maxToDisplay) {
            for (size_type k = 0; k < 2; ++k) {
              cerr << numEntriesPerRow[k] << " ";
            }
            cerr << "... ";
            for (size_type k = numRows - 2; k < numRows - 1; ++k) {
              cerr << numEntriesPerRow[k] << " ";
            }
          } else {
            for (size_type k = 0; k < numRows - 1; ++k) {
              cerr << numEntriesPerRow[k] << " ";
            }
          }
          cerr << numEntriesPerRow[numRows - 1];
        }  // numRows > 0
        cerr << "]" << endl;
 
        cerr << "----- Proc 0: rowPtr ";
        if (numRows > maxToDisplay) {
          cerr << "(only showing first and last few entries) ";
        }
        cerr << "= [";
        if (numRows > maxToDisplay) {
          for (size_type k = 0; k < 2; ++k) {
            cerr << rowPtr[k] << " ";
          }
          cerr << "... ";
          for (size_type k = numRows - 2; k < numRows; ++k) {
            cerr << rowPtr[k] << " ";
          }
        } else {
          for (size_type k = 0; k < numRows; ++k) {
            cerr << rowPtr[k] << " ";
          }
        }
        cerr << rowPtr[numRows] << "]" << endl;
      }
    }  // if myRank == rootRank
 
    // Broadcast success flag to all ranks
    int mergeAndConvertSucceededAll = mergeAndConvertSucceeded;
    broadcast(*pComm, rootRank, ptr(&mergeAndConvertSucceededAll));
    if (mergeAndConvertSucceededAll == 0) {
      TEUCHOS_TEST_FOR_EXCEPTION(
          true, std::runtime_error,
          "Failed to merge sparse matrix entries and convert to CSR format: "
              << errMsg.str());
    }
  }  // Done converting sparse matrix data to CSR format
 
  // Now we're done with the Adder, so we can release the
  // reference ("free" it) to save space.  This only actually
  // does anything on Rank 0, since pAdder is null on all the
  // other MPI processes.
  pAdder = null;
 
  if (debug && myRank == rootRank) {
    cerr << "-- Making range, domain, and row maps" << endl;
  }
 
  // Make the maps that describe the matrix's range and domain,
  // and the distribution of its rows.  Creating a Map is a
  // collective operation, so we don't have to do a broadcast of
  // a success Boolean.
  RCP<const map_type> pRangeMap = makeRangeMap(pComm, dims[0]);
  RCP<const map_type> pDomainMap =
      makeDomainMap(pRangeMap, dims[0], dims[1]);
  RCP<const map_type> pRowMap = makeRowMap(null, pComm, dims[0]);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Distributing the matrix data" << endl;
  }
 
  // Distribute the matrix data.  Each processor has to add the
  // rows that it owns.  If you try to make Proc 0 call
  // insertGlobalValues() for _all_ the rows, not just those it
  // owns, then fillComplete() will compute the number of
  // columns incorrectly.  That's why Proc 0 has to distribute
  // the matrix data and why we make all the processors (not
  // just Proc 0) call insertGlobalValues() on their own data.
  //
  // These arrays represent each processor's part of the matrix
  // data, in "CSR" format (sort of, since the row indices might
  // not be contiguous).
  ArrayRCP<size_t> myNumEntriesPerRow;
  ArrayRCP<size_t> myRowPtr;
  ArrayRCP<global_ordinal_type> myColInd;
  ArrayRCP<scalar_type> myValues;
  // Distribute the matrix data.  This is a collective operation.
  distribute(myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
             numEntriesPerRow, rowPtr, colInd, values, debug);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Inserting matrix entries on each process "
            "and calling fillComplete()"
         << endl;
  }
  // Each processor inserts its part of the matrix data, and
  // then they all call fillComplete().  This method invalidates
  // the my* distributed matrix data before calling
  // fillComplete(), in order to save space.  In general, we
  // never store more than two copies of the matrix's entries in
  // memory at once, which is no worse than what Tpetra
  // promises.
  Teuchos::RCP<sparse_matrix_type> pMatrix =
      makeMatrix(myNumEntriesPerRow, myRowPtr, myColInd, myValues,
                 pRowMap, pRangeMap, pDomainMap, constructorParams,
                 fillCompleteParams);
  // Only use a reduce-all in debug mode to check if pMatrix is
  // null.  Otherwise, just throw an exception.  We never expect
  // a null pointer here, so we can save a communication.
  if (debug) {
    int localIsNull  = pMatrix.is_null() ? 1 : 0;
    int globalIsNull = 0;
    reduceAll(*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr(&globalIsNull));
    TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
                               "MatrixMarketReader::makeMatrix() returned a null pointer on at least one "
                               "process.  Please report this bug to the Tpetra developers.");
  } else {
    TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
                               "MatrixMarketReader::makeMatrix() returned a null pointer.  "
                               "Please report this bug to the Tpetra developers.");
  }
 
  // Sanity check for dimensions (read from the Matrix Market
  // data, vs. dimensions reported by the CrsMatrix).
  //
  // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
  // what one might think it does, so you have to ask the range
  // resp. domain map for the number of rows resp. columns.
  if (extraDebug && debug) {
    const int numProcs = pComm->getSize();
    const global_size_t globalNumRows =
        pRangeMap->getGlobalNumElements();
    const global_size_t globalNumCols =
        pDomainMap->getGlobalNumElements();
    if (myRank == rootRank) {
      cerr << "-- Matrix is "
           << globalNumRows << " x " << globalNumCols
           << " with " << pMatrix->getGlobalNumEntries()
           << " entries, and index base "
           << pMatrix->getIndexBase() << "." << endl;
    }
    pComm->barrier();
    for (int p = 0; p < numProcs; ++p) {
      if (myRank == p) {
        cerr << "-- Proc " << p << " owns "
             << pMatrix->getLocalNumCols() << " columns, and "
             << pMatrix->getLocalNumEntries() << " entries." << endl;
      }
      pComm->barrier();
    }
  }  // if (extraDebug && debug)
 
  if (debug && myRank == rootRank) {
    cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
         << endl;
  }
  return pMatrix;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparse(std::istream& in,
                                                                          const Teuchos::RCP<const map_type>& rowMap,
                                                                          Teuchos::RCP<const map_type>& colMap,
                                                                          const Teuchos::RCP<const map_type>& domainMap,
                                                                          const Teuchos::RCP<const map_type>& rangeMap,
                                                                          const bool callFillComplete,
                                                                          const bool tolerant,
                                                                          const bool debug) {
  using std::cerr;
  using std::endl;
  using Teuchos::arcp;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::as;
  using Teuchos::broadcast;
  using Teuchos::Comm;
  using Teuchos::null;
  using Teuchos::ptr;
  using Teuchos::RCP;
  using Teuchos::reduceAll;
  using Teuchos::Tuple;
  using Teuchos::MatrixMarket::Banner;
  typedef Teuchos::ScalarTraits<scalar_type> STS;
 
  RCP<const Comm<int>> pComm = rowMap->getComm();
  const int myRank           = pComm->getRank();
  const int rootRank         = 0;
  const bool extraDebug      = false;
 
  // Fast checks for invalid input.  We can't check other
  // attributes of the Maps until we've read in the matrix
  // dimensions.
  TEUCHOS_TEST_FOR_EXCEPTION(
      rowMap.is_null(), std::invalid_argument,
      "Row Map must be nonnull.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      rangeMap.is_null(), std::invalid_argument,
      "Range Map must be nonnull.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      domainMap.is_null(), std::invalid_argument,
      "Domain Map must be nonnull.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      rowMap->getComm().getRawPtr() != pComm.getRawPtr(),
      std::invalid_argument,
      "The specified row Map's communicator (rowMap->getComm())"
      "differs from the given separately supplied communicator pComm.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      domainMap->getComm().getRawPtr() != pComm.getRawPtr(),
      std::invalid_argument,
      "The specified domain Map's communicator (domainMap->getComm())"
      "differs from the given separately supplied communicator pComm.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      rangeMap->getComm().getRawPtr() != pComm.getRawPtr(),
      std::invalid_argument,
      "The specified range Map's communicator (rangeMap->getComm())"
      "differs from the given separately supplied communicator pComm.");
 
  // Current line number in the input stream.  Various calls
  // will modify this depending on the number of lines that are
  // read from the input stream.  Only Rank 0 modifies this.
  size_t lineNumber = 1;
 
  if (debug && myRank == rootRank) {
    cerr << "Matrix Market reader: readSparse:" << endl
         << "-- Reading banner line" << endl;
  }
 
  // The "Banner" tells you whether the input stream represents
  // a sparse matrix, the symmetry type of the matrix, and the
  // type of the data it contains.
  //
  // pBanner will only be nonnull on MPI Rank 0.  It will be
  // null on all other MPI processes.
  RCP<const Banner> pBanner;
  {
    // We read and validate the Banner on Proc 0, but broadcast
    // the validation result to all processes.
    // Teuchos::broadcast doesn't currently work with bool, so
    // we use int (true -> 1, false -> 0).
    int bannerIsCorrect = 1;
    std::ostringstream errMsg;
 
    if (myRank == rootRank) {
      // Read the Banner line from the input stream.
      try {
        pBanner = readBanner(in, lineNumber, tolerant, debug);
      } catch (std::exception& e) {
        errMsg << "Attempt to read the Matrix Market file's Banner line "
                  "threw an exception: "
               << e.what();
        bannerIsCorrect = 0;
      }
 
      if (bannerIsCorrect) {
        // Validate the Banner for the case of a sparse matrix.
        // We validate on Proc 0, since it reads the Banner.
 
        // In intolerant mode, the matrix type must be "coordinate".
        if (!tolerant && pBanner->matrixType() != "coordinate") {
          bannerIsCorrect = 0;
          errMsg << "The Matrix Market input file must contain a "
                    "\"coordinate\"-format sparse matrix in order to create a "
                    "Tpetra::CrsMatrix object from it, but the file's matrix "
                    "type is \""
                 << pBanner->matrixType() << "\" instead.";
        }
        // In tolerant mode, we allow the matrix type to be
        // anything other than "array" (which would mean that
        // the file contains a dense matrix).
        if (tolerant && pBanner->matrixType() == "array") {
          bannerIsCorrect = 0;
          errMsg << "Matrix Market file must contain a \"coordinate\"-"
                    "format sparse matrix in order to create a Tpetra::CrsMatrix "
                    "object from it, but the file's matrix type is \"array\" "
                    "instead.  That probably means the file contains dense matrix "
                    "data.";
        }
      }
    }  // Proc 0: Done reading the Banner, hopefully successfully.
 
    // Broadcast from Proc 0 whether the Banner was read correctly.
    broadcast(*pComm, rootRank, ptr(&bannerIsCorrect));
 
    // If the Banner is invalid, all processes throw an
    // exception.  Only Proc 0 gets the exception message, but
    // that's OK, since the main point is to "stop the world"
    // (rather than throw an exception on one process and leave
    // the others hanging).
    TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
                               std::invalid_argument, errMsg.str());
  }  // Done reading the Banner line and broadcasting success.
  if (debug && myRank == rootRank) {
    cerr << "-- Reading dimensions line" << endl;
  }
 
  // Read the matrix dimensions from the Matrix Market metadata.
  // dims = (numRows, numCols, numEntries).  Proc 0 does the
  // reading, but it broadcasts the results to all MPI
  // processes.  Thus, readCoordDims() is a collective
  // operation.  It does a collective check for correctness too.
  Tuple<global_ordinal_type, 3> dims =
      readCoordDims(in, lineNumber, pBanner, pComm, tolerant, debug);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Making Adder for collecting matrix data" << endl;
  }
 
  // "Adder" object for collecting all the sparse matrix entries
  // from the input stream.  This is only nonnull on Proc 0.
  // The Adder internally converts the one-based indices (native
  // Matrix Market format) into zero-based indices.
  RCP<adder_type> pAdder =
      makeAdder(pComm, pBanner, dims, tolerant, debug);
 
  if (debug && myRank == rootRank) {
    cerr << "-- Reading matrix data" << endl;
  }
  //
  // Read the matrix entries from the input stream on Proc 0.
  //
  {
    // We use readSuccess to broadcast the results of the read
    // (succeeded or not) to all MPI processes.  Since
    // Teuchos::broadcast doesn't currently know how to send
    // bools, we convert to int (true -> 1, false -> 0).
    int readSuccess = 1;
    std::ostringstream errMsg;  // Exception message (only valid on Proc 0)
    if (myRank == rootRank) {
      try {
        // Reader for "coordinate" format sparse matrix data.
        typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
                                                       global_ordinal_type, scalar_type, STS::isComplex>
            reader_type;
        reader_type reader(pAdder);
 
        // Read the sparse matrix entries.
        std::pair<bool, std::vector<size_t>> results =
            reader.read(in, lineNumber, tolerant, debug);
        readSuccess = results.first ? 1 : 0;
      } catch (std::exception& e) {
        readSuccess = 0;
        errMsg << e.what();
      }
    }
    broadcast(*pComm, rootRank, ptr(&readSuccess));
 
    // It would be nice to add a "verbose" flag, so that in
    // tolerant mode, we could log any bad line number(s) on
    // Proc 0.  For now, we just throw if the read fails to
    // succeed.
    //
    // Question: If we're in tolerant mode, and if the read did
    // not succeed, should we attempt to call fillComplete()?
    TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
                               "Failed to read the Matrix Market sparse matrix file: "
                                   << errMsg.str());
  }  // Done reading the matrix entries (stored on Proc 0 for now)
 
  if (debug && myRank == rootRank) {
    cerr << "-- Successfully read the Matrix Market data" << endl;
  }
 
  // In tolerant mode, we need to rebroadcast the matrix
  // dimensions, since they may be different after reading the
  // actual matrix data.  We only need to broadcast the number
  // of rows and columns.  Only Rank 0 needs to know the actual
  // global number of entries, since (a) we need to merge
  // duplicates on Rank 0 first anyway, and (b) when we
  // distribute the entries, each rank other than Rank 0 will
  // only need to know how many entries it owns, not the total
  // number of entries.
  if (tolerant) {
    if (debug && myRank == rootRank) {
      cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
           << endl
           << "----- Dimensions before: "
           << dims[0] << " x " << dims[1]
           << endl;
    }
    // Packed coordinate matrix dimensions (numRows, numCols).
    Tuple<global_ordinal_type, 2> updatedDims;
    if (myRank == rootRank) {
      // If one or more bottom rows of the matrix contain no
      // entries, then the Adder will report that the number
      // of rows is less than that specified in the
      // metadata.  We allow this case, and favor the
      // metadata so that the zero row(s) will be included.
      updatedDims[0] =
          std::max(dims[0], pAdder->getAdder()->numRows());
      updatedDims[1] = pAdder->getAdder()->numCols();
    }
    broadcast(*pComm, rootRank, updatedDims);
    dims[0] = updatedDims[0];
    dims[1] = updatedDims[1];
    if (debug && myRank == rootRank) {
      cerr << "----- Dimensions after: " << dims[0] << " x "
           << dims[1] << endl;
    }
  } else {
    // In strict mode, we require that the matrix's metadata and
    // its actual data agree, at least somewhat.  In particular,
    // the number of rows must agree, since otherwise we cannot
    // distribute the matrix correctly.
 
    // Teuchos::broadcast() doesn't know how to broadcast bools,
    // so we use an int with the standard 1 == true, 0 == false
    // encoding.
    int dimsMatch = 1;
    if (myRank == rootRank) {
      // If one or more bottom rows of the matrix contain no
      // entries, then the Adder will report that the number of
      // rows is less than that specified in the metadata.  We
      // allow this case, and favor the metadata, but do not
      // allow the Adder to think there are more rows in the
      // matrix than the metadata says.
      if (dims[0] < pAdder->getAdder()->numRows()) {
        dimsMatch = 0;
      }
    }
    broadcast(*pComm, 0, ptr(&dimsMatch));
    if (dimsMatch == 0) {
      // We're in an error state anyway, so we might as well
      // work a little harder to print an informative error
      // message.
      //
      // Broadcast the Adder's idea of the matrix dimensions
      // from Proc 0 to all processes.
      Tuple<global_ordinal_type, 2> addersDims;
      if (myRank == rootRank) {
        addersDims[0] = pAdder->getAdder()->numRows();
        addersDims[1] = pAdder->getAdder()->numCols();
      }
      broadcast(*pComm, 0, addersDims);
      TEUCHOS_TEST_FOR_EXCEPTION(
          dimsMatch == 0, std::runtime_error,
          "The matrix metadata says that the matrix is " << dims[0] << " x "
                                                         << dims[1] << ", but the actual data says that the matrix is "
                                                         << addersDims[0] << " x " << addersDims[1] << ".  That means the "
                                                                                                       "data includes more rows than reported in the metadata.  This "
                                                                                                       "is not allowed when parsing in strict mode.  Parse the matrix in "
                                                                                                       "tolerant mode to ignore the metadata when it disagrees with the "
                                                                                                       "data.");
    }
  }  // Matrix dimensions (# rows, # cols, # entries) agree.
 
  if (debug && myRank == rootRank) {
    cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
  }
 
  // Now that we've read in all the matrix entries from the
  // input stream into the adder on Proc 0, post-process them
  // into CSR format (still on Proc 0).  This will facilitate
  // distributing them to all the processors.
  //
  // These arrays represent the global matrix data as a CSR
  // matrix (with numEntriesPerRow as redundant but convenient
  // metadata, since it's computable from rowPtr and vice
  // versa).  They are valid only on Proc 0.
  ArrayRCP<size_t> numEntriesPerRow;
  ArrayRCP<size_t> rowPtr;
  ArrayRCP<global_ordinal_type> colInd;
  ArrayRCP<scalar_type> values;
 
  // Proc 0 first merges duplicate entries, and then converts
  // the coordinate-format matrix data to CSR.
  {
    int mergeAndConvertSucceeded = 1;
    std::ostringstream errMsg;
 
    if (myRank == rootRank) {
      try {
        typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
                                                    global_ordinal_type>
            element_type;
 
        // Number of rows in the matrix.  If we are in tolerant
        // mode, we've already synchronized dims with the actual
        // matrix data.  If in strict mode, we should use dims
        // (as read from the file's metadata) rather than the
        // matrix data to determine the dimensions.  (The matrix
        // data will claim fewer rows than the metadata, if one
        // or more rows have no entries stored in the file.)
        const size_type numRows = dims[0];
 
        // Additively merge duplicate matrix entries.
        pAdder->getAdder()->merge();
 
        // Get a temporary const view of the merged matrix entries.
        const std::vector<element_type>& entries =
            pAdder->getAdder()->getEntries();
 
        // Number of matrix entries (after merging).
        const size_t numEntries = (size_t)entries.size();
 
        if (debug) {
          cerr << "----- Proc 0: Matrix has numRows=" << numRows
               << " rows and numEntries=" << numEntries
               << " entries." << endl;
        }
 
        // Make space for the CSR matrix data.  Converting to
        // CSR is easier if we fill numEntriesPerRow with zeros
        // at first.
        numEntriesPerRow = arcp<size_t>(numRows);
        std::fill(numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
        rowPtr = arcp<size_t>(numRows + 1);
        std::fill(rowPtr.begin(), rowPtr.end(), 0);
        colInd = arcp<global_ordinal_type>(numEntries);
        values = arcp<scalar_type>(numEntries);
 
        // Convert from array-of-structs coordinate format to CSR
        // (compressed sparse row) format.
        global_ordinal_type prvRow = 0;
        size_t curPos              = 0;
        rowPtr[0]                  = 0;
        for (curPos = 0; curPos < numEntries; ++curPos) {
          const element_type& curEntry     = entries[curPos];
          const global_ordinal_type curRow = curEntry.rowIndex();
          TEUCHOS_TEST_FOR_EXCEPTION(
              curRow < prvRow, std::logic_error,
              "Row indices are out of order, even though they are supposed "
              "to be sorted.  curRow = "
                  << curRow << ", prvRow = "
                  << prvRow << ", at curPos = " << curPos << ".  Please report "
                                                             "this bug to the Tpetra developers.");
          if (curRow > prvRow) {
            prvRow = curRow;
          }
          numEntriesPerRow[curRow]++;
          colInd[curPos] = curEntry.colIndex();
          values[curPos] = curEntry.value();
        }
 
        rowPtr[0] = 0;
        for (global_ordinal_type row = 1; row <= numRows; ++row) {
          rowPtr[row] = numEntriesPerRow[row - 1] + rowPtr[row - 1];
        }
      }  // Finished conversion to CSR format
      catch (std::exception& e) {
        mergeAndConvertSucceeded = 0;
        errMsg << "Failed to merge sparse matrix entries and convert to "
                  "CSR format: "
               << e.what();
      }
 
      if (debug && mergeAndConvertSucceeded) {
        // Number of rows in the matrix.
        const size_type numRows      = dims[0];
        const size_type maxToDisplay = 100;
 
        cerr << "----- Proc 0: numEntriesPerRow[0.."
             << (numEntriesPerRow.size() - 1) << "] ";
        if (numRows > maxToDisplay) {
          cerr << "(only showing first and last few entries) ";
        }
        cerr << "= [";
        if (numRows > 0) {
          if (numRows > maxToDisplay) {
            for (size_type k = 0; k < 2; ++k) {
              cerr << numEntriesPerRow[k] << " ";
            }
            cerr << "... ";
            for (size_type k = numRows - 2; k < numRows - 1; ++k) {
              cerr << numEntriesPerRow[k] << " ";
            }
          } else {
            for (size_type k = 0; k < numRows - 1; ++k) {
              cerr << numEntriesPerRow[k] << " ";
            }
          }
          cerr << numEntriesPerRow[numRows - 1];
        }  // numRows > 0
        cerr << "]" << endl;
 
        cerr << "----- Proc 0: rowPtr ";
        if (numRows > maxToDisplay) {
          cerr << "(only showing first and last few entries) ";
        }
        cerr << "= [";
        if (numRows > maxToDisplay) {
          for (size_type k = 0; k < 2; ++k) {
            cerr << rowPtr[k] << " ";
          }
          cerr << "... ";
          for (size_type k = numRows - 2; k < numRows; ++k) {
            cerr << rowPtr[k] << " ";
          }
        } else {
          for (size_type k = 0; k < numRows; ++k) {
            cerr << rowPtr[k] << " ";
          }
        }
        cerr << rowPtr[numRows] << "]" << endl;
 
        cerr << "----- Proc 0: colInd = [";
        for (size_t k = 0; k < rowPtr[numRows]; ++k) {
          cerr << colInd[k] << " ";
        }
        cerr << "]" << endl;
      }
    }  // if myRank == rootRank
  }    // Done converting sparse matrix data to CSR format
 
  // Now we're done with the Adder, so we can release the
  // reference ("free" it) to save space.  This only actually
  // does anything on Rank 0, since pAdder is null on all the
  // other MPI processes.
  pAdder = null;
 
  // Verify details of the Maps.  Don't count the global number
  // of entries in the row Map, since that number doesn't
  // correctly count overlap.
  if (debug && myRank == rootRank) {
    cerr << "-- Verifying Maps" << endl;
  }
  TEUCHOS_TEST_FOR_EXCEPTION(
      as<global_ordinal_type>(dims[0]) != rangeMap->getMaxAllGlobalIndex() + 1 - rangeMap->getIndexBase(),
      std::invalid_argument,
      "The range Map has " << rangeMap->getMaxAllGlobalIndex()
                           << " max entry, but the matrix has a global number of rows " << dims[0]
                           << ".");
  TEUCHOS_TEST_FOR_EXCEPTION(
      as<global_ordinal_type>(dims[1]) != domainMap->getMaxAllGlobalIndex() + 1 - domainMap->getIndexBase(),
      std::invalid_argument,
      "The domain Map has " << domainMap->getMaxAllGlobalIndex()
                            << " max entry, but the matrix has a global number of columns "
                            << dims[1] << ".");
 
  // Create a row Map which is entirely owned on Proc 0.
  RCP<Teuchos::FancyOStream> err = debug ? Teuchos::getFancyOStream(Teuchos::rcpFromRef(cerr)) : null;
 
  RCP<const map_type> gatherRowMap = Details::computeGatherMap(rowMap, err, debug);
  ArrayView<const global_ordinal_type> myRows =
      gatherRowMap->getLocalElementList();
  const size_type myNumRows           = myRows.size();
  const global_ordinal_type indexBase = gatherRowMap->getIndexBase();
 
  ArrayRCP<size_t> gatherNumEntriesPerRow = arcp<size_t>(myNumRows);
  for (size_type i_ = 0; i_ < myNumRows; i_++) {
    gatherNumEntriesPerRow[i_] = numEntriesPerRow[myRows[i_] - indexBase];
  }
 
  // Create a matrix using this Map, and fill in on Proc 0.  We
  // know how many entries there are in each row, so we can use
  // static profile.
  RCP<sparse_matrix_type> A_proc0 =
      rcp(new sparse_matrix_type(gatherRowMap, gatherNumEntriesPerRow()));
  if (myRank == rootRank) {
    if (debug) {
      cerr << "-- Proc 0: Filling gather matrix" << endl;
    }
    if (debug) {
      cerr << "---- Rows: " << Teuchos::toString(myRows) << endl;
    }
 
    // Add Proc 0's matrix entries to the CrsMatrix.
    for (size_type i_ = 0; i_ < myNumRows; ++i_) {
      size_type i = myRows[i_] - indexBase;
 
      const size_type curPos                 = as<size_type>(rowPtr[i]);
      const local_ordinal_type curNumEntries = numEntriesPerRow[i];
      ArrayView<global_ordinal_type> curColInd =
          colInd.view(curPos, curNumEntries);
      ArrayView<scalar_type> curValues =
          values.view(curPos, curNumEntries);
 
      // Modify the column indices in place to have the right index base.
      for (size_type k = 0; k < curNumEntries; ++k) {
        curColInd[k] += indexBase;
      }
      if (debug) {
        cerr << "------ Columns: " << Teuchos::toString(curColInd) << endl;
        cerr << "------ Values: " << Teuchos::toString(curValues) << endl;
      }
      // Avoid constructing empty views of ArrayRCP objects.
      if (curNumEntries > 0) {
        A_proc0->insertGlobalValues(myRows[i_], curColInd, curValues);
      }
    }
    // Now we can save space by deallocating numEntriesPerRow,
    // rowPtr, colInd, and values, since we've already put those
    // data in the matrix.
    numEntriesPerRow = null;
    rowPtr           = null;
    colInd           = null;
    values           = null;
  }  // if myRank == rootRank
 
  RCP<sparse_matrix_type> A;
  typedef Export<local_ordinal_type, global_ordinal_type, node_type> export_type;
  export_type exp(gatherRowMap, rowMap);
 
  // Communicate the precise number of nonzeros per row, which was already
  // calculated above.
  typedef local_ordinal_type LO;
  typedef global_ordinal_type GO;
  typedef Tpetra::MultiVector<GO, LO, GO, node_type> mv_type_go;
  mv_type_go target_nnzPerRow(rowMap, 1);
  mv_type_go source_nnzPerRow(gatherRowMap, 1);
  Teuchos::ArrayRCP<GO> srcData = source_nnzPerRow.getDataNonConst(0);
  for (int i = 0; i < myNumRows; i++)
    srcData[i] = gatherNumEntriesPerRow[i];
  srcData = Teuchos::null;
  target_nnzPerRow.doExport(source_nnzPerRow, exp, Tpetra::INSERT);
  Teuchos::ArrayRCP<GO> targetData   = target_nnzPerRow.getDataNonConst(0);
  ArrayRCP<size_t> targetData_size_t = arcp<size_t>(targetData.size());
  for (int i = 0; i < targetData.size(); i++)
    targetData_size_t[i] = targetData[i];
 
  if (colMap.is_null()) {
    A = rcp(new sparse_matrix_type(rowMap, targetData_size_t()));
  } else {
    A = rcp(new sparse_matrix_type(rowMap, colMap, targetData_size_t()));
  }
  A->doExport(*A_proc0, exp, INSERT);
  if (callFillComplete) {
    A->fillComplete(domainMap, rangeMap);
  }
 
  // We can't get the dimensions of the matrix until after
  // fillComplete() is called.  Thus, we can't do the sanity
  // check (dimensions read from the Matrix Market data,
  // vs. dimensions reported by the CrsMatrix) unless the user
  // asked us to call fillComplete().
  //
  // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
  // what one might think it does, so you have to ask the range
  // resp. domain map for the number of rows resp. columns.
  if (callFillComplete) {
    const int numProcs = pComm->getSize();
 
    if (extraDebug && debug) {
      const global_size_t globalNumRows = rangeMap->getGlobalNumElements();
      const global_size_t globalNumCols = domainMap->getGlobalNumElements();
      if (myRank == rootRank) {
        cerr << "-- Matrix is "
             << globalNumRows << " x " << globalNumCols
             << " with " << A->getGlobalNumEntries()
             << " entries, and index base "
             << A->getIndexBase() << "." << endl;
      }
      pComm->barrier();
      for (int p = 0; p < numProcs; ++p) {
        if (myRank == p) {
          cerr << "-- Proc " << p << " owns "
               << A->getLocalNumCols() << " columns, and "
               << A->getLocalNumEntries() << " entries." << endl;
        }
        pComm->barrier();
      }
    }  // if (extraDebug && debug)
  }    // if (callFillComplete)
 
  if (debug && myRank == rootRank) {
    cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
         << endl;
  }
  return A;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readDenseFile(const std::string& filename,
                                                                             const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t& comm,
                                                                             Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& map,
                                                                             const bool tolerant,
                                                                             const bool debug,
                                                                             const bool binary) {
  std::ifstream in = MatrixMarketReader::openInFileOnRankZero(comm, filename, true, binary ? std::ios::binary : std::ios::in);
 
  return readDense(in, comm, map, tolerant, debug, binary);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
std::ifstream MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::openInFileOnRankZero(
    const trcp_tcomm_t& comm,
    const std::string& filename, const bool safe,
    std::ios_base::openmode mode) {
  // Input stream.
  std::ifstream in;
 
  // Placeholder for broadcasting in-stream state.
  int all_should_stop = false;
 
  // Try to open the in-stream on root rank.
  if (comm->getRank() == 0) {
    in.open(filename, mode);
    all_should_stop = !in && safe;
  }
 
  // Broadcast state and possibly throw.
  if (comm) Teuchos::broadcast(*comm, 0, &all_should_stop);
 
  TEUCHOS_TEST_FOR_EXCEPTION(
      all_should_stop,
      std::runtime_error,
      "Could not open input file '" + filename + "' on root rank 0.");
 
  return in;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::vector_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readVectorFile(const std::string& filename,
                                                                              const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t& comm,
                                                                              Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& map,
                                                                              const bool tolerant,
                                                                              const bool debug) {
  std::ifstream in = MatrixMarketReader::openInFileOnRankZero(comm, filename, true);
 
  return readVector(in, comm, map, tolerant, debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readDense(std::istream& in,
                                                                         const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t& comm,
                                                                         Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& map,
                                                                         const bool tolerant,
                                                                         const bool debug,
                                                                         const bool binary) {
  Teuchos::RCP<Teuchos::FancyOStream> err =
      Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
  return readDenseImpl<scalar_type>(in, comm, map, err, tolerant, debug, binary);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::vector_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readVector(std::istream& in,
                                                                          const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t& comm,
                                                                          Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& map,
                                                                          const bool tolerant,
                                                                          const bool debug) {
  Teuchos::RCP<Teuchos::FancyOStream> err =
      Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
  return readVectorImpl<scalar_type>(in, comm, map, err, tolerant, debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readMapFile(const std::string& filename,
                                                                           const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t& comm,
                                                                           const bool tolerant,
                                                                           const bool debug,
                                                                           const bool binary) {
  std::ifstream in = MatrixMarketReader::openInFileOnRankZero(comm, filename, true, binary ? std::ios::binary : std::ios::in);
 
  return readMap(in, comm, tolerant, debug, binary);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
template <class MultiVectorScalarType>
Teuchos::RCP<Tpetra::MultiVector<MultiVectorScalarType,
                                 typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_ordinal_type,
                                 typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_ordinal_type,
                                 typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::node_type>>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readDenseImpl(std::istream& in,
                                                                             const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t& comm,
                                                                             Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& map,
                                                                             const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                             const bool tolerant,
                                                                             const bool debug,
                                                                             const bool binary) {
  using std::endl;
  using Teuchos::ArrayRCP;
  using Teuchos::as;
  using Teuchos::broadcast;
  using Teuchos::outArg;
  using Teuchos::RCP;
  using Teuchos::Tuple;
  using Teuchos::MatrixMarket::Banner;
  using Teuchos::MatrixMarket::checkCommentLine;
  typedef MultiVectorScalarType ST;
  typedef local_ordinal_type LO;
  typedef global_ordinal_type GO;
  typedef node_type NT;
  typedef Teuchos::ScalarTraits<ST> STS;
  typedef typename STS::magnitudeType MT;
  typedef Teuchos::ScalarTraits<MT> STM;
  typedef Tpetra::MultiVector<ST, LO, GO, NT> MV;
 
  // Rank 0 is the only (MPI) process allowed to read from the
  // input stream.
  const int myRank = comm->getRank();
 
  if (!err.is_null()) {
    err->pushTab();
  }
  if (debug) {
    *err << myRank << ": readDenseImpl" << endl;
  }
  if (!err.is_null()) {
    err->pushTab();
  }
 
  // mfh 17 Feb 2013: It's not strictly necessary that the Comm
  // instances be identical and that the Node instances be
  // identical.  The essential condition is more complicated to
  // test and isn't the same for all Node types.  Thus, we just
  // leave it up to the user.
 
  // // If map is nonnull, check the precondition that its
  // // communicator resp. node equal comm resp. node.  Checking
  // // now avoids doing a lot of file reading before we detect the
  // // violated precondition.
  // TEUCHOS_TEST_FOR_EXCEPTION(
  //   ! map.is_null() && (map->getComm() != comm || map->getNode () != node,
  //   std::invalid_argument, "If you supply a nonnull Map, the Map's "
  //   "communicator and node must equal the supplied communicator resp. "
  //   "node.");
 
  // Process 0 will read in the matrix dimensions from the file,
  // and broadcast them to all ranks in the given communicator.
  // There are only 2 dimensions in the matrix, but we use the
  // third element of the Tuple to encode the banner's reported
  // data type: "real" == 0, "complex" == 1, and "integer" == 0
  // (same as "real").  We don't allow pattern matrices (i.e.,
  // graphs) since they only make sense for sparse data.
  Tuple<GO, 3> dims;
  dims[0] = 0;
  dims[1] = 0;
  dims[2] = 0;
 
  // Current line number in the input stream.  Only valid on
  // Proc 0.  Various calls will modify this depending on the
  // number of lines that are read from the input stream.
  size_t lineNumber = 1;
 
  // Capture errors and their messages on Proc 0.
  std::ostringstream exMsg;
  int localBannerReadSuccess = 1;
  int localDimsReadSuccess   = 1;
 
  // Only Proc 0 gets to read matrix data from the input stream.
  if (myRank == 0) {
    if (!binary) {
      if (debug) {
        *err << myRank << ": readDenseImpl: Reading banner line (dense)" << endl;
      }
 
      // The "Banner" tells you whether the input stream
      // represents a dense matrix, the symmetry type of the
      // matrix, and the type of the data it contains.
      RCP<const Banner> pBanner;
      try {
        pBanner = readBanner(in, lineNumber, tolerant, debug);
      } catch (std::exception& e) {
        exMsg << e.what();
        localBannerReadSuccess = 0;
      }
      // Make sure the input stream is the right kind of data.
      if (localBannerReadSuccess) {
        if (pBanner->matrixType() != "array") {
          exMsg << "The Matrix Market file does not contain dense matrix "
                   "data.  Its banner (first) line says that its matrix type is \""
                << pBanner->matrixType() << "\", rather that the required "
                                            "\"array\".";
          localBannerReadSuccess = 0;
        } else if (pBanner->dataType() == "pattern") {
          exMsg << "The Matrix Market file's banner (first) "
                   "line claims that the matrix's data type is \"pattern\".  This does "
                   "not make sense for a dense matrix, yet the file reports the matrix "
                   "as dense.  The only valid data types for a dense matrix are "
                   "\"real\", \"complex\", and \"integer\".";
          localBannerReadSuccess = 0;
        } else {
          // Encode the data type reported by the Banner as the
          // third element of the dimensions Tuple.
          dims[2] = encodeDataType(pBanner->dataType());
        }
      }  // if we successfully read the banner line
 
      // At this point, we've successfully read the banner line.
      // Now read the dimensions line.
      if (localBannerReadSuccess) {
        if (debug) {
          *err << myRank << ": readDenseImpl: Reading dimensions line (dense)" << endl;
        }
        // Keep reading lines from the input stream until we find
        // a non-comment line, or until we run out of lines.  The
        // latter is an error, since every "array" format Matrix
        // Market file must have a dimensions line after the
        // banner (even if the matrix has zero rows or columns, or
        // zero entries).
        std::string line;
        bool commentLine = true;
 
        while (commentLine) {
          // Test whether it is even valid to read from the input
          // stream wrapping the line.
          if (in.eof() || in.fail()) {
            exMsg << "Unable to get array dimensions line (at all) from line "
                  << lineNumber << " of input stream.  The input stream "
                  << "claims that it is "
                  << (in.eof() ? "at end-of-file." : "in a failed state.");
            localDimsReadSuccess = 0;
          } else {
            // Try to get the next line from the input stream.
            if (getline(in, line)) {
              ++lineNumber;  // We did actually read a line.
            }
            // Is the current line a comment line?  Ignore start
            // and size; they are only useful for reading the
            // actual matrix entries.  (We could use them here as
            // an optimization, but we've chosen not to.)
            size_t start = 0, size = 0;
            commentLine = checkCommentLine(line, start, size, lineNumber, tolerant);
          }  // whether we failed to read the line at all
        }    // while the line we just read is a comment line
 
        //
        // Get <numRows> <numCols> from the line we just read.
        //
        std::istringstream istr(line);
 
        // Test whether it is even valid to read from the input
        // stream wrapping the line.
        if (istr.eof() || istr.fail()) {
          exMsg << "Unable to read any data from line " << lineNumber
                << " of input; the line should contain the matrix dimensions "
                << "\"<numRows> <numCols>\".";
          localDimsReadSuccess = 0;
        } else {  // It's valid to read from the line.
          GO theNumRows = 0;
          istr >> theNumRows;  // Read in the number of rows.
          if (istr.fail()) {
            exMsg << "Failed to get number of rows from line "
                  << lineNumber << " of input; the line should contains the "
                  << "matrix dimensions \"<numRows> <numCols>\".";
            localDimsReadSuccess = 0;
          } else {                 // We successfully read the number of rows
            dims[0] = theNumRows;  // Save the number of rows
            if (istr.eof()) {      // Do we still have data to read?
              exMsg << "No more data after number of rows on line "
                    << lineNumber << " of input; the line should contain the "
                    << "matrix dimensions \"<numRows> <numCols>\".";
              localDimsReadSuccess = 0;
            } else {  // Still data left to read; read in number of columns.
              GO theNumCols = 0;
              istr >> theNumCols;  // Read in the number of columns
              if (istr.fail()) {
                exMsg << "Failed to get number of columns from line "
                      << lineNumber << " of input; the line should contain "
                      << "the matrix dimensions \"<numRows> <numCols>\".";
                localDimsReadSuccess = 0;
              } else {                 // We successfully read the number of columns
                dims[1] = theNumCols;  // Save the number of columns
              }                        // if istr.fail ()
            }                          // if istr.eof ()
          }                            // if we read the number of rows
        }                              // if the input stream wrapping the dims line was (in)valid
      }                                // if we successfully read the banner line
    }                                  // not in binary format
    else {                             // in binary format
      global_size_t numRows, numCols;
      in.read(reinterpret_cast<char*>(&numRows), sizeof(numRows));
      in.read(reinterpret_cast<char*>(&numCols), sizeof(numCols));
      TEUCHOS_TEST_FOR_EXCEPTION(
          numRows == 0 || numCols == 0,
          std::runtime_error,
          "Matrix Market dense binary header has non-positive dimensions: "
          "numRows="
              << numRows << ", numCols=" << numCols);
      dims[0] = Teuchos::as<GO>(numRows);
      dims[1] = Teuchos::as<GO>(numCols);
      if ((typeid(ST) == typeid(double)) || Teuchos::ScalarTraits<ST>::isOrdinal) {
        dims[2] = 0;
      } else if (Teuchos::ScalarTraits<ST>::isComplex) {
        dims[2] = 1;
      } else {
        TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
                                   "Unrecognized Matrix Market data type. "
                                   "We should never get here.  "
                                   "Please report this bug to the Tpetra developers.");
      }
      localDimsReadSuccess = true;
      localDimsReadSuccess = true;
    }  // in binary format
  }    // if (myRank == 0)
 
  // Broadcast the matrix dimensions, the encoded data type, and
  // whether or not Proc 0 succeeded in reading the banner and
  // dimensions.
  Tuple<GO, 5> bannerDimsReadResult;
  if (myRank == 0) {
    bannerDimsReadResult[0] = dims[0];  // numRows
    bannerDimsReadResult[1] = dims[1];  // numCols
    bannerDimsReadResult[2] = dims[2];  // encoded data type
    bannerDimsReadResult[3] = localBannerReadSuccess;
    bannerDimsReadResult[4] = localDimsReadSuccess;
  }
  // Broadcast matrix dimensions and the encoded data type from
  // Proc 0 to all the MPI processes.
  broadcast(*comm, 0, bannerDimsReadResult);
 
  TEUCHOS_TEST_FOR_EXCEPTION(
      bannerDimsReadResult[3] == 0, std::runtime_error,
      "Failed to read banner line: " << exMsg.str());
  TEUCHOS_TEST_FOR_EXCEPTION(
      bannerDimsReadResult[4] == 0, std::runtime_error,
      "Failed to read matrix dimensions line: " << exMsg.str());
  if (myRank != 0) {
    dims[0] = bannerDimsReadResult[0];
    dims[1] = bannerDimsReadResult[1];
    dims[2] = bannerDimsReadResult[2];
  }
 
  // Tpetra objects want the matrix dimensions in these types.
  const global_size_t numRows = static_cast<global_size_t>(dims[0]);
  const size_t numCols        = static_cast<size_t>(dims[1]);
 
  // Make a "Proc 0 owns everything" Map that we will use to
  // read in the multivector entries in the correct order on
  // Proc 0.  This must be a collective
  RCP<const map_type> proc0Map;  // "Proc 0 owns everything" Map
  if (map.is_null()) {
    // The user didn't supply a Map.  Make a contiguous
    // distributed Map for them, using the read-in multivector
    // dimensions.
    map                       = createUniformContigMapWithNode<LO, GO, NT>(numRows, comm);
    const size_t localNumRows = (myRank == 0) ? numRows : 0;
    // At this point, map exists and has a nonnull node.
    proc0Map = createContigMapWithNode<LO, GO, NT>(numRows, localNumRows,
                                                   comm);
  } else {  // The user supplied a Map.
    proc0Map = Details::computeGatherMap<map_type>(map, err, debug);
  }
 
  // Make a multivector X owned entirely by Proc 0.
  RCP<MV> X = createMultiVector<ST, LO, GO, NT>(proc0Map, numCols);
 
  //
  // On Proc 0, read the Matrix Market data from the input
  // stream into the multivector X.
  //
  int localReadDataSuccess = 1;
  if (myRank == 0) {
    if (!binary) {
      try {
        if (debug) {
          *err << myRank << ": readDenseImpl: Reading matrix data (dense)"
               << endl;
        }
 
        // Make sure that we can get a 1-D view of X.
        TEUCHOS_TEST_FOR_EXCEPTION(
            !X->isConstantStride(), std::logic_error,
            "Can't get a 1-D view of the entries of the MultiVector X on "
            "Process 0, because the stride between the columns of X is not "
            "constant.  This shouldn't happen because we just created X and "
            "haven't filled it in yet.  Please report this bug to the Tpetra "
            "developers.");
 
        // Get a writeable 1-D view of the entries of X.  Rank 0
        // owns all of them.  The view will expire at the end of
        // scope, so (if necessary) it will be written back to X
        // at this time.
        ArrayRCP<ST> X_view = X->get1dViewNonConst();
        TEUCHOS_TEST_FOR_EXCEPTION(
            as<global_size_t>(X_view.size()) < numRows * numCols,
            std::logic_error,
            "The view of X has size " << X_view.size() << " which is not enough to "
                                                          "accommodate the expected number of entries numRows*numCols = "
                                      << numRows << "*" << numCols << " = " << numRows * numCols << ".  "
                                                                                                    "Please report this bug to the Tpetra developers.");
        const size_t stride = X->getStride();
 
        // The third element of the dimensions Tuple encodes the data
        // type reported by the Banner: "real" == 0, "complex" == 1,
        // "integer" == 0 (same as "real"), "pattern" == 2.  We do not
        // allow dense matrices to be pattern matrices, so dims[2] ==
        // 0 or 1.  We've already checked for this above.
        const bool isComplex = (dims[2] == 1);
        size_type count = 0, curRow = 0, curCol = 0;
 
        std::string line;
        while (getline(in, line)) {
          ++lineNumber;
          // Is the current line a comment line?  If it's not,
          // line.substr(start,size) contains the data.
          size_t start = 0, size = 0;
          const bool commentLine =
              checkCommentLine(line, start, size, lineNumber, tolerant);
          if (!commentLine) {
            // Make sure we have room in which to put the new matrix
            // entry.  We check this only after checking for a
            // comment line, because there may be one or more
            // comment lines at the end of the file.  In tolerant
            // mode, we simply ignore any extra data.
            if (count >= X_view.size()) {
              if (tolerant) {
                break;
              } else {
                TEUCHOS_TEST_FOR_EXCEPTION(
                    count >= X_view.size(),
                    std::runtime_error,
                    "The Matrix Market input stream has more data in it than "
                    "its metadata reported.  Current line number is "
                        << lineNumber << ".");
              }
            }
 
            // mfh 19 Dec 2012: Ignore everything up to the initial
            // colon.  writeDense() has the option to print out the
            // global row index in front of each entry, followed by
            // a colon and space.
            {
              const size_t pos = line.substr(start, size).find(':');
              if (pos != std::string::npos) {
                start = pos + 1;
              }
            }
            std::istringstream istr(line.substr(start, size));
            // Does the line contain anything at all?  Can we
            // safely read from the input stream wrapping the
            // line?
            if (istr.eof() || istr.fail()) {
              // In tolerant mode, simply ignore the line.
              if (tolerant) {
                break;
              }
              // We repeat the full test here so the exception
              // message is more informative.
              TEUCHOS_TEST_FOR_EXCEPTION(
                  !tolerant && (istr.eof() || istr.fail()),
                  std::runtime_error,
                  "Line " << lineNumber << " of the Matrix Market file is "
                                           "empty, or we cannot read from it for some other reason.");
            }
            // Current matrix entry to read in.
            ST val = STS::zero();
            // Real and imaginary parts of the current matrix entry.
            // The imaginary part is zero if the matrix is real-valued.
            MT real = STM::zero(), imag = STM::zero();
 
            // isComplex refers to the input stream's data, not to
            // the scalar type S.  It's OK to read real-valued
            // data into a matrix storing complex-valued data; in
            // that case, all entries' imaginary parts are zero.
            if (isComplex) {
              // STS::real() and STS::imag() return a copy of
              // their respective components, not a writeable
              // reference.  Otherwise we could just assign to
              // them using the istream extraction operator (>>).
              // That's why we have separate magnitude type "real"
              // and "imag" variables.
 
              // Attempt to read the real part of the current entry.
              istr >> real;
              if (istr.fail()) {
                TEUCHOS_TEST_FOR_EXCEPTION(
                    !tolerant && istr.eof(), std::runtime_error,
                    "Failed to get the real part of a complex-valued matrix "
                    "entry from line "
                        << lineNumber << " of the Matrix Market "
                                         "file.");
                // In tolerant mode, just skip bad lines.
                if (tolerant) {
                  break;
                }
              } else if (istr.eof()) {
                TEUCHOS_TEST_FOR_EXCEPTION(
                    !tolerant && istr.eof(), std::runtime_error,
                    "Missing imaginary part of a complex-valued matrix entry "
                    "on line "
                        << lineNumber << " of the Matrix Market file.");
                // In tolerant mode, let any missing imaginary part be 0.
              } else {
                // Attempt to read the imaginary part of the current
                // matrix entry.
                istr >> imag;
                TEUCHOS_TEST_FOR_EXCEPTION(
                    !tolerant && istr.fail(), std::runtime_error,
                    "Failed to get the imaginary part of a complex-valued "
                    "matrix entry from line "
                        << lineNumber << " of the "
                                         "Matrix Market file.");
                // In tolerant mode, let any missing or corrupted
                // imaginary part be 0.
              }
            } else {  // Matrix Market file contains real-valued data.
              // Attempt to read the current matrix entry.
              istr >> real;
              TEUCHOS_TEST_FOR_EXCEPTION(
                  !tolerant && istr.fail(), std::runtime_error,
                  "Failed to get a real-valued matrix entry from line "
                      << lineNumber << " of the Matrix Market file.");
              // In tolerant mode, simply ignore the line if
              // we failed to read a matrix entry.
              if (istr.fail() && tolerant) {
                break;
              }
            }
            // In tolerant mode, we simply let pass through whatever
            // data we got.
            TEUCHOS_TEST_FOR_EXCEPTION(
                !tolerant && istr.fail(), std::runtime_error,
                "Failed to read matrix data from line " << lineNumber
                                                        << " of the Matrix Market file.");
 
            // Assign val = ST(real, imag).
            Teuchos::MatrixMarket::details::assignScalar<ST>(val, real, imag);
 
            curRow                           = count % numRows;
            curCol                           = count / numRows;
            X_view[curRow + curCol * stride] = val;
            ++count;
          }  // if not a comment line
        }    // while there are still lines in the file, get the next one
 
        TEUCHOS_TEST_FOR_EXCEPTION(
            !tolerant && static_cast<global_size_t>(count) < numRows * numCols,
            std::runtime_error,
            "The Matrix Market metadata reports that the dense matrix is "
                << numRows << " x " << numCols << ", and thus has "
                << numRows * numCols << " total entries, but we only found " << count
                << " entr" << (count == 1 ? "y" : "ies") << " in the file.");
      } catch (std::exception& e) {
        exMsg << e.what();
        localReadDataSuccess = 0;
      }
    }  // not binary file
    else {
      if (debug) {
        *err << myRank << ": readDenseImpl: Reading matrix data (dense)"
             << endl;
      }
 
      // Make sure that we can get a 1-D view of X.
      TEUCHOS_TEST_FOR_EXCEPTION(
          !X->isConstantStride(), std::logic_error,
          "Can't get a 1-D view of the entries of the MultiVector X on "
          "Process 0, because the stride between the columns of X is not "
          "constant.  This shouldn't happen because we just created X and "
          "haven't filled it in yet.  Please report this bug to the Tpetra "
          "developers.");
 
      // Get a writeable 1-D view of the entries of X.  Rank 0
      // owns all of them.  The view will expire at the end of
      // scope, so (if necessary) it will be written back to X
      // at this time.
      auto X_view = X->getLocalViewHost(Access::OverwriteAll);
 
      TEUCHOS_TEST_FOR_EXCEPTION(
          as<global_size_t>(X_view.extent(0)) < numRows,
          std::logic_error,
          "The view of X has " << X_view.extent(0) << " rows which is not enough to "
                                                      "accommodate the expected number of entries numRows = "
                               << numRows << ".  "
                                             "Please report this bug to the Tpetra developers.");
      TEUCHOS_TEST_FOR_EXCEPTION(
          as<global_size_t>(X_view.extent(1)) < numCols,
          std::logic_error,
          "The view of X has " << X_view.extent(1) << " colums which is not enough to "
                                                      "accommodate the expected number of entries numRows = "
                               << numCols << ".  "
                                             "Please report this bug to the Tpetra developers.");
 
      for (size_t curRow = 0; curRow < numRows; ++curRow) {
        for (size_t curCol = 0; curCol < numCols; ++curCol) {
          if (Teuchos::ScalarTraits<ST>::isOrdinal) {
            global_size_t val;
            in.read(reinterpret_cast<char*>(&val), sizeof(val));
            X_view(curRow, curCol) = val;
          } else {
            double val;
            in.read(reinterpret_cast<char*>(&val), sizeof(val));
            X_view(curRow, curCol) = val;
          }
        }
      }
    }  // binary
  }    // if (myRank == 0)
 
  if (debug) {
    *err << myRank << ": readDenseImpl: done reading data" << endl;
  }
 
  // Synchronize on whether Proc 0 successfully read the data.
  int globalReadDataSuccess = localReadDataSuccess;
  broadcast(*comm, 0, outArg(globalReadDataSuccess));
  TEUCHOS_TEST_FOR_EXCEPTION(
      globalReadDataSuccess == 0, std::runtime_error,
      "Failed to read the multivector's data: " << exMsg.str());
 
  // If there's only one MPI process and the user didn't supply
  // a Map (i.e., pMap is null), we're done.  Set pMap to the
  // Map used to distribute X, and return X.
  if (comm->getSize() == 1 && map.is_null()) {
    map = proc0Map;
    if (!err.is_null()) {
      err->popTab();
    }
    if (debug) {
      *err << myRank << ": readDenseImpl: done" << endl;
    }
    if (!err.is_null()) {
      err->popTab();
    }
    return X;
  }
 
  if (debug) {
    *err << myRank << ": readDenseImpl: Creating target MV" << endl;
  }
 
  // Make a multivector Y with the distributed map pMap.
  RCP<MV> Y = createMultiVector<ST, LO, GO, NT>(map, numCols);
 
  if (debug) {
    *err << myRank << ": readDenseImpl: Creating Export" << endl;
  }
 
  // Make an Export object that will export X to Y.  First
  // argument is the source map, second argument is the target
  // map.
  Export<LO, GO, NT> exporter(proc0Map, map, err);
 
  if (debug) {
    *err << myRank << ": readDenseImpl: Exporting" << endl;
  }
  // Export X into Y.
  Y->doExport(*X, exporter, INSERT);
 
  if (!err.is_null()) {
    err->popTab();
  }
  if (debug) {
    *err << myRank << ": readDenseImpl: done" << endl;
  }
  if (!err.is_null()) {
    err->popTab();
  }
 
  // Y is distributed over all process(es) in the communicator.
  return Y;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
template <class VectorScalarType>
Teuchos::RCP<Tpetra::Vector<VectorScalarType,
                            typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_ordinal_type,
                            typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_ordinal_type,
                            typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::node_type>>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readVectorImpl(std::istream& in,
                                                                              const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t& comm,
                                                                              Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& map,
                                                                              const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                              const bool tolerant,
                                                                              const bool debug) {
  using std::endl;
  using Teuchos::as;
  using Teuchos::broadcast;
  using Teuchos::outArg;
  using Teuchos::RCP;
  using Teuchos::Tuple;
  using Teuchos::MatrixMarket::Banner;
  using Teuchos::MatrixMarket::checkCommentLine;
  typedef VectorScalarType ST;
  typedef local_ordinal_type LO;
  typedef global_ordinal_type GO;
  typedef node_type NT;
  typedef Teuchos::ScalarTraits<ST> STS;
  typedef typename STS::magnitudeType MT;
  typedef Teuchos::ScalarTraits<MT> STM;
  typedef Tpetra::Vector<ST, LO, GO, NT> MV;
 
  // Rank 0 is the only (MPI) process allowed to read from the
  // input stream.
  const int myRank = comm->getRank();
 
  if (!err.is_null()) {
    err->pushTab();
  }
  if (debug) {
    *err << myRank << ": readVectorImpl" << endl;
  }
  if (!err.is_null()) {
    err->pushTab();
  }
 
  // mfh 17 Feb 2013: It's not strictly necessary that the Comm
  // instances be identical and that the Node instances be
  // identical.  The essential condition is more complicated to
  // test and isn't the same for all Node types.  Thus, we just
  // leave it up to the user.
 
  // // If map is nonnull, check the precondition that its
  // // communicator resp. node equal comm resp. node.  Checking
  // // now avoids doing a lot of file reading before we detect the
  // // violated precondition.
  // TEUCHOS_TEST_FOR_EXCEPTION(
  //   ! map.is_null() && (map->getComm() != comm || map->getNode () != node,
  //   std::invalid_argument, "If you supply a nonnull Map, the Map's "
  //   "communicator and node must equal the supplied communicator resp. "
  //   "node.");
 
  // Process 0 will read in the matrix dimensions from the file,
  // and broadcast them to all ranks in the given communicator.
  // There are only 2 dimensions in the matrix, but we use the
  // third element of the Tuple to encode the banner's reported
  // data type: "real" == 0, "complex" == 1, and "integer" == 0
  // (same as "real").  We don't allow pattern matrices (i.e.,
  // graphs) since they only make sense for sparse data.
  Tuple<GO, 3> dims;
  dims[0] = 0;
  dims[1] = 0;
 
  // Current line number in the input stream.  Only valid on
  // Proc 0.  Various calls will modify this depending on the
  // number of lines that are read from the input stream.
  size_t lineNumber = 1;
 
  // Capture errors and their messages on Proc 0.
  std::ostringstream exMsg;
  int localBannerReadSuccess = 1;
  int localDimsReadSuccess   = 1;
 
  // Only Proc 0 gets to read matrix data from the input stream.
  if (myRank == 0) {
    if (debug) {
      *err << myRank << ": readVectorImpl: Reading banner line (dense)" << endl;
    }
 
    // The "Banner" tells you whether the input stream
    // represents a dense matrix, the symmetry type of the
    // matrix, and the type of the data it contains.
    RCP<const Banner> pBanner;
    try {
      pBanner = readBanner(in, lineNumber, tolerant, debug);
    } catch (std::exception& e) {
      exMsg << e.what();
      localBannerReadSuccess = 0;
    }
    // Make sure the input stream is the right kind of data.
    if (localBannerReadSuccess) {
      if (pBanner->matrixType() != "array") {
        exMsg << "The Matrix Market file does not contain dense matrix "
                 "data.  Its banner (first) line says that its matrix type is \""
              << pBanner->matrixType() << "\", rather that the required "
                                          "\"array\".";
        localBannerReadSuccess = 0;
      } else if (pBanner->dataType() == "pattern") {
        exMsg << "The Matrix Market file's banner (first) "
                 "line claims that the matrix's data type is \"pattern\".  This does "
                 "not make sense for a dense matrix, yet the file reports the matrix "
                 "as dense.  The only valid data types for a dense matrix are "
                 "\"real\", \"complex\", and \"integer\".";
        localBannerReadSuccess = 0;
      } else {
        // Encode the data type reported by the Banner as the
        // third element of the dimensions Tuple.
        dims[2] = encodeDataType(pBanner->dataType());
      }
    }  // if we successfully read the banner line
 
    // At this point, we've successfully read the banner line.
    // Now read the dimensions line.
    if (localBannerReadSuccess) {
      if (debug) {
        *err << myRank << ": readVectorImpl: Reading dimensions line (dense)" << endl;
      }
      // Keep reading lines from the input stream until we find
      // a non-comment line, or until we run out of lines.  The
      // latter is an error, since every "array" format Matrix
      // Market file must have a dimensions line after the
      // banner (even if the matrix has zero rows or columns, or
      // zero entries).
      std::string line;
      bool commentLine = true;
 
      while (commentLine) {
        // Test whether it is even valid to read from the input
        // stream wrapping the line.
        if (in.eof() || in.fail()) {
          exMsg << "Unable to get array dimensions line (at all) from line "
                << lineNumber << " of input stream.  The input stream "
                << "claims that it is "
                << (in.eof() ? "at end-of-file." : "in a failed state.");
          localDimsReadSuccess = 0;
        } else {
          // Try to get the next line from the input stream.
          if (getline(in, line)) {
            ++lineNumber;  // We did actually read a line.
          }
          // Is the current line a comment line?  Ignore start
          // and size; they are only useful for reading the
          // actual matrix entries.  (We could use them here as
          // an optimization, but we've chosen not to.)
          size_t start = 0, size = 0;
          commentLine = checkCommentLine(line, start, size, lineNumber, tolerant);
        }  // whether we failed to read the line at all
      }    // while the line we just read is a comment line
 
      //
      // Get <numRows> <numCols> from the line we just read.
      //
      std::istringstream istr(line);
 
      // Test whether it is even valid to read from the input
      // stream wrapping the line.
      if (istr.eof() || istr.fail()) {
        exMsg << "Unable to read any data from line " << lineNumber
              << " of input; the line should contain the matrix dimensions "
              << "\"<numRows> <numCols>\".";
        localDimsReadSuccess = 0;
      } else {  // It's valid to read from the line.
        GO theNumRows = 0;
        istr >> theNumRows;  // Read in the number of rows.
        if (istr.fail()) {
          exMsg << "Failed to get number of rows from line "
                << lineNumber << " of input; the line should contains the "
                << "matrix dimensions \"<numRows> <numCols>\".";
          localDimsReadSuccess = 0;
        } else {                 // We successfully read the number of rows
          dims[0] = theNumRows;  // Save the number of rows
          if (istr.eof()) {      // Do we still have data to read?
            exMsg << "No more data after number of rows on line "
                  << lineNumber << " of input; the line should contain the "
                  << "matrix dimensions \"<numRows> <numCols>\".";
            localDimsReadSuccess = 0;
          } else {  // Still data left to read; read in number of columns.
            GO theNumCols = 0;
            istr >> theNumCols;  // Read in the number of columns
            if (istr.fail()) {
              exMsg << "Failed to get number of columns from line "
                    << lineNumber << " of input; the line should contain "
                    << "the matrix dimensions \"<numRows> <numCols>\".";
              localDimsReadSuccess = 0;
            } else {                 // We successfully read the number of columns
              dims[1] = theNumCols;  // Save the number of columns
            }                        // if istr.fail ()
          }                          // if istr.eof ()
        }                            // if we read the number of rows
      }                              // if the input stream wrapping the dims line was (in)valid
    }                                // if we successfully read the banner line
  }                                  // if (myRank == 0)
 
  // Check if file has a Vector
  if (dims[1] != 1) {
    exMsg << "File does not contain a 1D Vector";
    localDimsReadSuccess = 0;
  }
 
  // Broadcast the matrix dimensions, the encoded data type, and
  // whether or not Proc 0 succeeded in reading the banner and
  // dimensions.
  Tuple<GO, 5> bannerDimsReadResult;
  if (myRank == 0) {
    bannerDimsReadResult[0] = dims[0];  // numRows
    bannerDimsReadResult[1] = dims[1];  // numCols
    bannerDimsReadResult[2] = dims[2];  // encoded data type
    bannerDimsReadResult[3] = localBannerReadSuccess;
    bannerDimsReadResult[4] = localDimsReadSuccess;
  }
 
  // Broadcast matrix dimensions and the encoded data type from
  // Proc 0 to all the MPI processes.
  broadcast(*comm, 0, bannerDimsReadResult);
 
  TEUCHOS_TEST_FOR_EXCEPTION(
      bannerDimsReadResult[3] == 0, std::runtime_error,
      "Failed to read banner line: " << exMsg.str());
  TEUCHOS_TEST_FOR_EXCEPTION(
      bannerDimsReadResult[4] == 0, std::runtime_error,
      "Failed to read matrix dimensions line: " << exMsg.str());
  if (myRank != 0) {
    dims[0] = bannerDimsReadResult[0];
    dims[1] = bannerDimsReadResult[1];
    dims[2] = bannerDimsReadResult[2];
  }
 
  // Tpetra objects want the matrix dimensions in these types.
  const global_size_t numRows = static_cast<global_size_t>(dims[0]);
  const size_t numCols        = static_cast<size_t>(dims[1]);
 
  // Make a "Proc 0 owns everything" Map that we will use to
  // read in the multivector entries in the correct order on
  // Proc 0.  This must be a collective
  RCP<const map_type> proc0Map;  // "Proc 0 owns everything" Map
  if (map.is_null()) {
    // The user didn't supply a Map.  Make a contiguous
    // distributed Map for them, using the read-in multivector
    // dimensions.
    map                       = createUniformContigMapWithNode<LO, GO, NT>(numRows, comm);
    const size_t localNumRows = (myRank == 0) ? numRows : 0;
    // At this point, map exists and has a nonnull node.
    proc0Map = createContigMapWithNode<LO, GO, NT>(numRows, localNumRows,
                                                   comm);
  } else {  // The user supplied a Map.
    proc0Map = Details::computeGatherMap<map_type>(map, err, debug);
  }
 
  // Make a multivector X owned entirely by Proc 0.
  RCP<MV> X = createVector<ST, LO, GO, NT>(proc0Map);
 
  //
  // On Proc 0, read the Matrix Market data from the input
  // stream into the multivector X.
  //
  int localReadDataSuccess = 1;
  if (myRank == 0) {
    try {
      if (debug) {
        *err << myRank << ": readVectorImpl: Reading matrix data (dense)"
             << endl;
      }
 
      // Make sure that we can get a 1-D view of X.
      TEUCHOS_TEST_FOR_EXCEPTION(
          !X->isConstantStride(), std::logic_error,
          "Can't get a 1-D view of the entries of the MultiVector X on "
          "Process 0, because the stride between the columns of X is not "
          "constant.  This shouldn't happen because we just created X and "
          "haven't filled it in yet.  Please report this bug to the Tpetra "
          "developers.");
 
      // Get a writeable 1-D view of the entries of X.  Rank 0
      // owns all of them.  The view will expire at the end of
      // scope, so (if necessary) it will be written back to X
      // at this time.
      Teuchos::ArrayRCP<ST> X_view = X->get1dViewNonConst();
      TEUCHOS_TEST_FOR_EXCEPTION(
          as<global_size_t>(X_view.size()) < numRows * numCols,
          std::logic_error,
          "The view of X has size " << X_view << " which is not enough to "
                                                 "accommodate the expected number of entries numRows*numCols = "
                                    << numRows << "*" << numCols << " = " << numRows * numCols << ".  "
                                                                                                  "Please report this bug to the Tpetra developers.");
      const size_t stride = X->getStride();
 
      // The third element of the dimensions Tuple encodes the data
      // type reported by the Banner: "real" == 0, "complex" == 1,
      // "integer" == 0 (same as "real"), "pattern" == 2.  We do not
      // allow dense matrices to be pattern matrices, so dims[2] ==
      // 0 or 1.  We've already checked for this above.
      const bool isComplex = (dims[2] == 1);
      size_type count = 0, curRow = 0, curCol = 0;
 
      std::string line;
      while (getline(in, line)) {
        ++lineNumber;
        // Is the current line a comment line?  If it's not,
        // line.substr(start,size) contains the data.
        size_t start = 0, size = 0;
        const bool commentLine =
            checkCommentLine(line, start, size, lineNumber, tolerant);
        if (!commentLine) {
          // Make sure we have room in which to put the new matrix
          // entry.  We check this only after checking for a
          // comment line, because there may be one or more
          // comment lines at the end of the file.  In tolerant
          // mode, we simply ignore any extra data.
          if (count >= X_view.size()) {
            if (tolerant) {
              break;
            } else {
              TEUCHOS_TEST_FOR_EXCEPTION(
                  count >= X_view.size(),
                  std::runtime_error,
                  "The Matrix Market input stream has more data in it than "
                  "its metadata reported.  Current line number is "
                      << lineNumber << ".");
            }
          }
 
          // mfh 19 Dec 2012: Ignore everything up to the initial
          // colon.  writeDense() has the option to print out the
          // global row index in front of each entry, followed by
          // a colon and space.
          {
            const size_t pos = line.substr(start, size).find(':');
            if (pos != std::string::npos) {
              start = pos + 1;
            }
          }
          std::istringstream istr(line.substr(start, size));
          // Does the line contain anything at all?  Can we
          // safely read from the input stream wrapping the
          // line?
          if (istr.eof() || istr.fail()) {
            // In tolerant mode, simply ignore the line.
            if (tolerant) {
              break;
            }
            // We repeat the full test here so the exception
            // message is more informative.
            TEUCHOS_TEST_FOR_EXCEPTION(
                !tolerant && (istr.eof() || istr.fail()),
                std::runtime_error,
                "Line " << lineNumber << " of the Matrix Market file is "
                                         "empty, or we cannot read from it for some other reason.");
          }
          // Current matrix entry to read in.
          ST val = STS::zero();
          // Real and imaginary parts of the current matrix entry.
          // The imaginary part is zero if the matrix is real-valued.
          MT real = STM::zero(), imag = STM::zero();
 
          // isComplex refers to the input stream's data, not to
          // the scalar type S.  It's OK to read real-valued
          // data into a matrix storing complex-valued data; in
          // that case, all entries' imaginary parts are zero.
          if (isComplex) {
            // STS::real() and STS::imag() return a copy of
            // their respective components, not a writeable
            // reference.  Otherwise we could just assign to
            // them using the istream extraction operator (>>).
            // That's why we have separate magnitude type "real"
            // and "imag" variables.
 
            // Attempt to read the real part of the current entry.
            istr >> real;
            if (istr.fail()) {
              TEUCHOS_TEST_FOR_EXCEPTION(
                  !tolerant && istr.eof(), std::runtime_error,
                  "Failed to get the real part of a complex-valued matrix "
                  "entry from line "
                      << lineNumber << " of the Matrix Market "
                                       "file.");
              // In tolerant mode, just skip bad lines.
              if (tolerant) {
                break;
              }
            } else if (istr.eof()) {
              TEUCHOS_TEST_FOR_EXCEPTION(
                  !tolerant && istr.eof(), std::runtime_error,
                  "Missing imaginary part of a complex-valued matrix entry "
                  "on line "
                      << lineNumber << " of the Matrix Market file.");
              // In tolerant mode, let any missing imaginary part be 0.
            } else {
              // Attempt to read the imaginary part of the current
              // matrix entry.
              istr >> imag;
              TEUCHOS_TEST_FOR_EXCEPTION(
                  !tolerant && istr.fail(), std::runtime_error,
                  "Failed to get the imaginary part of a complex-valued "
                  "matrix entry from line "
                      << lineNumber << " of the "
                                       "Matrix Market file.");
              // In tolerant mode, let any missing or corrupted
              // imaginary part be 0.
            }
          } else {  // Matrix Market file contains real-valued data.
            // Attempt to read the current matrix entry.
            istr >> real;
            TEUCHOS_TEST_FOR_EXCEPTION(
                !tolerant && istr.fail(), std::runtime_error,
                "Failed to get a real-valued matrix entry from line "
                    << lineNumber << " of the Matrix Market file.");
            // In tolerant mode, simply ignore the line if
            // we failed to read a matrix entry.
            if (istr.fail() && tolerant) {
              break;
            }
          }
          // In tolerant mode, we simply let pass through whatever
          // data we got.
          TEUCHOS_TEST_FOR_EXCEPTION(
              !tolerant && istr.fail(), std::runtime_error,
              "Failed to read matrix data from line " << lineNumber
                                                      << " of the Matrix Market file.");
 
          // Assign val = ST(real, imag).
          Teuchos::MatrixMarket::details::assignScalar<ST>(val, real, imag);
 
          curRow                           = count % numRows;
          curCol                           = count / numRows;
          X_view[curRow + curCol * stride] = val;
          ++count;
        }  // if not a comment line
      }    // while there are still lines in the file, get the next one
 
      TEUCHOS_TEST_FOR_EXCEPTION(
          !tolerant && static_cast<global_size_t>(count) < numRows * numCols,
          std::runtime_error,
          "The Matrix Market metadata reports that the dense matrix is "
              << numRows << " x " << numCols << ", and thus has "
              << numRows * numCols << " total entries, but we only found " << count
              << " entr" << (count == 1 ? "y" : "ies") << " in the file.");
    } catch (std::exception& e) {
      exMsg << e.what();
      localReadDataSuccess = 0;
    }
  }  // if (myRank == 0)
 
  if (debug) {
    *err << myRank << ": readVectorImpl: done reading data" << endl;
  }
 
  // Synchronize on whether Proc 0 successfully read the data.
  int globalReadDataSuccess = localReadDataSuccess;
  broadcast(*comm, 0, outArg(globalReadDataSuccess));
  TEUCHOS_TEST_FOR_EXCEPTION(
      globalReadDataSuccess == 0, std::runtime_error,
      "Failed to read the multivector's data: " << exMsg.str());
 
  // If there's only one MPI process and the user didn't supply
  // a Map (i.e., pMap is null), we're done.  Set pMap to the
  // Map used to distribute X, and return X.
  if (comm->getSize() == 1 && map.is_null()) {
    map = proc0Map;
    if (!err.is_null()) {
      err->popTab();
    }
    if (debug) {
      *err << myRank << ": readVectorImpl: done" << endl;
    }
    if (!err.is_null()) {
      err->popTab();
    }
    return X;
  }
 
  if (debug) {
    *err << myRank << ": readVectorImpl: Creating target MV" << endl;
  }
 
  // Make a multivector Y with the distributed map pMap.
  RCP<MV> Y = createVector<ST, LO, GO, NT>(map);
 
  if (debug) {
    *err << myRank << ": readVectorImpl: Creating Export" << endl;
  }
 
  // Make an Export object that will export X to Y.  First
  // argument is the source map, second argument is the target
  // map.
  Export<LO, GO, NT> exporter(proc0Map, map, err);
 
  if (debug) {
    *err << myRank << ": readVectorImpl: Exporting" << endl;
  }
  // Export X into Y.
  Y->doExport(*X, exporter, INSERT);
 
  if (!err.is_null()) {
    err->popTab();
  }
  if (debug) {
    *err << myRank << ": readVectorImpl: done" << endl;
  }
  if (!err.is_null()) {
    err->popTab();
  }
 
  // Y is distributed over all process(es) in the communicator.
  return Y;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readMap(std::istream& in,
                                                                       const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t& comm,
                                                                       const bool tolerant,
                                                                       const bool debug,
                                                                       const bool binary) {
  Teuchos::RCP<Teuchos::FancyOStream> err =
      Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
  return readMap(in, comm, err, tolerant, debug, binary);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readMap(std::istream& in,
                                                                       const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t& comm,
                                                                       const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                       const bool tolerant,
                                                                       const bool debug,
                                                                       const bool binary) {
  using std::endl;
  using Teuchos::arcp;
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::as;
  using Teuchos::broadcast;
  using Teuchos::Comm;
  using Teuchos::CommRequest;
  using Teuchos::inOutArg;
  using Teuchos::ireceive;
  using Teuchos::isend;
  using Teuchos::outArg;
  using Teuchos::RCP;
  using Teuchos::receive;
  using Teuchos::REDUCE_MIN;
  using Teuchos::reduceAll;
  using Teuchos::SerialComm;
  using Teuchos::toString;
  using Teuchos::wait;
  typedef Tpetra::global_size_t GST;
  typedef ptrdiff_t int_type;  // Can hold int and GO
  typedef local_ordinal_type LO;
  typedef global_ordinal_type GO;
  typedef node_type NT;
  typedef Tpetra::MultiVector<GO, LO, GO, NT> MV;
 
  const int numProcs = comm->getSize();
  const int myRank   = comm->getRank();
 
  if (!err.is_null()) {
    err->pushTab();
  }
  if (debug) {
    std::ostringstream os;
    os << myRank << ": readMap: " << endl;
    *err << os.str();
  }
  if (!err.is_null()) {
    err->pushTab();
  }
 
  // Tag for receive-size / send-size messages.  writeMap used
  // tags 1337 and 1338; we count up from there.
  const int sizeTag = 1339;
  // Tag for receive-data / send-data messages.
  const int dataTag = 1340;
 
  // These are for sends on Process 0, and for receives on all
  // other processes.  sizeReq is for the {receive,send}-size
  // message, and dataReq is for the message containing the
  // actual GIDs to belong to the receiving process.
  RCP<CommRequest<int>> sizeReq;
  RCP<CommRequest<int>> dataReq;
 
  // Each process will have to receive the number of GIDs to
  // expect.  Thus, we can post the receives now, and cancel
  // them if something should go wrong in the meantime.
  ArrayRCP<int_type> numGidsToRecv(1);
  numGidsToRecv[0] = 0;
  if (myRank != 0) {
    sizeReq = ireceive<int, int_type>(numGidsToRecv, 0, sizeTag, *comm);
  }
 
  int readSuccess = 1;
  std::ostringstream exMsg;
  RCP<MV> data;  // Will only be valid on Proc 0
  if (myRank == 0) {
    // If we want to reuse readDenseImpl, we have to make a
    // communicator that only contains Proc 0.  Otherwise,
    // readDenseImpl will redistribute the data to all
    // processes.  While we eventually want that, neither we nor
    // readDenseImpl know the correct Map to use at the moment.
    // That depends on the second column of the multivector.
    RCP<const Comm<int>> proc0Comm(new SerialComm<int>());
    try {
      RCP<const map_type> dataMap;
      // This is currently the only place where we use the
      // 'tolerant' argument.  Later, if we want to be clever,
      // we could have tolerant mode allow PIDs out of order.
      data = readDenseImpl<GO>(in, proc0Comm, dataMap, err, tolerant, debug, binary);
      (void)dataMap;  // Silence "unused" warnings
      if (data.is_null()) {
        readSuccess = 0;
        exMsg << "readDenseImpl() returned null." << endl;
      }
    } catch (std::exception& e) {
      readSuccess = 0;
      exMsg << e.what() << endl;
    }
  }
 
  // Map from PID to all the GIDs for that PID.
  // Only populated on Process 0.
  std::map<int, Array<GO>> pid2gids;
 
  // The index base must be the global minimum GID.
  // We will compute this on Process 0 and broadcast,
  // so that all processes can set up the Map.
  int_type globalNumGIDs = 0;
 
  // The index base must be the global minimum GID.
  // We will compute this on Process 0 and broadcast,
  // so that all processes can set up the Map.
  GO indexBase = 0;
 
  // Process 0: If the above read of the MultiVector succeeded,
  // extract the GIDs and PIDs into pid2gids, and find the
  // global min GID.
  if (myRank == 0 && readSuccess == 1) {
    if (data->getNumVectors() == 2) {  // Map format 1.0
      ArrayRCP<const GO> GIDs = data->getData(0);
      ArrayRCP<const GO> PIDs = data->getData(1);  // convert to int
      globalNumGIDs           = GIDs.size();
 
      // Start computing the global min GID, while collecting
      // the GIDs for each PID.
      if (globalNumGIDs > 0) {
        const int pid = static_cast<int>(PIDs[0]);
 
        if (pid < 0 || pid >= numProcs) {
          readSuccess = 0;
          exMsg << "Tpetra::MatrixMarket::readMap: "
                << "Encountered invalid PID " << pid << "." << endl;
        } else {
          const GO gid = GIDs[0];
          pid2gids[pid].push_back(gid);
          indexBase = gid;  // the current min GID
        }
      }
      if (readSuccess == 1) {
        // Collect the rest of the GIDs for each PID, and compute
        // the global min GID.
        for (size_type k = 1; k < globalNumGIDs; ++k) {
          const int pid = static_cast<int>(PIDs[k]);
          if (pid < 0 || pid >= numProcs) {
            readSuccess = 0;
            exMsg << "Tpetra::MatrixMarket::readMap: "
                  << "Encountered invalid PID " << pid << "." << endl;
          } else {
            const int_type gid = GIDs[k];
            pid2gids[pid].push_back(gid);
            if (gid < indexBase) {
              indexBase = gid;  // the current min GID
            }
          }
        }
      }
    } else if (data->getNumVectors() == 1) {  // Map format 2.0
      if (data->getGlobalLength() % 2 != static_cast<GST>(0)) {
        readSuccess = 0;
        exMsg << "Tpetra::MatrixMarket::readMap: Input data has the "
                 "wrong format (for Map format 2.0).  The global number of rows "
                 "in the MultiVector must be even (divisible by 2)."
              << endl;
      } else {
        ArrayRCP<const GO> theData = data->getData(0);
        globalNumGIDs              = static_cast<GO>(data->getGlobalLength()) /
                        static_cast<GO>(2);
 
        // Start computing the global min GID, while
        // collecting the GIDs for each PID.
        if (globalNumGIDs > 0) {
          const int pid = static_cast<int>(theData[1]);
          if (pid < 0 || pid >= numProcs) {
            readSuccess = 0;
            exMsg << "Tpetra::MatrixMarket::readMap: "
                  << "Encountered invalid PID " << pid << "." << endl;
          } else {
            const GO gid = theData[0];
            pid2gids[pid].push_back(gid);
            indexBase = gid;  // the current min GID
          }
        }
        // Collect the rest of the GIDs for each PID, and
        // compute the global min GID.
        for (int_type k = 1; k < globalNumGIDs; ++k) {
          const int pid = static_cast<int>(theData[2 * k + 1]);
          if (pid < 0 || pid >= numProcs) {
            readSuccess = 0;
            exMsg << "Tpetra::MatrixMarket::readMap: "
                  << "Encountered invalid PID " << pid << "." << endl;
          } else {
            const GO gid = theData[2 * k];
            pid2gids[pid].push_back(gid);
            if (gid < indexBase) {
              indexBase = gid;  // the current min GID
            }
          }
        }  // for each GID
      }    // if the amount of data is correct
    } else {
      readSuccess = 0;
      exMsg << "Tpetra::MatrixMarket::readMap: Input data must have "
               "either 1 column (for the new Map format 2.0) or 2 columns (for "
               "the old Map format 1.0).";
    }
  }  // myRank is zero
 
  // Broadcast the indexBase, the global number of GIDs, and the
  // current success status.  Use int_type for all of these.
  {
    int_type readResults[3];
    readResults[0] = static_cast<int_type>(indexBase);
    readResults[1] = static_cast<int_type>(globalNumGIDs);
    readResults[2] = static_cast<int_type>(readSuccess);
    broadcast<int, int_type>(*comm, 0, 3, readResults);
 
    indexBase     = static_cast<GO>(readResults[0]);
    globalNumGIDs = static_cast<int_type>(readResults[1]);
    readSuccess   = static_cast<int>(readResults[2]);
  }
 
  // Unwinding the stack will invoke sizeReq's destructor, which
  // will cancel the receive-size request on all processes that
  // posted it.
  TEUCHOS_TEST_FOR_EXCEPTION(
      readSuccess != 1, std::runtime_error,
      "Tpetra::MatrixMarket::readMap: Reading the Map failed with the "
      "following exception message: "
          << exMsg.str());
 
  if (myRank == 0) {
    // Proc 0: Send each process' number of GIDs to that process.
    for (int p = 1; p < numProcs; ++p) {
      ArrayRCP<int_type> numGidsToSend(1);
 
      auto it = pid2gids.find(p);
      if (it == pid2gids.end()) {
        numGidsToSend[0] = 0;
      } else {
        numGidsToSend[0] = it->second.size();
      }
      sizeReq = isend<int, int_type>(numGidsToSend, p, sizeTag, *comm);
      wait<int>(*comm, outArg(sizeReq));
    }
  } else {
    // Wait on the receive-size message to finish.
    wait<int>(*comm, outArg(sizeReq));
  }
 
  // Allocate / get the array for my GIDs.
  // Only Process 0 will have its actual GIDs at this point.
  ArrayRCP<GO> myGids;
  int_type myNumGids = 0;
  if (myRank == 0) {
    GO* myGidsRaw = NULL;
 
    typename std::map<int, Array<GO>>::iterator it = pid2gids.find(0);
    if (it != pid2gids.end()) {
      myGidsRaw = it->second.getRawPtr();
      myNumGids = it->second.size();
      // Nonowning ArrayRCP just views the Array.
      myGids = arcp<GO>(myGidsRaw, 0, myNumGids, false);
    }
  } else {  // myRank != 0
    myNumGids = numGidsToRecv[0];
    myGids    = arcp<GO>(myNumGids);
  }
 
  if (myRank != 0) {
    // Post receive for data, now that we know how much data we
    // will receive.  Only post receive if my process actually
    // has nonzero GIDs.
    if (myNumGids > 0) {
      dataReq = ireceive<int, GO>(myGids, 0, dataTag, *comm);
    }
  }
 
  for (int p = 1; p < numProcs; ++p) {
    if (myRank == 0) {
      ArrayRCP<GO> sendGids;  // to send to Process p
      GO* sendGidsRaw      = NULL;
      int_type numSendGids = 0;
 
      typename std::map<int, Array<GO>>::iterator it = pid2gids.find(p);
      if (it != pid2gids.end()) {
        numSendGids = it->second.size();
        sendGidsRaw = it->second.getRawPtr();
        sendGids    = arcp<GO>(sendGidsRaw, 0, numSendGids, false);
      }
      // Only send if that process actually has nonzero GIDs.
      if (numSendGids > 0) {
        dataReq = isend<int, GO>(sendGids, p, dataTag, *comm);
      }
      wait<int>(*comm, outArg(dataReq));
    } else if (myRank == p) {
      // Wait on my receive of GIDs to finish.
      wait<int>(*comm, outArg(dataReq));
    }
  }  // for each process rank p in 1, 2, ..., numProcs-1
 
  if (debug) {
    std::ostringstream os;
    os << myRank << ": readMap: creating Map" << endl;
    *err << os.str();
  }
  const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid();
  RCP<const map_type> newMap;
 
  // Create the Map; test whether the constructor threw.  This
  // avoids deadlock and makes error reporting more readable.
 
  int lclSuccess = 1;
  int gblSuccess = 0;  // output argument
  std::ostringstream errStrm;
  try {
    newMap = rcp(new map_type(INVALID, myGids(), indexBase, comm));
  } catch (std::exception& e) {
    lclSuccess = 0;
    errStrm << "Process " << comm->getRank() << " threw an exception: "
            << e.what() << std::endl;
  } catch (...) {
    lclSuccess = 0;
    errStrm << "Process " << comm->getRank() << " threw an exception "
                                                "not a subclass of std::exception"
            << std::endl;
  }
  Teuchos::reduceAll<int, int>(*comm, Teuchos::REDUCE_MIN,
                               lclSuccess, Teuchos::outArg(gblSuccess));
  if (gblSuccess != 1) {
    Tpetra::Details::gathervPrint(std::cerr, errStrm.str(), *comm);
  }
  TEUCHOS_TEST_FOR_EXCEPTION(gblSuccess != 1, std::runtime_error, "Map constructor failed!");
 
  if (!err.is_null()) {
    err->popTab();
  }
  if (debug) {
    std::ostringstream os;
    os << myRank << ": readMap: done" << endl;
    *err << os.str();
  }
  if (!err.is_null()) {
    err->popTab();
  }
  return newMap;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
int MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::encodeDataType(const std::string& dataType) {
  if (dataType == "real" || dataType == "integer") {
    return 0;
  } else if (dataType == "complex") {
    return 1;
  } else if (dataType == "pattern") {
    return 2;
  } else {
    // We should never get here, since Banner validates the
    // reported data type and ensures it is one of the accepted
    // values.
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
                               "Unrecognized Matrix Market data type \"" << dataType
                                                                         << "\".  We should never get here.  "
                                                                            "Please report this bug to the Tpetra developers.");
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type>
MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::readSparsePerRank(const std::string& filename_prefix,
                                                                                 const std::string& filename_suffix,
                                                                                 const Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& rowMap,
                                                                                 Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& colMap,
                                                                                 const Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& domainMap,
                                                                                 const Teuchos::RCP<const typename MatrixMarketReader<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type>& rangeMap,
                                                                                 const bool callFillComplete,
                                                                                 const bool tolerant,
                                                                                 const int ranksToReadAtOnce,
                                                                                 const bool debug) {
  using ST  = scalar_type;
  using LO  = local_ordinal_type;
  using GO  = global_ordinal_type;
  using STS = typename Teuchos::ScalarTraits<ST>;
  using Teuchos::arcp;
  using Teuchos::ArrayRCP;
  using Teuchos::RCP;
  using Teuchos::rcp;
 
  // Sanity Checks
  // Fast checks for invalid input.  We can't check other
  // attributes of the Maps until we've read in the matrix
  // dimensions.
  TEUCHOS_TEST_FOR_EXCEPTION(
      rowMap.is_null(), std::invalid_argument,
      "Row Map must be nonnull.");
  Teuchos::RCP<const Teuchos::Comm<int>> comm = rowMap->getComm();
  TEUCHOS_TEST_FOR_EXCEPTION(comm.is_null(), std::invalid_argument,
                             "The input row map's communicator (Teuchos::Comm object) is null.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      rangeMap.is_null(), std::invalid_argument,
      "Range Map must be nonnull.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      domainMap.is_null(), std::invalid_argument,
      "Domain Map must be nonnull.");
  TEUCHOS_TEST_FOR_EXCEPTION(
      domainMap->getComm().getRawPtr() != comm.getRawPtr(),
      std::invalid_argument,
      "The specified domain Map's communicator (domainMap->getComm())"
      "differs from row Map's communicator");
  TEUCHOS_TEST_FOR_EXCEPTION(
      rangeMap->getComm().getRawPtr() != comm.getRawPtr(),
      std::invalid_argument,
      "The specified range Map's communicator (rangeMap->getComm())"
      "differs from row Map's communicator");
 
  // Setup
  const int myRank     = comm->getRank();
  const int numProc    = comm->getSize();
  std::string filename = filename_prefix + std::to_string(myRank) + filename_suffix;
 
  // Bounds check the writing limits
  int rank_limit = std::min(std::max(ranksToReadAtOnce, 1), numProc);
 
  // Data structures for constructor
  ArrayRCP<size_t> numEntriesPerRow;
  ArrayRCP<size_t> rowPtr;
  ArrayRCP<global_ordinal_type> colInd;
  ArrayRCP<scalar_type> values;
  std::ostringstream errMsg;

  // Start the reading of the banners to get
  // local row / nnz counts and then read the
  // data.  We'll pack everything into big ol'
  // rowptr/colind/values ArrayRCPs
  bool success        = true;
  int bannerIsCorrect = 1, readSuccess = 1;
  LO numRows, numCols, numNonzeros;
  for (int base_rank = 0; base_rank < numProc; base_rank += rank_limit) {
    int stop = std::min(base_rank + rank_limit, numProc);
 
    // Is my rank in this batch?
    if (base_rank <= myRank && myRank < stop) {
      // My turn to read
      std::ifstream in(filename);
      using Teuchos::MatrixMarket::Banner;
      size_t lineNumber = 1;
      RCP<const Banner> pBanner;
      try {
        pBanner = readBanner(in, lineNumber, tolerant, debug);
      } catch (std::exception& e) {
        errMsg << "Attempt to read the Matrix Market file's Banner line "
                  "threw an exception: "
               << e.what();
        bannerIsCorrect = 0;
      }
      if (bannerIsCorrect) {
        // Validate the Banner for the case of a sparse matrix.
        // We validate on Proc 0, since it reads the Banner.
 
        // In intolerant mode, the matrix type must be "coordinate".
        if (!tolerant && pBanner->matrixType() != "coordinate") {
          bannerIsCorrect = 0;
          errMsg << "The Matrix Market input file must contain a "
                    "\"coordinate\"-format sparse matrix in order to create a "
                    "Tpetra::CrsMatrix object from it, but the file's matrix "
                    "type is \""
                 << pBanner->matrixType() << "\" instead.";
        }
        // In tolerant mode, we allow the matrix type to be
        // anything other than "array" (which would mean that
        // the file contains a dense matrix).
        if (tolerant && pBanner->matrixType() == "array") {
          bannerIsCorrect = 0;
          errMsg << "Matrix Market file must contain a \"coordinate\"-"
                    "format sparse matrix in order to create a Tpetra::CrsMatrix "
                    "object from it, but the file's matrix type is \"array\" "
                    "instead.  That probably means the file contains dense matrix "
                    "data.";
        }
      }
 
      // Unpacked coordinate matrix dimensions
      using Teuchos::MatrixMarket::readCoordinateDimensions;
      success = readCoordinateDimensions(in, numRows, numCols,
                                         numNonzeros, lineNumber,
                                         tolerant);
 
      // Sanity checking of headers
      TEUCHOS_TEST_FOR_EXCEPTION(numRows != (LO)rowMap->getLocalNumElements(), std::invalid_argument,
                                 "# rows in file does not match rowmap.");
      TEUCHOS_TEST_FOR_EXCEPTION(!colMap.is_null() && numCols != (LO)colMap->getLocalNumElements(), std::invalid_argument,
                                 "# rows in file does not match colmap.");
 
      // Read the data
      typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> raw_adder_type;
      bool tolerant_required = true;
      Teuchos::RCP<raw_adder_type> pRaw =
          Teuchos::rcp(new raw_adder_type(numRows, numCols, numNonzeros, tolerant_required, debug));
      RCP<adder_type> pAdder = Teuchos::rcp(new adder_type(pRaw, pBanner->symmType()));
 
      if (debug) {
        std::cerr << "-- Reading matrix data" << std::endl;
      }
 
      try {
        // Reader for "coordinate" format sparse matrix data.
        typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
                                                       global_ordinal_type, scalar_type, STS::isComplex>
            reader_type;
        reader_type reader(pAdder);
 
        // Read the sparse matrix entries.
        std::pair<bool, std::vector<size_t>> results = reader.read(in, lineNumber, tolerant_required, debug);
 
        readSuccess = results.first ? 1 : 0;
      } catch (std::exception& e) {
        readSuccess = 0;
        errMsg << e.what();
      }

      // Create the CSR Arrays
      typedef Teuchos::MatrixMarket::Raw::Element<scalar_type, global_ordinal_type> element_type;
 
      // Additively merge duplicate matrix entries.
      pAdder->getAdder()->merge();
 
      // Get a temporary const view of the merged matrix entries.
      const std::vector<element_type>& entries = pAdder->getAdder()->getEntries();
 
      // Number of matrix entries (after merging).
      const size_t numEntries = (size_t)entries.size();
 
      if (debug) {
        std::cerr << "----- Proc " << myRank << ": Matrix has numRows=" << numRows
                  << " rows and numEntries=" << numEntries
                  << " entries." << std::endl;
      }
 
      // Make space for the CSR matrix data.  Converting to
      // CSR is easier if we fill numEntriesPerRow with zeros
      // at first.
      numEntriesPerRow = arcp<size_t>(numRows);
      std::fill(numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
      rowPtr = arcp<size_t>(numRows + 1);
      std::fill(rowPtr.begin(), rowPtr.end(), 0);
      colInd = arcp<global_ordinal_type>(numEntries);
      values = arcp<scalar_type>(numEntries);
 
      // Convert from array-of-structs coordinate format to CSR
      // (compressed sparse row) format.
      global_ordinal_type l_prvRow = 0;
      size_t curPos                = 0;
      LO INVALID                   = Teuchos::OrdinalTraits<LO>::invalid();
      rowPtr[0]                    = 0;
      LO indexBase                 = rowMap->getIndexBase();
      for (curPos = 0; curPos < numEntries; ++curPos) {
        const element_type& curEntry     = entries[curPos];
        const global_ordinal_type curRow = curEntry.rowIndex() + indexBase;
        LO l_curRow                      = rowMap->getLocalElement(curRow);
 
        TEUCHOS_TEST_FOR_EXCEPTION(l_curRow == INVALID, std::logic_error,
                                   "Current global row " << curRow << " is invalid.");
 
        TEUCHOS_TEST_FOR_EXCEPTION(l_curRow < l_prvRow, std::logic_error,
                                   "Row indices are out of order, even though they are supposed "
                                   "to be sorted.  curRow = "
                                       << l_curRow << ", prvRow = "
                                       << l_prvRow << ", at curPos = " << curPos << ".  Please report "
                                                                                    "this bug to the Tpetra developers.");
        if (l_curRow > l_prvRow) {
          for (LO r = l_prvRow + 1; r <= l_curRow; ++r) {
            rowPtr[r] = curPos;
          }
          l_prvRow = l_curRow;
        }
        numEntriesPerRow[l_curRow]++;
        colInd[curPos] = curEntry.colIndex() + indexBase;
        values[curPos] = curEntry.value();
      }
      // rowPtr has one more entry than numEntriesPerRow.  The
      // last entry of rowPtr is the number of entries in
      // colInd and values.
      rowPtr[numRows] = numEntries;
 
    }  // end base_rank <= myRank < stop
 
    // Barrier between batches to keep the filesystem happy
    comm->barrier();
 
  }  // end outer rank loop
 
  // Call the matrix constructor and fill.  This isn't particularly efficient
  RCP<sparse_matrix_type> A;
  if (colMap.is_null()) {
    A = rcp(new sparse_matrix_type(rowMap, numEntriesPerRow()));
    for (size_t i = 0; i < rowMap->getLocalNumElements(); i++) {
      GO g_row     = rowMap->getGlobalElement(i);
      size_t start = rowPtr[i];
      size_t size  = rowPtr[i + 1] - rowPtr[i];
      if (size > 0) {
        A->insertGlobalValues(g_row, size, &values[start], &colInd[start]);
      }
    }
  } else {
    throw std::runtime_error("Reading with a column map is not yet implemented");
  }
  RCP<const map_type> myDomainMap = domainMap.is_null() ? rowMap : domainMap;
  RCP<const map_type> myRangeMap  = rangeMap.is_null() ? rowMap : rangeMap;
 
  A->fillComplete(myDomainMap, myRangeMap);
 
  if (!readSuccess)
    success = false;
  TEUCHOS_TEST_FOR_EXCEPTION(success == false, std::runtime_error,
                             "Read failed.");
 
  return A;
}  // end readSparsePerRank
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparseFile(const std::string& filename,
                                                                                    const sparse_matrix_type& matrix,
                                                                                    const std::string& matrixName,
                                                                                    const std::string& matrixDescription,
                                                                                    const bool debug) {
  trcp_tcomm_t comm = matrix.getComm();
  TEUCHOS_TEST_FOR_EXCEPTION(comm.is_null(), std::invalid_argument,
                             "The input matrix's communicator (Teuchos::Comm object) is null.");
  const int myRank = comm->getRank();
 
  auto out = MatrixMarketWriter::openOutFileOnRankZero(comm, filename, myRank, true);
 
  writeSparse(out, matrix, matrixName, matrixDescription, debug);
  // We can rely on the destructor of the output stream to close
  // the file on scope exit, even if writeSparse() throws an
  // exception.
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparseFile(const std::string& filename,
                                                                                    const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
                                                                                    const std::string& matrixName,
                                                                                    const std::string& matrixDescription,
                                                                                    const bool debug) {
  TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::invalid_argument,
                             "The input matrix is null.");
  writeSparseFile(filename, *pMatrix, matrixName,
                  matrixDescription, debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparseFile(const std::string& filename,
                                                                                    const sparse_matrix_type& matrix,
                                                                                    const bool debug) {
  writeSparseFile(filename, matrix, "", "", debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparseFile(const std::string& filename,
                                                                                    const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
                                                                                    const bool debug) {
  writeSparseFile(filename, *pMatrix, "", "", debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparse(std::ostream& out,
                                                                                const sparse_matrix_type& matrix,
                                                                                const std::string& matrixName,
                                                                                const std::string& matrixDescription,
                                                                                const bool debug) {
  using std::cerr;
  using std::endl;
  using Teuchos::ArrayView;
  using Teuchos::Comm;
  using Teuchos::FancyOStream;
  using Teuchos::getFancyOStream;
  using Teuchos::null;
  using Teuchos::RCP;
  using Teuchos::rcpFromRef;
  using ST  = scalar_type;
  using LO  = local_ordinal_type;
  using GO  = global_ordinal_type;
  using STS = typename Teuchos::ScalarTraits<ST>;
 
  // Make the output stream write floating-point numbers in
  // scientific notation.  It will politely put the output
  // stream back to its state on input, when this scope
  // terminates.
  Teuchos::SetScientific<ST> sci(out);
 
  // Get the matrix's communicator.
  trcp_tcomm_t comm = matrix.getComm();
  TEUCHOS_TEST_FOR_EXCEPTION(
      comm.is_null(), std::invalid_argument,
      "The input matrix's communicator (Teuchos::Comm object) is null.");
  const int myRank = comm->getRank();
 
  // Optionally, make a stream for debugging output.
  RCP<FancyOStream> err =
      debug ? getFancyOStream(rcpFromRef(std::cerr)) : null;
  if (debug) {
    std::ostringstream os;
    os << myRank << ": writeSparse" << endl;
    *err << os.str();
    comm->barrier();
    os << "-- " << myRank << ": past barrier" << endl;
    *err << os.str();
  }
 
  // Whether to print debugging output to stderr.
  const bool debugPrint = debug && myRank == 0;
 
  RCP<const map_type> rowMap    = matrix.getRowMap();
  RCP<const map_type> colMap    = matrix.getColMap();
  RCP<const map_type> domainMap = matrix.getDomainMap();
  RCP<const map_type> rangeMap  = matrix.getRangeMap();
 
  if (!rowMap->haveGlobalConstants())
    Teuchos::rcp_const_cast<map_type>(rowMap)->computeGlobalConstants();
  if (!colMap->haveGlobalConstants())
    Teuchos::rcp_const_cast<map_type>(colMap)->computeGlobalConstants();
  if (!domainMap->haveGlobalConstants())
    Teuchos::rcp_const_cast<map_type>(domainMap)->computeGlobalConstants();
  if (!rangeMap->haveGlobalConstants())
    Teuchos::rcp_const_cast<map_type>(rangeMap)->computeGlobalConstants();
 
  const global_size_t numRows = rangeMap->getMaxAllGlobalIndex() + 1 - rangeMap->getIndexBase();
  const global_size_t numCols = domainMap->getMaxAllGlobalIndex() + 1 - domainMap->getIndexBase();
 
  if (debug && myRank == 0) {
    std::ostringstream os;
    os << "-- Input sparse matrix is:"
       << "---- " << numRows << " x " << numCols << endl
       << "---- "
       << (matrix.isGloballyIndexed() ? "Globally" : "Locally")
       << " indexed." << endl
       << "---- Its row map has " << rowMap->getGlobalNumElements()
       << " elements." << endl
       << "---- Its col map has " << colMap->getGlobalNumElements()
       << " elements." << endl;
    *err << os.str();
  }
  // Make the "gather" row map, where Proc 0 owns all rows and
  // the other procs own no rows.
  const size_t localNumRows = (myRank == 0) ? numRows : 0;
  if (debug) {
    std::ostringstream os;
    os << "-- " << myRank << ": making gatherRowMap" << endl;
    *err << os.str();
  }
  RCP<const map_type> gatherRowMap =
      Details::computeGatherMap(rowMap, err, debug);
 
  // Since the matrix may in general be non-square, we need to
  // make a column map as well.  In this case, the column map
  // contains all the columns of the original matrix, because we
  // are gathering the whole matrix onto Proc 0.  We call
  // computeGatherMap to preserve the original order of column
  // indices over all the processes.
  const size_t localNumCols = (myRank == 0) ? numCols : 0;
  RCP<const map_type> gatherColMap =
      Details::computeGatherMap(colMap, err, debug);
 
  // Current map is the source map, gather map is the target map.
  typedef Import<LO, GO, node_type> import_type;
  import_type importer(rowMap, gatherRowMap);
 
  // Create a new CrsMatrix to hold the result of the import.
  // The constructor needs a column map as well as a row map,
  // for the case that the matrix is not square.
  RCP<sparse_matrix_type> newMatrix =
      rcp(new sparse_matrix_type(gatherRowMap, gatherColMap,
                                 static_cast<size_t>(0)));
  // Import the sparse matrix onto Proc 0.
  newMatrix->doImport(matrix, importer, INSERT);
 
  // fillComplete() needs the domain and range maps for the case
  // that the matrix is not square.
  {
    RCP<const map_type> gatherDomainMap =
        rcp(new map_type(numCols, localNumCols,
                         domainMap->getIndexBase(),
                         comm));
    RCP<const map_type> gatherRangeMap =
        rcp(new map_type(numRows, localNumRows,
                         rangeMap->getIndexBase(),
                         comm));
    newMatrix->fillComplete(gatherDomainMap, gatherRangeMap);
  }
 
  if (debugPrint) {
    cerr << "-- Output sparse matrix is:"
         << "---- " << newMatrix->getRangeMap()->getGlobalNumElements()
         << " x "
         << newMatrix->getDomainMap()->getGlobalNumElements()
         << " with "
         << newMatrix->getGlobalNumEntries() << " entries;" << endl
         << "---- "
         << (newMatrix->isGloballyIndexed() ? "Globally" : "Locally")
         << " indexed." << endl
         << "---- Its row map has "
         << newMatrix->getRowMap()->getGlobalNumElements()
         << " elements, with index base "
         << newMatrix->getRowMap()->getIndexBase() << "." << endl
         << "---- Its col map has "
         << newMatrix->getColMap()->getGlobalNumElements()
         << " elements, with index base "
         << newMatrix->getColMap()->getIndexBase() << "." << endl
         << "---- Element count of output matrix's column Map may differ "
         << "from that of the input matrix's column Map, if some columns "
         << "of the matrix contain no entries." << endl;
  }
 
  //
  // Print the metadata and the matrix entries on Rank 0.
  //
  if (myRank == 0) {
    // Print the Matrix Market banner line.  CrsMatrix stores
    // data nonsymmetrically ("general").  This implies that
    // readSparse() on a symmetrically stored input file,
    // followed by writeSparse() on the resulting sparse matrix,
    // will result in an output file with a different banner
    // line than the original input file.
    out << "%%MatrixMarket matrix coordinate "
        << (STS::isComplex ? "complex" : "real")
        << " general" << endl;
 
    // Print comments (the matrix name and / or description).
    if (matrixName != "") {
      printAsComment(out, matrixName);
    }
    if (matrixDescription != "") {
      printAsComment(out, matrixDescription);
    }
 
    // Print the Matrix Market header (# rows, # columns, #
    // nonzeros).  Use the range resp. domain map for the number
    // of rows resp. columns, since Tpetra::CrsMatrix uses the
    // column map for the number of columns.  That only
    // corresponds to the "linear-algebraic" number of columns
    // when the column map is uniquely owned (a.k.a. one-to-one),
    // which only happens if the matrix is (block) diagonal.
    out << newMatrix->getRangeMap()->getMaxAllGlobalIndex() + 1 - newMatrix->getRangeMap()->getIndexBase() << " "
        << newMatrix->getDomainMap()->getMaxAllGlobalIndex() + 1 - newMatrix->getDomainMap()->getIndexBase() << " "
        << newMatrix->getGlobalNumEntries() << endl;
 
    // The Matrix Market format expects one-based row and column
    // indices.  We'll convert the indices on output from
    // whatever index base they use to one-based indices.
    const GO rowIndexBase = gatherRowMap->getIndexBase();
    const GO colIndexBase = newMatrix->getColMap()->getIndexBase();
    //
    // Print the entries of the matrix.
    //
    // newMatrix can never be globally indexed, since we called
    // fillComplete() on it.  We include code for both cases
    // (globally or locally indexed) just in case that ever
    // changes.
    if (newMatrix->isGloballyIndexed()) {
      // We know that the "gather" row Map is contiguous, so we
      // don't need to get the list of GIDs.
      const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex();
      const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex();
      for (GO globalRowIndex = minAllGlobalIndex;
           globalRowIndex <= maxAllGlobalIndex;  // inclusive range
           ++globalRowIndex) {
        typename sparse_matrix_type::global_inds_host_view_type ind;
        typename sparse_matrix_type::values_host_view_type val;
        newMatrix->getGlobalRowView(globalRowIndex, ind, val);
        for (size_t ii = 0; ii < ind.extent(0); ii++) {
          const GO globalColIndex = ind(ii);
          // Convert row and column indices to 1-based.
          // This works because the global index type is signed.
          out << (globalRowIndex + 1 - rowIndexBase) << " "
              << (globalColIndex + 1 - colIndexBase) << " ";
          if (STS::isComplex) {
            out << STS::real(val(ii)) << " " << STS::imag(val(ii));
          } else {
            out << val(ii);
          }
          out << endl;
        }     // For each entry in the current row
      }       // For each row of the "gather" matrix
    } else {  // newMatrix is locally indexed
      using OTG = Teuchos::OrdinalTraits<GO>;
      for (LO localRowIndex = gatherRowMap->getMinLocalIndex();
           localRowIndex <= gatherRowMap->getMaxLocalIndex();
           ++localRowIndex) {
        // Convert from local to global row index.
        const GO globalRowIndex =
            gatherRowMap->getGlobalElement(localRowIndex);
        TEUCHOS_TEST_FOR_EXCEPTION(
            globalRowIndex == OTG::invalid(), std::logic_error,
            "Failed to convert the supposed local row index "
                << localRowIndex << " into a global row index.  "
                                    "Please report this bug to the Tpetra developers.");
        typename sparse_matrix_type::local_inds_host_view_type ind;
        typename sparse_matrix_type::values_host_view_type val;
        newMatrix->getLocalRowView(localRowIndex, ind, val);
        for (size_t ii = 0; ii < ind.extent(0); ii++) {
          // Convert the column index from local to global.
          const GO globalColIndex =
              newMatrix->getColMap()->getGlobalElement(ind(ii));
          TEUCHOS_TEST_FOR_EXCEPTION(
              globalColIndex == OTG::invalid(), std::logic_error,
              "On local row " << localRowIndex << " of the sparse matrix: "
                                                  "Failed to convert the supposed local column index "
                              << ind(ii) << " into a global column index.  Please report "
                                            "this bug to the Tpetra developers.");
          // Convert row and column indices to 1-based.
          // This works because the global index type is signed.
          out << (globalRowIndex + 1 - rowIndexBase) << " "
              << (globalColIndex + 1 - colIndexBase) << " ";
          if (STS::isComplex) {
            out << STS::real(val(ii)) << " " << STS::imag(val(ii));
          } else {
            out << val(ii);
          }
          out << endl;
        }  // For each entry in the current row
      }    // For each row of the "gather" matrix
    }      // Whether the "gather" matrix is locally or globally indexed
  }        // If my process' rank is 0
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparse(std::ostream& out,
                                                                                const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
                                                                                const std::string& matrixName,
                                                                                const std::string& matrixDescription,
                                                                                const bool debug) {
  TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::invalid_argument,
                             "The input matrix is null.");
  writeSparse(out, *pMatrix, matrixName, matrixDescription, debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparseGraph(std::ostream& out,
                                                                                     const crs_graph_type& graph,
                                                                                     const std::string& graphName,
                                                                                     const std::string& graphDescription,
                                                                                     const bool debug) {
  using std::cerr;
  using std::endl;
  using Teuchos::ArrayView;
  using Teuchos::Comm;
  using Teuchos::FancyOStream;
  using Teuchos::getFancyOStream;
  using Teuchos::null;
  using Teuchos::RCP;
  using Teuchos::rcpFromRef;
  typedef local_ordinal_type LO;
  typedef global_ordinal_type GO;
 
  // Get the graph's communicator.  Processes on which the
  // graph's Map or communicator is null don't participate in
  // this operation.  This function shouldn't even be called on
  // those processes.
  auto rowMap = graph.getRowMap();
  if (rowMap.is_null()) {
    return;
  }
  auto comm = rowMap->getComm();
  if (comm.is_null()) {
    return;
  }
  const int myRank = comm->getRank();
 
  // Optionally, make a stream for debugging output.
  RCP<FancyOStream> err =
      debug ? getFancyOStream(rcpFromRef(std::cerr)) : null;
  if (debug) {
    std::ostringstream os;
    os << myRank << ": writeSparseGraph" << endl;
    *err << os.str();
    comm->barrier();
    os << "-- " << myRank << ": past barrier" << endl;
    *err << os.str();
  }
 
  // Whether to print debugging output to stderr.
  const bool debugPrint = debug && myRank == 0;
 
  // We've already gotten the rowMap above.
  auto colMap    = graph.getColMap();
  auto domainMap = graph.getDomainMap();
  auto rangeMap  = graph.getRangeMap();
 
  if (!rowMap->haveGlobalConstants())
    Teuchos::rcp_const_cast<map_type>(rowMap)->computeGlobalConstants();
  if (!colMap->haveGlobalConstants())
    Teuchos::rcp_const_cast<map_type>(colMap)->computeGlobalConstants();
  if (!domainMap->haveGlobalConstants())
    Teuchos::rcp_const_cast<map_type>(domainMap)->computeGlobalConstants();
  if (!rangeMap->haveGlobalConstants())
    Teuchos::rcp_const_cast<map_type>(rangeMap)->computeGlobalConstants();
 
  const global_size_t numRows = rangeMap->getGlobalNumElements();
  const global_size_t numCols = domainMap->getGlobalNumElements();
 
  if (debug && myRank == 0) {
    std::ostringstream os;
    os << "-- Input sparse graph is:"
       << "---- " << numRows << " x " << numCols << " with "
       << graph.getGlobalNumEntries() << " entries;" << endl
       << "---- "
       << (graph.isGloballyIndexed() ? "Globally" : "Locally")
       << " indexed." << endl
       << "---- Its row Map has " << rowMap->getGlobalNumElements()
       << " elements." << endl
       << "---- Its col Map has " << colMap->getGlobalNumElements()
       << " elements." << endl;
    *err << os.str();
  }
  // Make the "gather" row map, where Proc 0 owns all rows and
  // the other procs own no rows.
  const size_t localNumRows = (myRank == 0) ? numRows : 0;
  if (debug) {
    std::ostringstream os;
    os << "-- " << myRank << ": making gatherRowMap" << endl;
    *err << os.str();
  }
  auto gatherRowMap = Details::computeGatherMap(rowMap, err, debug);
 
  // Since the graph may in general be non-square, we need to
  // make a column map as well.  In this case, the column map
  // contains all the columns of the original graph, because we
  // are gathering the whole graph onto Proc 0.  We call
  // computeGatherMap to preserve the original order of column
  // indices over all the processes.
  const size_t localNumCols = (myRank == 0) ? numCols : 0;
  auto gatherColMap         = Details::computeGatherMap(colMap, err, debug);
 
  // Current map is the source map, gather map is the target map.
  Import<LO, GO, node_type> importer(rowMap, gatherRowMap);
 
  // Create a new CrsGraph to hold the result of the import.
  // The constructor needs a column map as well as a row map,
  // for the case that the graph is not square.
  crs_graph_type newGraph(gatherRowMap, gatherColMap,
                          static_cast<size_t>(0));
  // Import the sparse graph onto Proc 0.
  newGraph.doImport(graph, importer, INSERT);
 
  // fillComplete() needs the domain and range maps for the case
  // that the graph is not square.
  {
    RCP<const map_type> gatherDomainMap =
        rcp(new map_type(numCols, localNumCols,
                         domainMap->getIndexBase(),
                         comm));
    RCP<const map_type> gatherRangeMap =
        rcp(new map_type(numRows, localNumRows,
                         rangeMap->getIndexBase(),
                         comm));
    newGraph.fillComplete(gatherDomainMap, gatherRangeMap);
  }
 
  if (debugPrint) {
    cerr << "-- Output sparse graph is:"
         << "---- " << newGraph.getRangeMap()->getGlobalNumElements()
         << " x "
         << newGraph.getDomainMap()->getGlobalNumElements()
         << " with "
         << newGraph.getGlobalNumEntries() << " entries;" << endl
         << "---- "
         << (newGraph.isGloballyIndexed() ? "Globally" : "Locally")
         << " indexed." << endl
         << "---- Its row map has "
         << newGraph.getRowMap()->getGlobalNumElements()
         << " elements, with index base "
         << newGraph.getRowMap()->getIndexBase() << "." << endl
         << "---- Its col map has "
         << newGraph.getColMap()->getGlobalNumElements()
         << " elements, with index base "
         << newGraph.getColMap()->getIndexBase() << "." << endl
         << "---- Element count of output graph's column Map may differ "
         << "from that of the input matrix's column Map, if some columns "
         << "of the matrix contain no entries." << endl;
  }
 
  //
  // Print the metadata and the graph entries on Process 0 of
  // the graph's communicator.
  //
  if (myRank == 0) {
    // Print the Matrix Market banner line.  CrsGraph stores
    // data nonsymmetrically ("general").  This implies that
    // readSparseGraph() on a symmetrically stored input file,
    // followed by writeSparseGraph() on the resulting sparse
    // graph, will result in an output file with a different
    // banner line than the original input file.
    out << "%%MatrixMarket matrix coordinate pattern general" << endl;
 
    // Print comments (the graph name and / or description).
    if (graphName != "") {
      printAsComment(out, graphName);
    }
    if (graphDescription != "") {
      printAsComment(out, graphDescription);
    }
 
    // Print the Matrix Market header (# rows, # columns, #
    // stored entries).  Use the range resp. domain map for the
    // number of rows resp. columns, since Tpetra::CrsGraph uses
    // the column map for the number of columns.  That only
    // corresponds to the "linear-algebraic" number of columns
    // when the column map is uniquely owned
    // (a.k.a. one-to-one), which only happens if the graph is
    // block diagonal (one block per process).
    out << newGraph.getRangeMap()->getMaxAllGlobalIndex() + 1 - newGraph.getRangeMap()->getIndexBase() << " "
        << newGraph.getDomainMap()->getMaxAllGlobalIndex() + 1 - newGraph.getDomainMap()->getIndexBase() << " "
        << newGraph.getGlobalNumEntries() << endl;
 
    // The Matrix Market format expects one-based row and column
    // indices.  We'll convert the indices on output from
    // whatever index base they use to one-based indices.
    const GO rowIndexBase = gatherRowMap->getIndexBase();
    const GO colIndexBase = newGraph.getColMap()->getIndexBase();
    //
    // Print the entries of the graph.
    //
    // newGraph can never be globally indexed, since we called
    // fillComplete() on it.  We include code for both cases
    // (globally or locally indexed) just in case that ever
    // changes.
    if (newGraph.isGloballyIndexed()) {
      // We know that the "gather" row Map is contiguous, so we
      // don't need to get the list of GIDs.
      const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex();
      const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex();
      for (GO globalRowIndex = minAllGlobalIndex;
           globalRowIndex <= maxAllGlobalIndex;  // inclusive range
           ++globalRowIndex) {
        typename crs_graph_type::global_inds_host_view_type ind;
        newGraph.getGlobalRowView(globalRowIndex, ind);
        for (size_t ii = 0; ii < ind.extent(0); ii++) {
          const GO globalColIndex = ind(ii);
          // Convert row and column indices to 1-based.
          // This works because the global index type is signed.
          out << (globalRowIndex + 1 - rowIndexBase) << " "
              << (globalColIndex + 1 - colIndexBase) << " ";
          out << endl;
        }     // For each entry in the current row
      }       // For each row of the "gather" graph
    } else {  // newGraph is locally indexed
      typedef Teuchos::OrdinalTraits<GO> OTG;
      for (LO localRowIndex = gatherRowMap->getMinLocalIndex();
           localRowIndex <= gatherRowMap->getMaxLocalIndex();
           ++localRowIndex) {
        // Convert from local to global row index.
        const GO globalRowIndex =
            gatherRowMap->getGlobalElement(localRowIndex);
        TEUCHOS_TEST_FOR_EXCEPTION(globalRowIndex == OTG::invalid(), std::logic_error,
                                   "Failed "
                                   "to convert the supposed local row index "
                                       << localRowIndex << " into a global row index.  Please report this bug to the "
                                                           "Tpetra developers.");
        typename crs_graph_type::local_inds_host_view_type ind;
        newGraph.getLocalRowView(localRowIndex, ind);
        for (size_t ii = 0; ii < ind.extent(0); ii++) {
          // Convert the column index from local to global.
          const GO globalColIndex =
              newGraph.getColMap()->getGlobalElement(ind(ii));
          TEUCHOS_TEST_FOR_EXCEPTION(
              globalColIndex == OTG::invalid(), std::logic_error,
              "On local row " << localRowIndex << " of the sparse graph: "
                                                  "Failed to convert the supposed local column index "
                              << ind(ii) << " into a global column index.  Please report "
                                            "this bug to the Tpetra developers.");
          // Convert row and column indices to 1-based.
          // This works because the global index type is signed.
          out << (globalRowIndex + 1 - rowIndexBase) << " "
              << (globalColIndex + 1 - colIndexBase) << " ";
          out << endl;
        }  // For each entry in the current row
      }    // For each row of the "gather" graph
    }      // Whether the "gather" graph is locally or globally indexed
  }        // If my process' rank is 0
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparseGraph(std::ostream& out,
                                                                                     const crs_graph_type& graph,
                                                                                     const bool debug) {
  writeSparseGraph(out, graph, "", "", debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparseGraphFile(const std::string& filename,
                                                                                         const crs_graph_type& graph,
                                                                                         const std::string& graphName,
                                                                                         const std::string& graphDescription,
                                                                                         const bool debug) {
  auto comm = graph.getComm();
  if (comm.is_null()) {
    // Processes on which the communicator is null shouldn't
    // even call this function.  The convention is that
    // processes on which the object's communicator is null do
    // not participate in collective operations involving the
    // object.
    return;
  }
  const int myRank = comm->getRank();
 
  auto out = MatrixMarketWriter::openOutFileOnRankZero(comm, filename, myRank, true);
 
  writeSparseGraph(out, graph, graphName, graphDescription, debug);
  // We can rely on the destructor of the output stream to close
  // the file on scope exit, even if writeSparseGraph() throws
  // an exception.
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparseGraphFile(const std::string& filename,
                                                                                         const crs_graph_type& graph,
                                                                                         const bool debug) {
  writeSparseGraphFile(filename, graph, "", "", debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparseGraphFile(const std::string& filename,
                                                                                         const Teuchos::RCP<const crs_graph_type>& pGraph,
                                                                                         const std::string& graphName,
                                                                                         const std::string& graphDescription,
                                                                                         const bool debug) {
  writeSparseGraphFile(filename, *pGraph, graphName, graphDescription, debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparseGraphFile(const std::string& filename,
                                                                                         const Teuchos::RCP<const crs_graph_type>& pGraph,
                                                                                         const bool debug) {
  writeSparseGraphFile(filename, *pGraph, "", "", debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparse(std::ostream& out,
                                                                                const sparse_matrix_type& matrix,
                                                                                const bool debug) {
  writeSparse(out, matrix, "", "", debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparse(std::ostream& out,
                                                                                const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
                                                                                const bool debug) {
  writeSparse(out, *pMatrix, "", "", debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeDenseFile(const std::string& filename,
                                                                                   const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type& X,
                                                                                   const std::string& matrixName,
                                                                                   const std::string& matrixDescription,
                                                                                   const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                                   const Teuchos::RCP<Teuchos::FancyOStream>& dbg) {
  trcp_tcomm_t comm = MatrixMarketWriter::getComm(X.getMap());
  const int myRank  = MatrixMarketWriter::getRank(comm);
 
  auto out = MatrixMarketWriter::openOutFileOnRankZero(comm, filename, myRank, true);
 
  writeDense(out, X, matrixName, matrixDescription, err, dbg);
  // We can rely on the destructor of the output stream to close
  // the file on scope exit, even if writeDense() throws an
  // exception.
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeDenseFile(const std::string& filename,
                                                                                   const Teuchos::RCP<const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type>& X,
                                                                                   const std::string& matrixName,
                                                                                   const std::string& matrixDescription,
                                                                                   const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                                   const Teuchos::RCP<Teuchos::FancyOStream>& dbg) {
  TEUCHOS_TEST_FOR_EXCEPTION(
      X.is_null(), std::invalid_argument,
      "Tpetra::MatrixMarket::"
      "writeDenseFile: The input MultiVector X is null.");
  writeDenseFile(filename, *X, matrixName, matrixDescription, err, dbg);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeDenseFile(const std::string& filename,
                                                                                   const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type& X,
                                                                                   const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                                   const Teuchos::RCP<Teuchos::FancyOStream>& dbg) {
  writeDenseFile(filename, X, "", "", err, dbg);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeDenseFile(const std::string& filename,
                                                                                   const Teuchos::RCP<const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type>& X,
                                                                                   const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                                   const Teuchos::RCP<Teuchos::FancyOStream>& dbg) {
  TEUCHOS_TEST_FOR_EXCEPTION(
      X.is_null(), std::invalid_argument,
      "Tpetra::MatrixMarket::"
      "writeDenseFile: The input MultiVector X is null.");
  writeDenseFile(filename, *X, err, dbg);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeDense(std::ostream& out,
                                                                               const multivector_type& X,
                                                                               const std::string& matrixName,
                                                                               const std::string& matrixDescription,
                                                                               const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                               const Teuchos::RCP<Teuchos::FancyOStream>& dbg) {
  using std::endl;
  using Teuchos::Comm;
  using Teuchos::outArg;
  using Teuchos::REDUCE_MAX;
  using Teuchos::reduceAll;
 
  trcp_tcomm_t comm = MatrixMarketWriter::getComm(X.getMap());
  const int myRank  = MatrixMarketWriter::getRank(comm);
 
  // If the caller provides a nonnull debug output stream, we
  // print debugging output to it.  This is a local thing; we
  // don't have to check across processes.
  const bool debug = !dbg.is_null();
  if (debug) {
    dbg->pushTab();
    std::ostringstream os;
    os << myRank << ": writeDense" << endl;
    *dbg << os.str();
    dbg->pushTab();
  }
  // Print the Matrix Market header.
  writeDenseHeader(out, X, matrixName, matrixDescription, err, dbg);
 
  // Print each column one at a time.  This is a (perhaps)
  // temporary fix for Bug 6288.
  const size_t numVecs = X.getNumVectors();
  for (size_t j = 0; j < numVecs; ++j) {
    writeDenseColumn(out, *(X.getVector(j)), err, dbg);
  }
 
  if (debug) {
    dbg->popTab();
    std::ostringstream os;
    os << myRank << ": writeDense: Done" << endl;
    *dbg << os.str();
    dbg->popTab();
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
std::ofstream MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::openOutFileOnRankZero(
    const trcp_tcomm_t& comm,
    const std::string& filename, const int rank, const bool safe,
    const std::ios_base::openmode mode) {
  // Placeholder for the output stream.
  std::ofstream out;
 
  // State that will make all ranks throw if the root rank wasn't able to open the stream (using @c int for broadcasting).
  int all_should_stop = 0;
 
  // Try to open the file and update the state.
  if (rank == 0) {
    out.open(filename, mode);
    all_should_stop = !out && safe;
  }
 
  // Broadcast the stream state and throw from all ranks if needed.
  if (comm) Teuchos::broadcast(*comm, 0, &all_should_stop);
 
  TEUCHOS_TEST_FOR_EXCEPTION(
      all_should_stop,
      std::runtime_error,
      "Could not open output file '" + filename + "' on root rank 0.");
 
  return out;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeDenseHeader(std::ostream& out,
                                                                                     const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type& X,
                                                                                     const std::string& matrixName,
                                                                                     const std::string& matrixDescription,
                                                                                     const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                                     const Teuchos::RCP<Teuchos::FancyOStream>& dbg) {
  using std::endl;
  using Teuchos::Comm;
  using Teuchos::outArg;
  using Teuchos::RCP;
  using Teuchos::REDUCE_MAX;
  using Teuchos::reduceAll;
  typedef Teuchos::ScalarTraits<scalar_type> STS;
  const char prefix[] = "Tpetra::MatrixMarket::writeDenseHeader: ";
 
  trcp_tcomm_t comm = MatrixMarketWriter::getComm(X.getMap());
  const int myRank  = MatrixMarketWriter::getRank(comm);
  int lclErr        = 0;  // whether this MPI process has seen an error
  int gblErr        = 0;  // whether we know if some MPI process has seen an error
 
  // If the caller provides a nonnull debug output stream, we
  // print debugging output to it.  This is a local thing; we
  // don't have to check across processes.
  const bool debug = !dbg.is_null();
 
  if (debug) {
    dbg->pushTab();
    std::ostringstream os;
    os << myRank << ": writeDenseHeader" << endl;
    *dbg << os.str();
    dbg->pushTab();
  }
 
  //
  // Process 0: Write the MatrixMarket header.
  //
  if (myRank == 0) {
    try {
      // Print the Matrix Market header.  MultiVector stores data
      // nonsymmetrically, hence "general" in the banner line.
      // Print first to a temporary string output stream, and then
      // write it to the main output stream, so that at least the
      // header output has transactional semantics.  We can't
      // guarantee transactional semantics for the whole output,
      // since that would not be memory scalable.  (This could be
      // done in the file system by using a temporary file; we
      // don't do this, but users could.)
      std::ostringstream hdr;
      {
        std::string dataType;
        if (STS::isComplex) {
          dataType = "complex";
        } else if (STS::isOrdinal) {
          dataType = "integer";
        } else {
          dataType = "real";
        }
        hdr << "%%MatrixMarket matrix array " << dataType << " general"
            << endl;
      }
 
      // Print comments (the matrix name and / or description).
      if (matrixName != "") {
        printAsComment(hdr, matrixName);
      }
      if (matrixDescription != "") {
        printAsComment(hdr, matrixDescription);
      }
      // Print the Matrix Market dimensions header for dense matrices.
      hdr << X.getGlobalLength() << " " << X.getNumVectors() << endl;
 
      // Write the MatrixMarket header to the output stream.
      out << hdr.str();
    } catch (std::exception& e) {
      if (!err.is_null()) {
        *err << prefix << "While writing the Matrix Market header, "
                          "Process 0 threw an exception: "
             << e.what() << endl;
      }
      lclErr = 1;
    }
  }  // if I am Process 0
 
  // Establish global agreement on the error state.  It wouldn't
  // be good for other processes to keep going, if Process 0
  // finds out that it can't write to the given output stream.
  reduceAll<int, int>(*comm, REDUCE_MAX, lclErr, outArg(gblErr));
  TEUCHOS_TEST_FOR_EXCEPTION(
      gblErr == 1, std::runtime_error, prefix << "Some error occurred "
                                                 "which prevented this method from completing.");
 
  if (debug) {
    dbg->popTab();
    *dbg << myRank << ": writeDenseHeader: Done" << endl;
    dbg->popTab();
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeDenseColumn(std::ostream& out,
                                                                                     const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type& X,
                                                                                     const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                                     const Teuchos::RCP<Teuchos::FancyOStream>& dbg) {
  using std::endl;
  using Teuchos::arcp;
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::Comm;
  using Teuchos::CommRequest;
  using Teuchos::ireceive;
  using Teuchos::isend;
  using Teuchos::outArg;
  using Teuchos::RCP;
  using Teuchos::REDUCE_MAX;
  using Teuchos::reduceAll;
  using Teuchos::TypeNameTraits;
  using Teuchos::wait;
  typedef Teuchos::ScalarTraits<scalar_type> STS;
 
  const Comm<int>& comm = *(X.getMap()->getComm());
  const int myRank      = comm.getRank();
  const int numProcs    = comm.getSize();
  int lclErr            = 0;  // whether this MPI process has seen an error
  int gblErr            = 0;  // whether we know if some MPI process has seen an error
 
  // If the caller provides a nonnull debug output stream, we
  // print debugging output to it.  This is a local thing; we
  // don't have to check across processes.
  const bool debug = !dbg.is_null();
 
  if (debug) {
    dbg->pushTab();
    std::ostringstream os;
    os << myRank << ": writeDenseColumn" << endl;
    *dbg << os.str();
    dbg->pushTab();
  }
 
  // Make the output stream write floating-point numbers in
  // scientific notation.  It will politely put the output
  // stream back to its state on input, when this scope
  // terminates.
  Teuchos::SetScientific<scalar_type> sci(out);
 
  const size_t myNumRows = X.getLocalLength();
  const size_t numCols   = X.getNumVectors();
  // Use a different tag for the "size" messages than for the
  // "data" messages, in order to help us debug any mix-ups.
  const int sizeTag = 1337;
  const int dataTag = 1338;
 
  // Process 0 pipelines nonblocking receives with file output.
  //
  // Constraints:
  //   - Process 0 can't post a receive for another process'
  //     actual data, until it posts and waits on the receive
  //     from that process with the amount of data to receive.
  //     (We could just post receives with a max data size, but
  //     I feel uncomfortable about that.)
  //   - The C++ standard library doesn't allow nonblocking
  //     output to an std::ostream.  (Thus, we have to start a
  //     receive or send before starting the write, and hope
  //     that MPI completes it in the background.)
  //
  // Process 0: Post receive-size receives from Processes 1 and 2.
  // Process 1: Post send-size send to Process 0.
  // Process 2: Post send-size send to Process 0.
  //
  // All processes: Pack my entries.
  //
  // Process 1:
  //   - Post send-data send to Process 0.
  //   - Wait on my send-size send to Process 0.
  //
  // Process 0:
  //   - Print MatrixMarket header.
  //   - Print my entries.
  //   - Wait on receive-size receive from Process 1.
  //   - Post receive-data receive from Process 1.
  //
  // For each process p = 1, 2, ... numProcs-1:
  //   If I am Process 0:
  //     - Post receive-size receive from Process p + 2
  //     - Wait on receive-size receive from Process p + 1
  //     - Post receive-data receive from Process p + 1
  //     - Wait on receive-data receive from Process p
  //     - Write data from Process p.
  //   Else if I am Process p:
  //     - Wait on my send-data send.
  //   Else if I am Process p+1:
  //     - Post send-data send to Process 0.
  //     - Wait on my send-size send.
  //   Else if I am Process p+2:
  //     - Post send-size send to Process 0.
  //
  // Pipelining has three goals here:
  //   1. Overlap communication (the receives) with file I/O
  //   2. Give Process 0 a chance to prepost some receives,
  //      before sends show up, by packing local data before
  //      posting sends
  //   3. Don't post _all_ receives or _all_ sends, because that
  //      wouldn't be memory scalable.  (Just because we can't
  //      see how much memory MPI consumes, doesn't mean that it
  //      doesn't consume any!)
 
  // These are used on every process.  sendReqSize[0] holds the
  // number of rows on this process, and sendReqBuf holds this
  // process' data.  Process 0 packs into sendReqBuf, but
  // doesn't send; it only uses that for printing.  All other
  // processes send both of these to Process 0.
  RCP<CommRequest<int>> sendReqSize, sendReqData;
 
  // These are used only on Process 0, for received data.  Keep
  // 3 of each, and treat the arrays as circular buffers.  When
  // receiving from Process p, the corresponding array index
  // here is p % 3.
  Array<ArrayRCP<size_t>> recvSizeBufs(3);
  Array<ArrayRCP<scalar_type>> recvDataBufs(3);
  Array<RCP<CommRequest<int>>> recvSizeReqs(3);
  Array<RCP<CommRequest<int>>> recvDataReqs(3);
 
  // Buffer for nonblocking send of the "send size."
  ArrayRCP<size_t> sendDataSize(1);
  sendDataSize[0] = myNumRows;
 
  if (myRank == 0) {
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Post receive-size receives from "
                      "Procs 1 and 2: tag = "
         << sizeTag << endl;
      *dbg << os.str();
    }
    // Process 0: Post receive-size receives from Processes 1 and 2.
    recvSizeBufs[0].resize(1);
    // Set these three to an invalid value as a flag.  If we
    // don't get these messages, then the invalid value will
    // remain, so we can test for it.
    (recvSizeBufs[0])[0] = Teuchos::OrdinalTraits<size_t>::invalid();
    recvSizeBufs[1].resize(1);
    (recvSizeBufs[1])[0] = Teuchos::OrdinalTraits<size_t>::invalid();
    recvSizeBufs[2].resize(1);
    (recvSizeBufs[2])[0] = Teuchos::OrdinalTraits<size_t>::invalid();
    if (numProcs > 1) {
      recvSizeReqs[1] =
          ireceive<int, size_t>(recvSizeBufs[1], 1, sizeTag, comm);
    }
    if (numProcs > 2) {
      recvSizeReqs[2] =
          ireceive<int, size_t>(recvSizeBufs[2], 2, sizeTag, comm);
    }
  } else if (myRank == 1 || myRank == 2) {
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Post send-size send: size = "
         << sendDataSize[0] << ", tag = " << sizeTag << endl;
      *dbg << os.str();
    }
    // Prime the pipeline by having Processes 1 and 2 start
    // their send-size sends.  We don't want _all_ the processes
    // to start their send-size sends, because that wouldn't be
    // memory scalable.
    sendReqSize = isend<int, size_t>(sendDataSize, 0, sizeTag, comm);
  } else {
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Not posting my send-size send yet" << endl;
      *dbg << os.str();
    }
  }
 
  //
  // Pack my entries, in column-major order.
  //
  if (debug) {
    std::ostringstream os;
    os << myRank << ": Pack my entries" << endl;
    *dbg << os.str();
  }
  ArrayRCP<scalar_type> sendDataBuf;
  try {
    sendDataBuf = arcp<scalar_type>(myNumRows * numCols);
    X.get1dCopy(sendDataBuf(), myNumRows);
  } catch (std::exception& e) {
    lclErr = 1;
    if (!err.is_null()) {
      std::ostringstream os;
      os << "Process " << myRank << ": Attempt to pack my MultiVector "
                                    "entries threw an exception: "
         << e.what() << endl;
      *err << os.str();
    }
  }
  if (debug) {
    std::ostringstream os;
    os << myRank << ": Done packing my entries" << endl;
    *dbg << os.str();
  }
 
  //
  // Process 1: post send-data send to Process 0.
  //
  if (myRank == 1) {
    if (debug) {
      *dbg << myRank << ": Post send-data send: tag = " << dataTag
           << endl;
    }
    sendReqData = isend<int, scalar_type>(sendDataBuf, 0, dataTag, comm);
  }
 
  //
  // Process 0: Write my entries.
  //
  if (myRank == 0) {
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Write my entries" << endl;
      *dbg << os.str();
    }
 
    // Write Process 0's data to the output stream.
    // Matrix Market prints dense matrices in column-major order.
    const size_t printNumRows              = myNumRows;
    ArrayView<const scalar_type> printData = sendDataBuf();
    const size_t printStride               = printNumRows;
    if (static_cast<size_t>(printData.size()) < printStride * numCols) {
      lclErr = 1;
      if (!err.is_null()) {
        std::ostringstream os;
        os << "Process " << myRank << ": My MultiVector data's size "
           << printData.size() << " does not match my local dimensions "
           << printStride << " x " << numCols << "." << endl;
        *err << os.str();
      }
    } else {
      // Matrix Market dense format wants one number per line.
      // It wants each complex number as two real numbers (real
      // resp. imaginary parts) with a space between.
      for (size_t col = 0; col < numCols; ++col) {
        for (size_t row = 0; row < printNumRows; ++row) {
          if (STS::isComplex) {
            out << STS::real(printData[row + col * printStride]) << " "
                << STS::imag(printData[row + col * printStride]) << endl;
          } else {
            out << printData[row + col * printStride] << endl;
          }
        }
      }
    }
  }
 
  if (myRank == 0) {
    // Wait on receive-size receive from Process 1.
    const int recvRank   = 1;
    const int circBufInd = recvRank % 3;
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Wait on receive-size receive from Process "
         << recvRank << endl;
      *dbg << os.str();
    }
    if (numProcs > 1) {
      wait<int>(comm, outArg(recvSizeReqs[circBufInd]));
 
      // We received the number of rows of data.  (The data
      // come in two columns.)
      size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
      if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid()) {
        lclErr = 1;
        if (!err.is_null()) {
          std::ostringstream os;
          os << myRank << ": Result of receive-size receive from Process "
             << recvRank << " is Teuchos::OrdinalTraits<size_t>::invalid() "
             << "= " << Teuchos::OrdinalTraits<size_t>::invalid() << ".  "
                                                                     "This should never happen, and suggests that the receive never "
                                                                     "got posted.  Please report this bug to the Tpetra developers."
             << endl;
          *err << os.str();
        }
 
        // If we're going to continue after error, set the
        // number of rows to receive to a reasonable size.  This
        // may cause MPI_ERR_TRUNCATE if the sending process is
        // sending more than 0 rows, but that's better than MPI
        // overflowing due to the huge positive value that is
        // Teuchos::OrdinalTraits<size_t>::invalid().
        recvNumRows = 0;
      }
 
      // Post receive-data receive from Process 1.
      recvDataBufs[circBufInd].resize(recvNumRows * numCols);
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Post receive-data receive from Process "
           << recvRank << ": tag = " << dataTag << ", buffer size = "
           << recvDataBufs[circBufInd].size() << endl;
        *dbg << os.str();
      }
      if (!recvSizeReqs[circBufInd].is_null()) {
        lclErr = 1;
        if (!err.is_null()) {
          std::ostringstream os;
          os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
                                                             "null, before posting the receive-data receive from Process "
             << recvRank << ".  This should never happen.  Please report "
                            "this bug to the Tpetra developers."
             << endl;
          *err << os.str();
        }
      }
      recvDataReqs[circBufInd] =
          ireceive<int, scalar_type>(recvDataBufs[circBufInd],
                                     recvRank, dataTag, comm);
    }  // numProcs > 1
  } else if (myRank == 1) {
    // Wait on my send-size send.
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Wait on my send-size send" << endl;
      *dbg << os.str();
    }
    wait<int>(comm, outArg(sendReqSize));
  }
 
  //
  // Pipeline loop
  //
  for (int p = 1; p < numProcs; ++p) {
    if (myRank == 0) {
      if (p + 2 < numProcs) {
        // Post receive-size receive from Process p + 2.
        const int recvRank   = p + 2;
        const int circBufInd = recvRank % 3;
        if (debug) {
          std::ostringstream os;
          os << myRank << ": Post receive-size receive from Process "
             << recvRank << ": tag = " << sizeTag << endl;
          *dbg << os.str();
        }
        if (!recvSizeReqs[circBufInd].is_null()) {
          lclErr = 1;
          if (!err.is_null()) {
            std::ostringstream os;
            os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
               << "null, for the receive-size receive from Process "
               << recvRank << "!  This may mean that this process never "
               << "finished waiting for the receive from Process "
               << (recvRank - 3) << "." << endl;
            *err << os.str();
          }
        }
        recvSizeReqs[circBufInd] =
            ireceive<int, size_t>(recvSizeBufs[circBufInd],
                                  recvRank, sizeTag, comm);
      }
 
      if (p + 1 < numProcs) {
        const int recvRank   = p + 1;
        const int circBufInd = recvRank % 3;
 
        // Wait on receive-size receive from Process p + 1.
        if (debug) {
          std::ostringstream os;
          os << myRank << ": Wait on receive-size receive from Process "
             << recvRank << endl;
          *dbg << os.str();
        }
        wait<int>(comm, outArg(recvSizeReqs[circBufInd]));
 
        // We received the number of rows of data.  (The data
        // come in two columns.)
        size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
        if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid()) {
          lclErr = 1;
          if (!err.is_null()) {
            std::ostringstream os;
            os << myRank << ": Result of receive-size receive from Process "
               << recvRank << " is Teuchos::OrdinalTraits<size_t>::invalid() "
               << "= " << Teuchos::OrdinalTraits<size_t>::invalid() << ".  "
                                                                       "This should never happen, and suggests that the receive never "
                                                                       "got posted.  Please report this bug to the Tpetra developers."
               << endl;
            *err << os.str();
          }
          // If we're going to continue after error, set the
          // number of rows to receive to a reasonable size.
          // This may cause MPI_ERR_TRUNCATE if the sending
          // process sends more than 0 rows, but that's better
          // than MPI overflowing due to the huge positive value
          // Teuchos::OrdinalTraits<size_t>::invalid().
          recvNumRows = 0;
        }
 
        // Post receive-data receive from Process p + 1.
        recvDataBufs[circBufInd].resize(recvNumRows * numCols);
        if (debug) {
          std::ostringstream os;
          os << myRank << ": Post receive-data receive from Process "
             << recvRank << ": tag = " << dataTag << ", buffer size = "
             << recvDataBufs[circBufInd].size() << endl;
          *dbg << os.str();
        }
        if (!recvDataReqs[circBufInd].is_null()) {
          lclErr = 1;
          if (!err.is_null()) {
            std::ostringstream os;
            os << myRank << ": recvDataReqs[" << circBufInd << "] is not "
               << "null, for the receive-data receive from Process "
               << recvRank << "!  This may mean that this process never "
               << "finished waiting for the receive from Process "
               << (recvRank - 3) << "." << endl;
            *err << os.str();
          }
        }
        recvDataReqs[circBufInd] =
            ireceive<int, scalar_type>(recvDataBufs[circBufInd],
                                       recvRank, dataTag, comm);
      }
 
      // Wait on receive-data receive from Process p.
      const int recvRank   = p;
      const int circBufInd = recvRank % 3;
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Wait on receive-data receive from Process "
           << recvRank << endl;
        *dbg << os.str();
      }
      wait<int>(comm, outArg(recvDataReqs[circBufInd]));
 
      // Write Process p's data.  Number of rows lives in
      // recvSizeBufs[circBufInd], and the actual data live in
      // recvDataBufs[circBufInd].  Do this after posting receives,
      // in order to expose overlap of comm. with file I/O.
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Write entries from Process " << recvRank
           << endl;
        *dbg << os.str() << endl;
      }
      size_t printNumRows = (recvSizeBufs[circBufInd])[0];
      if (printNumRows == Teuchos::OrdinalTraits<size_t>::invalid()) {
        lclErr = 1;
        if (!err.is_null()) {
          std::ostringstream os;
          os << myRank << ": Result of receive-size receive from Process "
             << recvRank << " was Teuchos::OrdinalTraits<size_t>::"
                            "invalid() = "
             << Teuchos::OrdinalTraits<size_t>::invalid()
             << ".  This should never happen, and suggests that its "
                "receive-size receive was never posted.  "
                "Please report this bug to the Tpetra developers."
             << endl;
          *err << os.str();
        }
        // If we're going to continue after error, set the
        // number of rows to print to a reasonable size.
        printNumRows = 0;
      }
      if (printNumRows > 0 && recvDataBufs[circBufInd].is_null()) {
        lclErr = 1;
        if (!err.is_null()) {
          std::ostringstream os;
          os << myRank << ": Result of receive-size receive from Proc "
             << recvRank << " was " << printNumRows << " > 0, but "
                                                       "recvDataBufs["
             << circBufInd << "] is null.  This should "
                              "never happen.  Please report this bug to the Tpetra "
                              "developers."
             << endl;
          *err << os.str();
        }
        // If we're going to continue after error, set the
        // number of rows to print to a reasonable size.
        printNumRows = 0;
      }
 
      // Write the received data to the output stream.
      // Matrix Market prints dense matrices in column-major order.
      ArrayView<const scalar_type> printData = (recvDataBufs[circBufInd])();
      const size_t printStride               = printNumRows;
      // Matrix Market dense format wants one number per line.
      // It wants each complex number as two real numbers (real
      // resp. imaginary parts) with a space between.
      for (size_t col = 0; col < numCols; ++col) {
        for (size_t row = 0; row < printNumRows; ++row) {
          if (STS::isComplex) {
            out << STS::real(printData[row + col * printStride]) << " "
                << STS::imag(printData[row + col * printStride]) << endl;
          } else {
            out << printData[row + col * printStride] << endl;
          }
        }
      }
    } else if (myRank == p) {  // Process p
      // Wait on my send-data send.
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Wait on my send-data send" << endl;
        *dbg << os.str();
      }
      wait<int>(comm, outArg(sendReqData));
    } else if (myRank == p + 1) {  // Process p + 1
      // Post send-data send to Process 0.
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Post send-data send: tag = " << dataTag
           << endl;
        *dbg << os.str();
      }
      sendReqData = isend<int, scalar_type>(sendDataBuf, 0, dataTag, comm);
      // Wait on my send-size send.
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Wait on my send-size send" << endl;
        *dbg << os.str();
      }
      wait<int>(comm, outArg(sendReqSize));
    } else if (myRank == p + 2) {  // Process p + 2
      // Post send-size send to Process 0.
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Post send-size send: size = "
           << sendDataSize[0] << ", tag = " << sizeTag << endl;
        *dbg << os.str();
      }
      sendReqSize = isend<int, size_t>(sendDataSize, 0, sizeTag, comm);
    }
  }
 
  // Establish global agreement on the error state.
  reduceAll<int, int>(comm, REDUCE_MAX, lclErr, outArg(gblErr));
  TEUCHOS_TEST_FOR_EXCEPTION(
      gblErr == 1, std::runtime_error,
      "Tpetra::MatrixMarket::writeDense "
      "experienced some kind of error and was unable to complete.");
 
  if (debug) {
    dbg->popTab();
    *dbg << myRank << ": writeDenseColumn: Done" << endl;
    dbg->popTab();
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeDense(std::ostream& out,
                                                                               const Teuchos::RCP<const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type>& X,
                                                                               const std::string& matrixName,
                                                                               const std::string& matrixDescription,
                                                                               const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                               const Teuchos::RCP<Teuchos::FancyOStream>& dbg) {
  TEUCHOS_TEST_FOR_EXCEPTION(
      X.is_null(), std::invalid_argument,
      "Tpetra::MatrixMarket::"
      "writeDense: The input MultiVector X is null.");
  writeDense(out, *X, matrixName, matrixDescription, err, dbg);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeDense(std::ostream& out,
                                                                               const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type& X,
                                                                               const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                               const Teuchos::RCP<Teuchos::FancyOStream>& dbg) {
  writeDense(out, X, "", "", err, dbg);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeDense(std::ostream& out,
                                                                               const Teuchos::RCP<const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::multivector_type>& X,
                                                                               const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                               const Teuchos::RCP<Teuchos::FancyOStream>& dbg) {
  TEUCHOS_TEST_FOR_EXCEPTION(
      X.is_null(), std::invalid_argument,
      "Tpetra::MatrixMarket::"
      "writeDense: The input MultiVector X is null.");
  writeDense(out, *X, "", "", err, dbg);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeMap(std::ostream& out, const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type& map, const bool debug) {
  Teuchos::RCP<Teuchos::FancyOStream> err =
      Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
  writeMap(out, map, err, debug);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeMap(std::ostream& out,
                                                                             const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type& map,
                                                                             const Teuchos::RCP<Teuchos::FancyOStream>& err,
                                                                             const bool debug) {
  using std::endl;
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::Comm;
  using Teuchos::CommRequest;
  using Teuchos::ireceive;
  using Teuchos::isend;
  using Teuchos::RCP;
  using Teuchos::TypeNameTraits;
  using Teuchos::wait;
  typedef global_ordinal_type GO;
  typedef int pid_type;
 
  // Treat the Map as a 1-column "multivector."  This differs
  // from the previous two-column format, in which column 0 held
  // the GIDs, and column 1 held the corresponding PIDs.  It
  // differs because printing that format requires knowing the
  // entire first column -- that is, all the GIDs -- in advance.
  // Sending messages from each process one at a time saves
  // memory, but it means that Process 0 doesn't ever have all
  // the GIDs at once.
  //
  // We pack the entries as ptrdiff_t, since this should be the
  // biggest signed built-in integer type that can hold any GO
  // or pid_type (= int) quantity without overflow.  Test this
  // assumption at run time.
  typedef ptrdiff_t int_type;
  TEUCHOS_TEST_FOR_EXCEPTION(
      sizeof(GO) > sizeof(int_type), std::logic_error,
      "The global ordinal type GO=" << TypeNameTraits<GO>::name()
                                    << " is too big for ptrdiff_t.  sizeof(GO) = " << sizeof(GO)
                                    << " > sizeof(ptrdiff_t) = " << sizeof(ptrdiff_t) << ".");
  TEUCHOS_TEST_FOR_EXCEPTION(
      sizeof(pid_type) > sizeof(int_type), std::logic_error,
      "The (MPI) process rank type pid_type=" << TypeNameTraits<pid_type>::name() << " is too big for ptrdiff_t.  "
                                                                                     "sizeof(pid_type) = "
                                              << sizeof(pid_type) << " > sizeof(ptrdiff_t)"
                                                                     " = "
                                              << sizeof(ptrdiff_t) << ".");
 
  const Comm<int>& comm = *(map.getComm());
  const int myRank      = comm.getRank();
  const int numProcs    = comm.getSize();
 
  if (!err.is_null()) {
    err->pushTab();
  }
  if (debug) {
    std::ostringstream os;
    os << myRank << ": writeMap" << endl;
    *err << os.str();
  }
  if (!err.is_null()) {
    err->pushTab();
  }
 
  const size_t myNumRows = map.getLocalNumElements();
  // Use a different tag for the "size" messages than for the
  // "data" messages, in order to help us debug any mix-ups.
  const int sizeTag = 1337;
  const int dataTag = 1338;
 
  // Process 0 pipelines nonblocking receives with file output.
  //
  // Constraints:
  //   - Process 0 can't post a receive for another process'
  //     actual data, until it posts and waits on the receive
  //     from that process with the amount of data to receive.
  //     (We could just post receives with a max data size, but
  //     I feel uncomfortable about that.)
  //   - The C++ standard library doesn't allow nonblocking
  //     output to an std::ostream.
  //
  // Process 0: Post receive-size receives from Processes 1 and 2.
  // Process 1: Post send-size send to Process 0.
  // Process 2: Post send-size send to Process 0.
  //
  // All processes: Pack my GIDs and PIDs.
  //
  // Process 1:
  //   - Post send-data send to Process 0.
  //   - Wait on my send-size send to Process 0.
  //
  // Process 0:
  //   - Print MatrixMarket header.
  //   - Print my GIDs and PIDs.
  //   - Wait on receive-size receive from Process 1.
  //   - Post receive-data receive from Process 1.
  //
  // For each process p = 1, 2, ... numProcs-1:
  //   If I am Process 0:
  //     - Post receive-size receive from Process p + 2
  //     - Wait on receive-size receive from Process p + 1
  //     - Post receive-data receive from Process p + 1
  //     - Wait on receive-data receive from Process p
  //     - Write data from Process p.
  //   Else if I am Process p:
  //     - Wait on my send-data send.
  //   Else if I am Process p+1:
  //     - Post send-data send to Process 0.
  //     - Wait on my send-size send.
  //   Else if I am Process p+2:
  //     - Post send-size send to Process 0.
  //
  // Pipelining has three goals here:
  //   1. Overlap communication (the receives) with file I/O
  //   2. Give Process 0 a chance to prepost some receives,
  //      before sends show up, by packing local data before
  //      posting sends
  //   3. Don't post _all_ receives or _all_ sends, because that
  //      wouldn't be memory scalable.  (Just because we can't
  //      see how much memory MPI consumes, doesn't mean that it
  //      doesn't consume any!)
 
  // These are used on every process.  sendReqSize[0] holds the
  // number of rows on this process, and sendReqBuf holds this
  // process' data.  Process 0 packs into sendReqBuf, but
  // doesn't send; it only uses that for printing.  All other
  // processes send both of these to Process 0.
  RCP<CommRequest<int>> sendReqSize, sendReqData;
 
  // These are used only on Process 0, for received data.  Keep
  // 3 of each, and treat the arrays as circular buffers.  When
  // receiving from Process p, the corresponding array index
  // here is p % 3.
  Array<ArrayRCP<int_type>> recvSizeBufs(3);
  Array<ArrayRCP<int_type>> recvDataBufs(3);
  Array<RCP<CommRequest<int>>> recvSizeReqs(3);
  Array<RCP<CommRequest<int>>> recvDataReqs(3);
 
  // Buffer for nonblocking send of the "send size."
  ArrayRCP<int_type> sendDataSize(1);
  sendDataSize[0] = myNumRows;
 
  if (myRank == 0) {
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Post receive-size receives from "
                      "Procs 1 and 2: tag = "
         << sizeTag << endl;
      *err << os.str();
    }
    // Process 0: Post receive-size receives from Processes 1 and 2.
    recvSizeBufs[0].resize(1);
    (recvSizeBufs[0])[0] = -1;  // error flag
    recvSizeBufs[1].resize(1);
    (recvSizeBufs[1])[0] = -1;  // error flag
    recvSizeBufs[2].resize(1);
    (recvSizeBufs[2])[0] = -1;  // error flag
    if (numProcs > 1) {
      recvSizeReqs[1] =
          ireceive<int, int_type>(recvSizeBufs[1], 1, sizeTag, comm);
    }
    if (numProcs > 2) {
      recvSizeReqs[2] =
          ireceive<int, int_type>(recvSizeBufs[2], 2, sizeTag, comm);
    }
  } else if (myRank == 1 || myRank == 2) {
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Post send-size send: size = "
         << sendDataSize[0] << ", tag = " << sizeTag << endl;
      *err << os.str();
    }
    // Prime the pipeline by having Processes 1 and 2 start
    // their send-size sends.  We don't want _all_ the processes
    // to start their send-size sends, because that wouldn't be
    // memory scalable.
    sendReqSize = isend<int, int_type>(sendDataSize, 0, sizeTag, comm);
  } else {
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Not posting my send-size send yet" << endl;
      *err << os.str();
    }
  }
 
  //
  // Pack my GIDs and PIDs.  Each (GID,PID) pair gets packed
  // consecutively, for better locality.
  //
 
  if (debug) {
    std::ostringstream os;
    os << myRank << ": Pack my GIDs and PIDs" << endl;
    *err << os.str();
  }
 
  ArrayRCP<int_type> sendDataBuf(myNumRows * 2);
 
  if (map.isContiguous()) {
    const int_type myMinGblIdx =
        static_cast<int_type>(map.getMinGlobalIndex());
    for (size_t k = 0; k < myNumRows; ++k) {
      const int_type gid     = myMinGblIdx + static_cast<int_type>(k);
      const int_type pid     = static_cast<int_type>(myRank);
      sendDataBuf[2 * k]     = gid;
      sendDataBuf[2 * k + 1] = pid;
    }
  } else {
    ArrayView<const GO> myGblInds = map.getLocalElementList();
    for (size_t k = 0; k < myNumRows; ++k) {
      const int_type gid     = static_cast<int_type>(myGblInds[k]);
      const int_type pid     = static_cast<int_type>(myRank);
      sendDataBuf[2 * k]     = gid;
      sendDataBuf[2 * k + 1] = pid;
    }
  }
 
  if (debug) {
    std::ostringstream os;
    os << myRank << ": Done packing my GIDs and PIDs" << endl;
    *err << os.str();
  }
 
  if (myRank == 1) {
    // Process 1: post send-data send to Process 0.
    if (debug) {
      *err << myRank << ": Post send-data send: tag = " << dataTag
           << endl;
    }
    sendReqData = isend<int, int_type>(sendDataBuf, 0, dataTag, comm);
  }
 
  if (myRank == 0) {
    if (debug) {
      *err << myRank << ": Write MatrixMarket header" << endl;
    }
 
    // Process 0: Write the MatrixMarket header.
    // Description section explains each column.
    std::ostringstream hdr;
 
    // Print the Matrix Market header.  MultiVector stores data
    // nonsymmetrically, hence "general" in the banner line.
    hdr << "%%MatrixMarket matrix array integer general" << endl
        << "% Format: Version 2.0" << endl
        << "%" << endl
        << "% This file encodes a Tpetra::Map." << endl
        << "% It is stored as a dense vector, with twice as many " << endl
        << "% entries as the global number of GIDs (global indices)." << endl
        << "% (GID, PID) pairs are stored contiguously, where the PID " << endl
        << "% is the rank of the process owning that GID." << endl
        << (2 * map.getGlobalNumElements()) << " " << 1 << endl;
    out << hdr.str();
 
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Write my GIDs and PIDs" << endl;
      *err << os.str();
    }
 
    // Write Process 0's data to the output stream.
    // Matrix Market prints dense matrices in column-major order.
    const int_type printNumRows         = myNumRows;
    ArrayView<const int_type> printData = sendDataBuf();
    for (int_type k = 0; k < printNumRows; ++k) {
      const int_type gid = printData[2 * k];
      const int_type pid = printData[2 * k + 1];
      out << gid << endl
          << pid << endl;
    }
  }
 
  if (myRank == 0) {
    // Wait on receive-size receive from Process 1.
    const int recvRank   = 1;
    const int circBufInd = recvRank % 3;
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Wait on receive-size receive from Process "
         << recvRank << endl;
      *err << os.str();
    }
    if (numProcs > 1) {
      wait<int>(comm, outArg(recvSizeReqs[circBufInd]));
 
      // We received the number of rows of data.  (The data
      // come in two columns.)
      const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
      if (debug && recvNumRows == -1) {
        std::ostringstream os;
        os << myRank << ": Result of receive-size receive from Process "
           << recvRank << " is -1.  This should never happen, and "
                          "suggests that the receive never got posted.  Please report "
                          "this bug to the Tpetra developers."
           << endl;
        *err << os.str();
      }
 
      // Post receive-data receive from Process 1.
      recvDataBufs[circBufInd].resize(recvNumRows * 2);
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Post receive-data receive from Process "
           << recvRank << ": tag = " << dataTag << ", buffer size = "
           << recvDataBufs[circBufInd].size() << endl;
        *err << os.str();
      }
      if (!recvSizeReqs[circBufInd].is_null()) {
        std::ostringstream os;
        os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
                                                           "null, before posting the receive-data receive from Process "
           << recvRank << ".  This should never happen.  Please report "
                          "this bug to the Tpetra developers."
           << endl;
        *err << os.str();
      }
      recvDataReqs[circBufInd] =
          ireceive<int, int_type>(recvDataBufs[circBufInd],
                                  recvRank, dataTag, comm);
    }  // numProcs > 1
  } else if (myRank == 1) {
    // Wait on my send-size send.
    if (debug) {
      std::ostringstream os;
      os << myRank << ": Wait on my send-size send" << endl;
      *err << os.str();
    }
    wait<int>(comm, outArg(sendReqSize));
  }
 
  //
  // Pipeline loop
  //
  for (int p = 1; p < numProcs; ++p) {
    if (myRank == 0) {
      if (p + 2 < numProcs) {
        // Post receive-size receive from Process p + 2.
        const int recvRank   = p + 2;
        const int circBufInd = recvRank % 3;
        if (debug) {
          std::ostringstream os;
          os << myRank << ": Post receive-size receive from Process "
             << recvRank << ": tag = " << sizeTag << endl;
          *err << os.str();
        }
        if (!recvSizeReqs[circBufInd].is_null()) {
          std::ostringstream os;
          os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
             << "null, for the receive-size receive from Process "
             << recvRank << "!  This may mean that this process never "
             << "finished waiting for the receive from Process "
             << (recvRank - 3) << "." << endl;
          *err << os.str();
        }
        recvSizeReqs[circBufInd] =
            ireceive<int, int_type>(recvSizeBufs[circBufInd],
                                    recvRank, sizeTag, comm);
      }
 
      if (p + 1 < numProcs) {
        const int recvRank   = p + 1;
        const int circBufInd = recvRank % 3;
 
        // Wait on receive-size receive from Process p + 1.
        if (debug) {
          std::ostringstream os;
          os << myRank << ": Wait on receive-size receive from Process "
             << recvRank << endl;
          *err << os.str();
        }
        wait<int>(comm, outArg(recvSizeReqs[circBufInd]));
 
        // We received the number of rows of data.  (The data
        // come in two columns.)
        const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
        if (debug && recvNumRows == -1) {
          std::ostringstream os;
          os << myRank << ": Result of receive-size receive from Process "
             << recvRank << " is -1.  This should never happen, and "
                            "suggests that the receive never got posted.  Please report "
                            "this bug to the Tpetra developers."
             << endl;
          *err << os.str();
        }
 
        // Post receive-data receive from Process p + 1.
        recvDataBufs[circBufInd].resize(recvNumRows * 2);
        if (debug) {
          std::ostringstream os;
          os << myRank << ": Post receive-data receive from Process "
             << recvRank << ": tag = " << dataTag << ", buffer size = "
             << recvDataBufs[circBufInd].size() << endl;
          *err << os.str();
        }
        if (!recvDataReqs[circBufInd].is_null()) {
          std::ostringstream os;
          os << myRank << ": recvDataReqs[" << circBufInd << "] is not "
             << "null, for the receive-data receive from Process "
             << recvRank << "!  This may mean that this process never "
             << "finished waiting for the receive from Process "
             << (recvRank - 3) << "." << endl;
          *err << os.str();
        }
        recvDataReqs[circBufInd] =
            ireceive<int, int_type>(recvDataBufs[circBufInd],
                                    recvRank, dataTag, comm);
      }
 
      // Wait on receive-data receive from Process p.
      const int recvRank   = p;
      const int circBufInd = recvRank % 3;
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Wait on receive-data receive from Process "
           << recvRank << endl;
        *err << os.str();
      }
      wait<int>(comm, outArg(recvDataReqs[circBufInd]));
 
      // Write Process p's data.  Number of rows lives in
      // recvSizeBufs[circBufInd], and the actual data live in
      // recvDataBufs[circBufInd].  Do this after posting receives,
      // in order to expose overlap of comm. with file I/O.
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Write GIDs and PIDs from Process "
           << recvRank << endl;
        *err << os.str() << endl;
      }
      const int_type printNumRows = (recvSizeBufs[circBufInd])[0];
      if (debug && printNumRows == -1) {
        std::ostringstream os;
        os << myRank << ": Result of receive-size receive from Process "
           << recvRank << " was -1.  This should never happen, and "
                          "suggests that its receive-size receive was never posted.  "
                          "Please report this bug to the Tpetra developers."
           << endl;
        *err << os.str();
      }
      if (debug && printNumRows > 0 && recvDataBufs[circBufInd].is_null()) {
        std::ostringstream os;
        os << myRank << ": Result of receive-size receive from Proc "
           << recvRank << " was " << printNumRows << " > 0, but "
                                                     "recvDataBufs["
           << circBufInd << "] is null.  This should "
                            "never happen.  Please report this bug to the Tpetra "
                            "developers."
           << endl;
        *err << os.str();
      }
      ArrayView<const int_type> printData = (recvDataBufs[circBufInd])();
      for (int_type k = 0; k < printNumRows; ++k) {
        const int_type gid = printData[2 * k];
        const int_type pid = printData[2 * k + 1];
        out << gid << endl
            << pid << endl;
      }
    } else if (myRank == p) {  // Process p
      // Wait on my send-data send.
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Wait on my send-data send" << endl;
        *err << os.str();
      }
      wait<int>(comm, outArg(sendReqData));
    } else if (myRank == p + 1) {  // Process p + 1
      // Post send-data send to Process 0.
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Post send-data send: tag = " << dataTag
           << endl;
        *err << os.str();
      }
      sendReqData = isend<int, int_type>(sendDataBuf, 0, dataTag, comm);
      // Wait on my send-size send.
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Wait on my send-size send" << endl;
        *err << os.str();
      }
      wait<int>(comm, outArg(sendReqSize));
    } else if (myRank == p + 2) {  // Process p + 2
      // Post send-size send to Process 0.
      if (debug) {
        std::ostringstream os;
        os << myRank << ": Post send-size send: size = "
           << sendDataSize[0] << ", tag = " << sizeTag << endl;
        *err << os.str();
      }
      sendReqSize = isend<int, int_type>(sendDataSize, 0, sizeTag, comm);
    }
  }
 
  if (!err.is_null()) {
    err->popTab();
  }
  if (debug) {
    *err << myRank << ": writeMap: Done" << endl;
  }
  if (!err.is_null()) {
    err->popTab();
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeMapFile(const std::string& filename,
                                                                                 const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::map_type& map) {
  const int myRank = map.getComm()->getRank();
 
  auto out = MatrixMarketWriter::openOutFileOnRankZero(map.getComm(), filename, myRank, true);
 
  writeMap(out, map);
  // We can rely on the destructor of the output stream to close
  // the file on scope exit, even if writeDense() throws an
  // exception.
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::printAsComment(std::ostream& out, const std::string& str) {
  using std::endl;
  std::istringstream inpstream(str);
  std::string line;
 
  while (getline(inpstream, line)) {
    if (!line.empty()) {
      // Note that getline() doesn't store '\n', so we have to
      // append the endline ourselves.
      if (line[0] == '%') {  // Line starts with a comment character.
        out << line << endl;
      } else {  // Line doesn't start with a comment character.
        out << "%% " << line << endl;
      }
    }
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeOperator(const std::string& fileName, typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::operator_type const& A) {
  Teuchos::ParameterList pl;
  writeOperator(fileName, A, pl);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeOperator(std::ostream& out, const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::operator_type& A) {
  Teuchos::ParameterList pl;
  writeOperator(out, A, pl);
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeOperator(const std::string& fileName,
                                                                                  const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::operator_type& A,
                                                                                  const Teuchos::ParameterList& params) {
  std::ofstream out;
  std::string tmpFile   = "__TMP__" + fileName;
  const int myRank      = A.getDomainMap()->getComm()->getRank();
  bool precisionChanged = false;
  int oldPrecision      = 0;
  // The number of nonzero entries in a Tpetra::Operator is
  // unknown until probing is completed.  In order to write a
  // MatrixMarket header, we write the matrix to a temporary
  // file.
  //
  // FIXME (mfh 23 May 2015) IT WASN'T MY IDEA TO WRITE TO A
  // TEMPORARY FILE.
  if (myRank == 0) {
    if (std::ifstream(tmpFile))
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
                                 "writeOperator: temporary file " << tmpFile << " already exists");
    out.open(tmpFile.c_str());
    if (params.isParameter("precision")) {
      oldPrecision     = out.precision(params.get<int>("precision"));
      precisionChanged = true;
    }
  }
 
  const std::string header = writeOperatorImpl(out, A, params);
 
  if (myRank == 0) {
    if (precisionChanged)
      out.precision(oldPrecision);
    out.close();
    out.open(fileName.c_str(), std::ios::binary);
    bool printMatrixMarketHeader = true;
    if (params.isParameter("print MatrixMarket header"))
      printMatrixMarketHeader = params.get<bool>("print MatrixMarket header");
    if (printMatrixMarketHeader && myRank == 0) {
      // Write header to final file.
      out << header;
    }
    // Append matrix from temporary to final file.
    std::ifstream src(tmpFile, std::ios_base::binary);
    out << src.rdbuf();
    src.close();
    // Delete the temporary file.
    remove(tmpFile.c_str());
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeOperator(std::ostream& out,
                                                                                  const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::operator_type& A,
                                                                                  const Teuchos::ParameterList& params) {
  const int myRank = A.getDomainMap()->getComm()->getRank();
 
  // The number of nonzero entries in a Tpetra::Operator is
  // unknown until probing is completed.  In order to write a
  // MatrixMarket header, we write the matrix to a temporary
  // output stream.
  //
  // NOTE (mfh 23 May 2015): Writing to a temporary output
  // stream may double the memory usage, depending on whether
  // 'out' is a file stream or an in-memory output stream (e.g.,
  // std::ostringstream).  It might be wise to use a temporary
  // file instead.  However, please look carefully at POSIX
  // functions for safe creation of temporary files.  Don't just
  // prepend "__TMP__" to the filename and hope for the best.
  // Furthermore, it should be valid to call the std::ostream
  // overload of this method even when Process 0 does not have
  // access to a file system.
  std::ostringstream tmpOut;
  if (myRank == 0) {
    if (params.isParameter("precision") && params.isType<int>("precision")) {
      (void)tmpOut.precision(params.get<int>("precision"));
    }
  }
 
  const std::string header = writeOperatorImpl(tmpOut, A, params);
 
  if (myRank == 0) {
    bool printMatrixMarketHeader = true;
    if (params.isParameter("print MatrixMarket header") &&
        params.isType<bool>("print MatrixMarket header")) {
      printMatrixMarketHeader = params.get<bool>("print MatrixMarket header");
    }
    if (printMatrixMarketHeader && myRank == 0) {
      out << header;  // write header to final output stream
    }
    // Append matrix from temporary output stream to final output stream.
    //
    // NOTE (mfh 23 May 2015) This might use a lot of memory.
    // However, we should not use temporary files in this
    // method.  Since it does not access the file system (unlike
    // the overload that takes a file name), it should not
    // require the file system at all.
    //
    // If memory usage becomes a problem, one thing we could do
    // is write the entries of the Operator one column (or a few
    // columns) at a time.  The Matrix Market sparse format does
    // not impose an order on its entries, so it would be OK to
    // write them in that order.
    out << tmpOut.str();
  }
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
std::string
MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeOperatorImpl(std::ostream& os,
                                                                                 const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::operator_type& A,
                                                                                 const Teuchos::ParameterList& params) {
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::RCP;
  using Teuchos::rcp;
 
  typedef local_ordinal_type LO;
  typedef global_ordinal_type GO;
  typedef Teuchos::OrdinalTraits<LO> TLOT;
  typedef Teuchos::OrdinalTraits<GO> TGOT;
  typedef Tpetra::Import<LO, GO, node_type> import_type;
  typedef Tpetra::MultiVector<GO, LO, GO, node_type> mv_type_go;
 
  const map_type& domainMap    = *(A.getDomainMap());
  RCP<const map_type> rangeMap = A.getRangeMap();
  trcp_tcomm_t comm            = rangeMap->getComm();
  const int myRank             = comm->getRank();
  const size_t numProcs        = comm->getSize();
 
  size_t numMVs = 10;
  if (params.isParameter("probing size"))
    numMVs = params.get<int>("probing size");
 
  GO globalNnz = 0;
  GO minColGid = domainMap.getMinAllGlobalIndex();
  GO maxColGid = domainMap.getMaxAllGlobalIndex();
  // Rather than replicating the domainMap on all processors, we instead
  // iterate from the min GID to the max GID.  If the map is gappy,
  // there will be invalid GIDs, i.e., GIDs no one has.  This will require
  // unnecessary matvecs  against potentially zero vectors.
  GO numGlobElts = maxColGid - minColGid + TGOT::one();
  GO numChunks   = numGlobElts / numMVs;
  GO rem         = numGlobElts % numMVs;
  GO indexBase   = rangeMap->getIndexBase();
 
  int offsetToUseInPrinting = 1 - indexBase;  // default is 1-based indexing
  if (params.isParameter("zero-based indexing")) {
    if (params.get<bool>("zero-based indexing") == true)
      offsetToUseInPrinting = -indexBase;  // If 0-based, use as-is. If 1-based, subtract 1.
  }
 
  // Create map that replicates the range map on pid 0 and is empty for all other pids
  size_t numLocalRangeEntries = rangeMap->getLocalNumElements();
 
  // Create contiguous source map
  RCP<const map_type> allGidsMap = rcp(new map_type(TGOT::invalid(), numLocalRangeEntries,
                                                    indexBase, comm));
  // Create vector based on above map.  Populate it with GIDs corresponding to this pid's GIDs in rangeMap.
  mv_type_go allGids(allGidsMap, 1);
  Teuchos::ArrayRCP<GO> allGidsData = allGids.getDataNonConst(0);
 
  for (size_t i = 0; i < numLocalRangeEntries; i++)
    allGidsData[i] = rangeMap->getGlobalElement(i);
  allGidsData = Teuchos::null;
 
  // Create target map that is nontrivial only on pid 0
  GO numTargetMapEntries = TGOT::zero();
  Teuchos::Array<GO> importGidList;
  if (myRank == 0) {
    numTargetMapEntries = rangeMap->getGlobalNumElements();
    importGidList.reserve(numTargetMapEntries);
    for (GO j = 0; j < numTargetMapEntries; ++j) importGidList.push_back(j + indexBase);
  } else {
    importGidList.reserve(numTargetMapEntries);
  }
  RCP<map_type> importGidMap = rcp(new map_type(TGOT::invalid(), importGidList(), indexBase, comm));
 
  // Import all rangeMap GIDs to pid 0
  import_type gidImporter(allGidsMap, importGidMap);
  mv_type_go importedGids(importGidMap, 1);
  importedGids.doImport(allGids, gidImporter, INSERT);
 
  // The following import map will be non-trivial only on pid 0.
  ArrayRCP<const GO> importedGidsData = importedGids.getData(0);
  RCP<const map_type> importMap       = rcp(new map_type(TGOT::invalid(), importedGidsData(), indexBase, comm));
 
  // Importer from original range map to pid 0
  import_type importer(rangeMap, importMap);
  // Target vector on pid 0
  RCP<mv_type> colsOnPid0 = rcp(new mv_type(importMap, numMVs));
 
  RCP<mv_type> ei    = rcp(new mv_type(A.getDomainMap(), numMVs));  // probing vector
  RCP<mv_type> colsA = rcp(new mv_type(A.getRangeMap(), numMVs));   // columns of A revealed by probing
 
  Array<GO> globalColsArray, localColsArray;
  globalColsArray.reserve(numMVs);
  localColsArray.reserve(numMVs);
 
  ArrayRCP<ArrayRCP<Scalar>> eiData(numMVs);
  for (size_t i = 0; i < numMVs; ++i)
    eiData[i] = ei->getDataNonConst(i);
 
  // //////////////////////////////////////
  // Discover A by chunks
  // //////////////////////////////////////
  for (GO k = 0; k < numChunks; ++k) {
    for (size_t j = 0; j < numMVs; ++j) {
      // GO curGlobalCol = maxColGid - numMVs + j + TGOT::one();
      GO curGlobalCol = minColGid + k * numMVs + j;
      globalColsArray.push_back(curGlobalCol);
      // TODO  extract the g2l map outside of this loop loop
      LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
      if (curLocalCol != TLOT::invalid()) {
        eiData[j][curLocalCol] = TGOT::one();
        localColsArray.push_back(curLocalCol);
      }
    }
 
    // drop host views before apply
    for (size_t i = 0; i < numMVs; ++i)
      eiData[i] = Teuchos::null;
    // probe
    A.apply(*ei, *colsA);
 
    colsOnPid0->doImport(*colsA, importer, INSERT);
 
    if (myRank == 0)
      globalNnz += writeColumns(os, *colsOnPid0, numMVs, importedGidsData(),
                                globalColsArray, offsetToUseInPrinting);
 
    // reconstruct dropped eiData
    for (size_t i = 0; i < numMVs; ++i)
      eiData[i] = ei->getDataNonConst(i);
    for (size_t j = 0; j < numMVs; ++j) {
      GO curGlobalCol = minColGid + k * numMVs + j;
      // TODO  extract the g2l map outside of this loop loop
      LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
      if (curLocalCol != TLOT::invalid()) {
        eiData[j][curLocalCol] = TGOT::one();
      }
    }
 
    // zero out the ei's
    for (size_t j = 0; j < numMVs; ++j) {
      for (int i = 0; i < localColsArray.size(); ++i)
        eiData[j][localColsArray[i]] = TGOT::zero();
    }
    globalColsArray.clear();
    localColsArray.clear();
  }
 
  // //////////////////////////////////////
  // Handle leftover part of A
  // //////////////////////////////////////
  if (rem > 0) {
    for (int j = 0; j < rem; ++j) {
      GO curGlobalCol = maxColGid - rem + j + TGOT::one();
      globalColsArray.push_back(curGlobalCol);
      // TODO  extract the g2l map outside of this loop loop
      LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
      if (curLocalCol != TLOT::invalid()) {
        eiData[j][curLocalCol] = TGOT::one();
        localColsArray.push_back(curLocalCol);
      }
    }
 
    // drop host views before apply
    for (size_t i = 0; i < numMVs; ++i)
      eiData[i] = Teuchos::null;
    // probe
    A.apply(*ei, *colsA);
 
    colsOnPid0->doImport(*colsA, importer, INSERT);
    if (myRank == 0)
      globalNnz += writeColumns(os, *colsOnPid0, rem, importedGidsData(),
                                globalColsArray, offsetToUseInPrinting);
 
    // reconstruct dropped eiData
    for (size_t i = 0; i < numMVs; ++i)
      eiData[i] = ei->getDataNonConst(i);
    for (int j = 0; j < rem; ++j) {
      GO curGlobalCol = maxColGid - rem + j + TGOT::one();
      // TODO  extract the g2l map outside of this loop loop
      LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
      if (curLocalCol != TLOT::invalid()) {
        eiData[j][curLocalCol] = TGOT::one();
      }
    }
 
    // zero out the ei's
    for (int j = 0; j < rem; ++j) {
      for (int i = 0; i < localColsArray.size(); ++i)
        eiData[j][localColsArray[i]] = TGOT::zero();
    }
    globalColsArray.clear();
    localColsArray.clear();
  }
 
  // Return the Matrix Market header.  It includes the header
  // line (that starts with "%%"), some comments, and the triple
  // of matrix dimensions and number of nonzero entries.  We
  // don't actually print this here, because we don't know the
  // number of nonzero entries until after probing.
  std::ostringstream oss;
  if (myRank == 0) {
    oss << "%%MatrixMarket matrix coordinate ";
    if (Teuchos::ScalarTraits<typename operator_type::scalar_type>::isComplex) {
      oss << "complex";
    } else {
      oss << "real";
    }
    oss << " general" << std::endl;
    oss << "% Tpetra::Operator" << std::endl;
    std::time_t now = std::time(NULL);
    oss << "% time stamp: " << ctime(&now);
    oss << "% written from " << numProcs << " processes" << std::endl;
    size_t numRows = rangeMap->getGlobalNumElements();
    size_t numCols = domainMap.getGlobalNumElements();
    oss << numRows << " " << numCols << " " << globalNnz << std::endl;
  }
 
  return oss.str();
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_ordinal_type
MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeColumns(std::ostream& os, mv_type const& colsA, size_t const& numCols,
                                                                            Teuchos::ArrayView<const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_ordinal_type> const& rowGids,
                                                                            Teuchos::Array<global_ordinal_type> const& colsArray,
                                                                            typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_ordinal_type const& indexBase) {
  typedef global_ordinal_type GO;
  typedef Teuchos::ScalarTraits<Scalar> STS;
 
  GO nnz               = 0;
  const Scalar zero    = STS::zero();
  const size_t numRows = colsA.getGlobalLength();
  for (size_t j = 0; j < numCols; ++j) {
    Teuchos::ArrayRCP<const Scalar> const curCol = colsA.getData(j);
    const GO J                                   = colsArray[j];
    for (size_t i = 0; i < numRows; ++i) {
      const Scalar val = curCol[i];
      if (val != zero) {
        os << rowGids[i] + indexBase << " " << J + indexBase << " " << val << std::endl;
        ++nnz;
      }
    }
  }
 
  return nnz;
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::writeSparsePerRank(const std::string& filename_prefix,
                                                                                       const std::string& filename_suffix,
                                                                                       const typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sparse_matrix_type& matrix,
                                                                                       const std::string& matrixName,
                                                                                       const std::string& matrixDescription,
                                                                                       const int ranksToWriteAtOnce,
                                                                                       const bool debug) {
  using ST = scalar_type;
  // using LO = local_ordinal_type;
  using GO  = global_ordinal_type;
  using STS = typename Teuchos::ScalarTraits<ST>;
  using Teuchos::RCP;
 
  // Sanity Checks
  trcp_tcomm_t comm = matrix.getComm();
  TEUCHOS_TEST_FOR_EXCEPTION(comm.is_null(), std::invalid_argument,
                             "The input matrix's communicator (Teuchos::Comm object) is null.");
  TEUCHOS_TEST_FOR_EXCEPTION(matrix.isGloballyIndexed() || !matrix.isFillComplete(), std::invalid_argument,
                             "The input matrix must not be GloballyIndexed and must be fillComplete.");
 
  // Setup
  const int myRank           = comm->getRank();
  const int numProc          = comm->getSize();
  std::string filename       = filename_prefix + std::to_string(myRank) + filename_suffix;
  RCP<const map_type> rowMap = matrix.getRowMap();
  RCP<const map_type> colMap = matrix.getColMap();
  size_t local_nnz           = matrix.getLocalNumEntries();
  size_t local_num_rows      = rowMap->getLocalNumElements();
  size_t local_num_cols      = colMap->getLocalNumElements();
  const GO rowIndexBase      = rowMap->getIndexBase();
  const GO colIndexBase      = colMap->getIndexBase();
 
  // Bounds check the writing limits
  int rank_limit = std::min(std::max(ranksToWriteAtOnce, 1), numProc);
 
  // Start the writing
  for (int base_rank = 0; base_rank < numProc; base_rank += rank_limit) {
    int stop = std::min(base_rank + rank_limit, numProc);
 
    if (base_rank <= myRank && myRank < stop) {
      // My turn to write
      std::ofstream out(filename);
 
      // MatrixMarket Header
      out << "%%MatrixMarket matrix coordinate "
          << (STS::isComplex ? "complex" : "real")
          << " general" << std::endl;
 
      // Print comments (the matrix name and / or description).
      if (matrixName != "") {
        printAsComment(out, matrixName);
      }
      if (matrixDescription != "") {
        printAsComment(out, matrixDescription);
      }
 
      // Print the Matrix Market header (# local rows, # local columns, #
      // local enonzeros).  This will *not* be read correctly by a generic matrix
      // market reader since we'll be writing out GIDs here and local row/col counts
      out << local_num_rows << " " << local_num_cols << " " << local_nnz << std::endl;
 
      {
        // Make the output stream write floating-point numbers in
        // scientific notation.  It will politely put the output
        // stream back to its state on input, when this scope
        // terminates.
        Teuchos::SetScientific<ST> sci(out);
 
        for (size_t l_row = 0; l_row < local_num_rows; l_row++) {
          GO g_row = rowMap->getGlobalElement(l_row);
 
          typename sparse_matrix_type::local_inds_host_view_type indices;
          typename sparse_matrix_type::values_host_view_type values;
          matrix.getLocalRowView(l_row, indices, values);
          for (size_t ii = 0; ii < indices.extent(0); ii++) {
            const GO g_col = colMap->getGlobalElement(indices(ii));
            // Convert row and column indices to 1-based.
            // This works because the global index type is signed.
            out << (g_row + 1 - rowIndexBase) << " "
                << (g_col + 1 - colIndexBase) << " ";
            if (STS::isComplex) {
              out << STS::real(values(ii)) << " " << STS::imag(values(ii));
            } else {
              out << values(ii);
            }
            out << std::endl;
          }  // For each entry in the current row
        }    // For each row of the matrix
      }      // end Teuchos::SetScientfic scoping
 
      out.close();
    }  // end if base_rank <= myRank < stop
 
    // Barrier after each writing "batch" to make sure we're not hammering the file system
    // too aggressively
    comm->barrier();
 
  }  // end outer loop
 
}  // end writeSparsePerRank
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
template <typename T>
typename MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::trcp_tcomm_t MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getComm(const Teuchos::RCP<T>& obj) {
  return obj.is_null() ? Teuchos::null : obj->getComm();
}
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
int MatrixMarketWriter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getRank(const trcp_tcomm_t& comm) {
  return comm.is_null() ? 0 : comm->getRank();
}
 
}  // namespace Tpetra
 
namespace Tpetra {
#define TPETRA_MATRIXMARKET_INSTANT(SCALAR, LO, GO, NODE)  \
  template class MatrixMarketReader<SCALAR, LO, GO, NODE>; \
  template class MatrixMarketWriter<SCALAR, LO, GO, NODE>;
}  // namespace Tpetra
 
#endif