MatrixMarket_TpetraNew.hpp

#ifndef __MATRIXMARKET_TPETRA_NEW_HPP
#define __MATRIXMARKET_TPETRA_NEW_HPP

// Functions to read a MatrixMarket file and load it into a matrix.
// Adapted from 
//    Trilinos/packages/epetraext/src/inout/EpetraExt_CrsMatrixIn.cpp
// Modified by Jon Berry and Karen Devine to make matrix reallocations 
// more efficient; rather than insert each non-zero separately, we 
// collect rows of non-zeros for insertion.
// Modified by Karen Devine and Steve Plimpton to prevent all processors 
// from reading the same file at the same time (ouch!); chunks of the file 
// are read and broadcast by processor zero; each processor then extracts 
// its nonzeros from the chunk, sorts them by row, and inserts them into 
// the matrix.
// The variable "chunkSize" can be changed to modify the size of the chunks
// read from the file.  Larger values of chunkSize lead to faster execution
// and greater memory use.

private:

template <typename global_ordinal_type, typename scalar_type>
static 
Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > 
buildDistribution(
  std::string &distribution,  // string indicating whether to use 1D, 2D or
                              // file-based matrix distribution
  size_t nRow,                // Number of global matrix rows
  size_t nCol,                // Number of global matrix rows
  const Teuchos::ParameterList &params, // Parameters to the file reading
  const Teuchos::RCP<const Teuchos::Comm<int> > &comm
                            // communicator to be used in maps
)
{
  // Function to build the sets of GIDs for 1D or 2D matrix distribution.
  // Output is a Distribution object.

  int me = comm->getRank();

  using basedist_t = Distribution<global_ordinal_type,scalar_type>;
  basedist_t *retval;

  bool randomize = false;   // Randomly permute rows and columns before storing
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("randomize");
  if (pe != NULL) 
    randomize = pe->getValue<bool>(&randomize);
  }

  std::string partitionFile = "";  // File to reassign rows to parts
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("partitionFile");
  if (pe != NULL) 
    partitionFile = pe->getValue<std::string>(&partitionFile);
  }

  if (distribution == "2D") {   // Generate 2D distribution
    if (partitionFile != "") {
      // Generate 2D distribution with vector partition file
      TEUCHOS_TEST_FOR_EXCEPTION(randomize, std::logic_error,
              "Randomization not available for 2DVec distributions.");
      Distribution2DVec<global_ordinal_type,scalar_type> *dist =
              new Distribution2DVec<global_ordinal_type, scalar_type>(
                                    nRow, comm, params, partitionFile);
      retval = dynamic_cast<basedist_t *>(dist);
    }
    else if (randomize) {
      // Randomly assign rows/columns to processors
      Distribution2DRandom<global_ordinal_type,scalar_type> *dist =
              new Distribution2DRandom<global_ordinal_type,scalar_type>(
                                       nRow, comm, params);
      retval = dynamic_cast<basedist_t *>(dist);
    }
    else {
      // The vector will be distributed linearly, with extra vector entries
      // spread among processors to maintain balance in rows and columns.
      Distribution2DLinear<global_ordinal_type,scalar_type> *dist =
              new Distribution2DLinear<global_ordinal_type,scalar_type>(
                                       nRow, comm, params);
      retval = dynamic_cast<basedist_t *>(dist);
    }
  }
  else if (distribution == "1D") {
    if (partitionFile != "") {
      // Generate 1D row-based distribution with vector partition file
      Distribution1DVec<global_ordinal_type,scalar_type> *dist =
              new Distribution1DVec<global_ordinal_type,scalar_type>(
                                    nRow, comm, params, partitionFile);
      retval = dynamic_cast<basedist_t*>(dist);
    }
    else if (randomize) {
      // Randomly assign rows to processors.
      Distribution1DRandom<global_ordinal_type,scalar_type> *dist =
              new Distribution1DRandom<global_ordinal_type,scalar_type>(
                                       nRow, comm, params);
      retval = dynamic_cast<basedist_t *>(dist);
    }
    else {
      // Linear map similar to Trilinos default.
      Distribution1DLinear<global_ordinal_type,scalar_type> *dist =
              new Distribution1DLinear<global_ordinal_type,scalar_type>(
                                       nRow, comm, params);
      retval = dynamic_cast<basedist_t *>(dist);
    }
  }
  else if (distribution == "LowerTriangularBlock") {
    if (randomize && me == 0) {
      std::cout << "WARNING: Randomization not available for "
                << "LowerTriangularBlock distributions." << std::endl;
    }

    DistributionLowerTriangularBlock<global_ordinal_type,scalar_type> *dist = 
          new DistributionLowerTriangularBlock<global_ordinal_type,scalar_type>(
                                       nRow, comm, params);
    retval = dynamic_cast<basedist_t*>(dist);
  }
  else if (distribution == "MMFile") {
    // Generate user-defined distribution from Matrix-Market file
    if (randomize && me == 0) {
      std::cout << "WARNING: Randomization not available for MMFile "
                << "distributions." << std::endl;
    }
    DistributionMMFile<global_ordinal_type,scalar_type> *dist =
            new DistributionMMFile<global_ordinal_type,scalar_type>(
                                   nRow, comm, params);
    retval = dynamic_cast<basedist_t*>(dist);
  }

  else {
    std::cout << "ERROR:  Invalid distribution method " << distribution 
              << std::endl;
    exit(-1);
  }
  return Teuchos::rcp<basedist_t>(retval);
}

static
void 
readMatrixMarket(  
  const std::string &filename,    // MatrixMarket file to read
  const Teuchos::RCP<const Teuchos::Comm<int> > &comm,  
  const Teuchos::ParameterList &params,
  size_t &nRow,
  size_t &nCol,
  typename Distribution<global_ordinal_type,scalar_type>::LocalNZmap_t &localNZ,
  Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > &dist
)
{
  bool useTimers = false;   // should we time various parts of the reader?
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("useTimers");
  if (pe != NULL) 
    useTimers = pe->getValue<bool>(&useTimers);
  }

  Teuchos::RCP<Teuchos::TimeMonitor> timer = Teuchos::null;
  if (useTimers) 
    timer = rcp(new Teuchos::TimeMonitor(
                    *Teuchos::TimeMonitor::getNewTimer("RMM parameterSetup")));

  int me = comm->getRank();

  bool verbose = false;   // Print status as reading
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("verbose");
  if (pe != NULL) 
    verbose = pe->getValue<bool>(&verbose);
  }

  size_t chunkSize = 500000;   // Number of lines to read / broadcast at once
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("chunkSize");
  if (pe != NULL) 
    chunkSize = pe->getValue<size_t>(&chunkSize);
  }

  bool symmetrize = false;  // Symmetrize the matrix
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("symmetrize");
  if (pe != NULL) 
    symmetrize = pe->getValue<bool>(&symmetrize);
  }

  bool transpose = false;   // Store the transpose
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("transpose");
  if (pe != NULL) 
    transpose = pe->getValue<bool>(&transpose);
  }

  std::string diagonal = "";  // Are diagonal entries required (so add them)
                              // or ignored (so delete them) or neither?
                              // Default is neither.
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("diagonal");
  if (pe != NULL) 
    diagonal = pe->getValue<std::string>(&diagonal);
  }
  bool ignoreDiagonal = (diagonal == "exclude");
  bool requireDiagonal = (diagonal == "require");

  std::string distribution = "1D";  // Default distribution is 1D row-based
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("distribution");
  if (pe != NULL) 
    distribution = pe->getValue<std::string>(&distribution);
  }

  if (useTimers) {
    timer = Teuchos::null;
    timer = rcp(new Teuchos::TimeMonitor(
                    *Teuchos::TimeMonitor::getNewTimer("RMM header")));
  }

  if (verbose && me == 0)
    std::cout << "Reading matrix market file... " << filename << std::endl;

  FILE *fp = NULL;
  size_t dim[3] = {0,0,0};  // nRow, nCol, nNz as read from MatrixMarket
  MM_typecode mmcode;

  if (me == 0) {

    fp = fopen(filename.c_str(), "r");

    if (fp == NULL) {
      std::cout << "Error:  cannot open file " << filename << std::endl;
    }
    else {
      // Read MatrixMarket banner to get type of matrix
      if (mm_read_banner(fp, &mmcode) != 0) {
        std::cout << "Error:  bad MatrixMarket banner " << std::endl;
      }
      else {
        // Error check for file types that this function can read
        if (!mm_is_valid(mmcode) ||
             mm_is_dense(mmcode) || mm_is_array(mmcode) ||
             mm_is_complex(mmcode) || 
             mm_is_skew(mmcode) || mm_is_hermitian(mmcode)) {
          std::cout << "Error:  matrix type is not supported by this reader\n "
                    << "This reader supports sparse, coordinate, "
                    << "real, pattern, integer, symmetric, general"
                    << std::endl;
        }
        else {
          // Read nRow, nCol, nNz from MatrixMarket file
          if (mm_read_mtx_crd_size(fp, &dim[0], &dim[1], &dim[2]) != 0) {
            std::cout << "Error:  bad matrix dimensions " << std::endl;
            dim[0] = dim[1] = dim[2] = 0;
          }
        }
      }
    }
  }

  // Broadcast relevant info
  // Bad input if dim[0] or dim[1] still is zero after broadcast; 
  // all procs throw together
  Teuchos::broadcast<int, size_t>(*comm, 0, 3, dim);
  if (dim[0] == 0 || dim[1] == 0) {
    throw std::runtime_error("Error:  bad matrix header information");
  }
  Teuchos::broadcast<int, char>(*comm, 0, sizeof(MM_typecode), mmcode);

  nRow = dim[0];
  nCol = dim[1];

  TEUCHOS_TEST_FOR_EXCEPTION(nRow != nCol, std::logic_error,
        "This overload of readSparseFile requires nRow == nCol " 
         << "(nRow = " << nRow << ", nCol = " << nCol << "); "
         << "For now, use a different overload of readSparseFile until this "
         << "overload is fixed to read rectangular matrices. "
         << "See Trilinos github issues #7045 and #8472.");
  
  size_t nNz = dim[2];
  bool patternInput = mm_is_pattern(mmcode);
  bool symmetricInput = mm_is_symmetric(mmcode);
  if (symmetricInput) symmetrize = true;

  if (verbose && me == 0)
    std::cout << "Matrix market file... " 
              << "\n  symmetrize          = " << symmetrize
              << "\n  transpose           = " << transpose
              << "\n  change diagonal     = " << diagonal
              << "\n  distribution        = " << distribution
              << std::endl;

  if (useTimers) {
    timer = Teuchos::null;
    timer = rcp(new Teuchos::TimeMonitor(
                    *Teuchos::TimeMonitor::getNewTimer("RMM distribution")));
  }
  
  // Create distribution based on nRow, nCol, npRow, npCol
  dist = buildDistribution<global_ordinal_type,scalar_type>(distribution,
                  nRow, nCol, params,
                  comm);
  if (useTimers) {
    timer = Teuchos::null;
    timer = rcp(new Teuchos::TimeMonitor(
                    *Teuchos::TimeMonitor::getNewTimer("RMM readChunks")));
  }

  std::set<global_ordinal_type> diagset;  
                            // If diagonal == require, this set keeps track of 
                            // which diagonal entries already existing so we can
                            // add those that don't

  using nzindex_t = 
        typename Distribution<global_ordinal_type,scalar_type>::NZindex_t;

  // Chunk information and buffers
  const int maxLineLength = 81;
  char *buffer = new char[chunkSize*maxLineLength+1];
  size_t nChunk; 
  size_t nMillion = 0;
  size_t nRead = 0;
  size_t rlen;

  auto timerRead = Teuchos::TimeMonitor::getNewTimer("RMM    readNonzeros");
  auto timerSelect = Teuchos::TimeMonitor::getNewTimer("RMM    selectNonzeros");
  // Read chunks until the entire file is read
  Teuchos::RCP<Teuchos::TimeMonitor> innerTimer = Teuchos::null;
  while (nRead < nNz) {

    innerTimer = rcp(new Teuchos::TimeMonitor(*timerRead));

    if (nNz-nRead > chunkSize) nChunk = chunkSize;
    else nChunk = (nNz - nRead);

    // Processor 0 reads a chunk
    if (me == 0) {
      char *eof;
      rlen = 0;
      for (size_t i = 0; i < nChunk; i++) {
        eof = fgets(&buffer[rlen],maxLineLength,fp);
        if (eof == NULL) {
          std::cout << "Unexpected end of matrix file." << std::endl;
          std::cout.flush();
          exit(-1);
        }
        rlen += strlen(&buffer[rlen]);
      }
      buffer[rlen++]='\n';
    }

    // Processor 0 broadcasts the chunk
    Teuchos::broadcast<int, size_t>(*comm, 0, 1, &rlen);
    Teuchos::broadcast<int, char>(*comm, 0, rlen, buffer);

    buffer[rlen++] = '\0';
    nRead += nChunk;

    innerTimer = Teuchos::null;
    innerTimer = rcp(new Teuchos::TimeMonitor(*timerSelect));

    // All processors check the received data, saving nonzeros belonging to them
    char *lineptr = buffer;
    for (rlen = 0; rlen < nChunk; rlen++) {

      char *next = strchr(lineptr,'\n');
      global_ordinal_type I = atol(strtok(lineptr," \t\n")) 
              - 1; // since MatrixMarket files are one-based
      global_ordinal_type J = atol(strtok(NULL," \t\n")) 
              - 1; // since MatrixMarket files are one-based

      scalar_type V = (patternInput ? -1. : atof(strtok(NULL," \t\n")));
      lineptr = next + 1;

      // Special processing of nonzero
      if ((I == J) && ignoreDiagonal) continue;

      if (transpose) std::swap<global_ordinal_type>(I,J);

      // Add nonzero (I,J) to the map if it should be on this processor
      // Some file-based distributions have processor assignment stored as 
      // the non-zero's value, so pass the value to Mine.
      if (dist->Mine(I,J,int(V))) {
        nzindex_t idx = std::make_pair(I,J);
        localNZ[idx] = V;   
        if (requireDiagonal && (I == J)) diagset.insert(I);
      }

      // If symmetrizing, add (J,I) to the map if it should be on this processor
      // Some file-based distributions have processor assignment stored as 
      // the non-zero's value, so pass the value to Mine.
      if (symmetrize && (I != J) && dist->Mine(J,I,int(V))) {
        //  Add entry (J, I) if need to symmetrize
        //  This processor keeps this non-zero.
        nzindex_t idx = std::make_pair(J,I);
        localNZ[idx] = V;   
      }
    }

    // Status check
    if (verbose) {
      if (nRead / 1000000 > nMillion) {
        nMillion++;
        if (me == 0) std::cout << nMillion << "M ";
      }
    }

    innerTimer = Teuchos::null;
  }

  if (verbose && me == 0) 
    std::cout << std::endl << nRead << " nonzeros read " << std::endl;

  if (fp != NULL) fclose(fp);
  delete [] buffer;

  if (useTimers) {
    timer = Teuchos::null;
    timer = rcp(new Teuchos::TimeMonitor(
                    *Teuchos::TimeMonitor::getNewTimer("RMM diagonal")));
  }

  //  Add diagonal entries if they are required.
  //  Check to make sure they are all here; add them if they are not.
  if (requireDiagonal) {
    if (dist->DistType() == MMFile) {
      // All diagonal entries should be present in the file; we cannot create
      // them for file-based data distributions.  
      // Do an error check to ensure they are all there.
      size_t localss = diagset.size();
      size_t globalss;
      Teuchos::reduceAll<int, size_t>(*comm, Teuchos::REDUCE_SUM, 1,
                                      &localss, &globalss);
      TEUCHOS_TEST_FOR_EXCEPTION(globalss != nRow, std::logic_error,
        "File-based nonzero distribution requires all diagonal "
         << "entries to be present in the file.  \n" 
         << "Number of diagonal entries in file = " << globalss << "\n"
         << "Number of matrix rows = " << nRow << "\n");
    }
    else {
      for (global_ordinal_type i = 0;
                               i < static_cast<global_ordinal_type>(nRow); i++) 
      {
        if (dist->Mine(i,i)) {
          if (diagset.find(i) == diagset.end()) {
            nzindex_t idx = std::make_pair(i,i);
            localNZ[idx] = 0;   
          }
        }
      }
    }
  }
  // Done with diagset; free its memory
  std::set<global_ordinal_type>().swap(diagset);

  if (useTimers) 
    timer = Teuchos::null;
}

//        (This format is used by the upcoming readBinary function)         //
// 
// File format:
//   #vertices #edges src_1 dst_1 src_2 dst_2 ... src_#edges dst_#edges
// 
// Types:
//   #edges:            unsigned long long
//   everything else:   unsigned int
//
// Base of indexing: 
//   1
// 
//
// A code example that converts MatrixMarket format into this format:
//
//   typedef unsigned int ord_type;
//   typedef unsigned long long size_type;
//
//   std::string line;
//   std::ifstream in(infilename);
//   std::ofstream out(outfilename,  std::ios::out | std::ios::binary);
//
//  ord_type nv, dummy, edge[2];
//  size_type ne;
//
//   do
//     std::getline (in, line);
//   while(line[0] == '%');
//
//   std::stringstream ss(line);
//   ss >> nv >> dummy >> ne;
//   out.write((char *)&nv, sizeof(ord_type));
//   out.write((char *)&ne, sizeof(size_type));
//
//   while(std::getline(in, line)) {
//     std::stringstream ss2(line);
//     ss2 >> edge[0] >> edge[1];
//     out.write((char *)edge, sizeof(ord_type)*2);
//
//    }
//   in.close();
//   out.close();
//
  
static
void 
readBinary(  
  const std::string &filename,    // MatrixMarket file to read
  const Teuchos::RCP<const Teuchos::Comm<int> > &comm,  
  const Teuchos::ParameterList &params,
  size_t &nRow,
  size_t &nCol,
  typename Distribution<global_ordinal_type,scalar_type>::LocalNZmap_t &localNZ,
  Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > &dist
)
{

  int me = comm->getRank();

  bool verbose = false;   // Print status as reading
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("verbose");
  if (pe != NULL) 
    verbose = pe->getValue<bool>(&verbose);
  }

  size_t chunkSize = 500000;   // Number of lines to read / broadcast at once
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("chunkSize");
  if (pe != NULL) 
    chunkSize = pe->getValue<size_t>(&chunkSize);
  }

  bool symmetrize = false;  // Symmetrize the matrix
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("symmetrize");
  if (pe != NULL) 
    symmetrize = pe->getValue<bool>(&symmetrize);
  }

  bool transpose = false;   // Store the transpose
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("transpose");
  if (pe != NULL) 
    transpose = pe->getValue<bool>(&transpose);
  }

  std::string diagonal = "";  // Are diagonal entries required (so add them)
                              // or ignored (so delete them) or neither?
                              // Default is neither.
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("diagonal");
  if (pe != NULL) 
    diagonal = pe->getValue<std::string>(&diagonal);
  }
  bool ignoreDiagonal = (diagonal == "exclude");
  bool requireDiagonal = (diagonal == "require");

  std::string distribution = "1D";  // Default distribution is 1D row-based
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("distribution");
  if (pe != NULL) 
    distribution = pe->getValue<std::string>(&distribution);
  }

  if (verbose && me == 0)
    std::cout << "Reading binary file... " << filename << std::endl;

  FILE *fp = NULL;
  size_t dim[3] = {0,0,0};  // Expected to read nRow and nNz, nCol = nRow

  if (me == 0) {

    fp = fopen(filename.c_str(), "rb");

    if (fp == NULL) {
      std::cout << "Error:  cannot open file " << filename << std::endl;
    }
    else {
      // The header in the binary file contains only numVertices and numEdges
      unsigned int nv = 0;
      unsigned long long ne = 0;
      if (fread(&nv, sizeof(unsigned int), 1, fp) != 1)
           throw std::runtime_error("Error: bad number of read values.");
      if (fread(&ne, sizeof(unsigned long long), 1, fp) != 1)
           throw std::runtime_error("Error: bad number of read values.");
      dim[0] = nv;
      dim[1] = dim[0];
      dim[2] = ne;
    }
  }

  // Broadcast relevant info
  // Bad input if dim[0] or dim[1] still is zero after broadcast; 
  // all procs throw together
  Teuchos::broadcast<int, size_t>(*comm, 0, 3, dim);
  if (dim[0] == 0 || dim[1] == 0) {
    throw std::runtime_error("Error:  bad matrix header information");
  }

  nRow = dim[0];
  nCol = dim[1];
  size_t nNz = dim[2];

  if (verbose && me == 0)
    std::cout << "Binary file... " 
              << "\n  symmetrize          = " << symmetrize
              << "\n  transpose           = " << transpose
              << "\n  change diagonal     = " << diagonal
              << "\n  distribution        = " << distribution
              << std::endl;
  
  // Create distribution based on nRow, nCol, npRow, npCol
  dist = buildDistribution<global_ordinal_type,scalar_type>(distribution,
                  nRow, nCol, params,
                  comm);

  std::set<global_ordinal_type> diagset;  
                            // If diagonal == require, this set keeps track of 
                            // which diagonal entries already existing so we can
                            // add those that don't

  using nzindex_t = 
        typename Distribution<global_ordinal_type,scalar_type>::NZindex_t;

  // Chunk information and buffers
  unsigned int *buffer = new unsigned int[chunkSize*2];
  size_t nChunk; 
  size_t nMillion = 0;
  size_t nRead = 0;
  size_t rlen;
  const scalar_type ONE = Teuchos::ScalarTraits<scalar_type>::one();


  // Read chunks until the entire file is read
  while (nRead < nNz) {
    if (nNz-nRead > chunkSize) nChunk = chunkSize;
    else nChunk = (nNz - nRead);

    // Processor 0 reads a chunk
    if (me == 0) {
      size_t ret = 0;
      rlen = 0;
      for (size_t i = 0; i < nChunk; i++) {
  ret = fread(&buffer[rlen], sizeof(unsigned int), 2, fp);
        if (ret == 0) {
          std::cout << "Unexpected end of matrix file." << std::endl;
          std::cout.flush();
          exit(-1);
        }
        rlen += 2;
      }
    }

    // Processor 0 broadcasts the chunk
    Teuchos::broadcast<int, size_t>(*comm, 0, 1, &rlen);
    Teuchos::broadcast<int, unsigned int>(*comm, 0, rlen, buffer);

    nRead += nChunk;

    // All processors check the received data, saving nonzeros belonging to them
    for (rlen = 0; rlen < nChunk; rlen++) {

      global_ordinal_type I = buffer[2*rlen]-1;
      global_ordinal_type J = buffer[2*rlen+1]-1;
      
      // Special processing of nonzero
      if ((I == J) && ignoreDiagonal) continue;

      if (transpose) std::swap<global_ordinal_type>(I,J);

      // Add nonzero (I,J) to the map if it should be on this processor
      // Some file-based distributions have processor assignment stored as 
      // the non-zero's value, so pass the value to Mine.
      if (dist->Mine(I,J,ONE)) {
        nzindex_t idx = std::make_pair(I,J);
        localNZ[idx] = ONE;  // For now, the input binary format does not 
                       // support numeric values, so we insert one.  
        if (requireDiagonal && (I == J)) diagset.insert(I);
      }

      // If symmetrizing, add (J,I) to the map if it should be on this processor
      // Some file-based distributions have processor assignment stored as 
      // the non-zero's value, so pass the value to Mine.
      if (symmetrize && (I != J) && dist->Mine(J,I,ONE)) {
        //  Add entry (J, I) if need to symmetrize
        //  This processor keeps this non-zero.
        nzindex_t idx = std::make_pair(J,I);
        localNZ[idx] = ONE;   
      }
    }

    // Status check
    if (verbose) {
      if (nRead / 1000000 > nMillion) {
        nMillion++;
        if (me == 0) std::cout << nMillion << "M ";
      }
    }
  }

  if (verbose && me == 0) 
    std::cout << std::endl << nRead << " nonzeros read " << std::endl;

  if (fp != NULL) fclose(fp);
  delete [] buffer;

  //  Add diagonal entries if they are required.
  //  Check to make sure they are all here; add them if they are not.
  if (requireDiagonal) {
    if (dist->DistType() == MMFile) {
      // All diagonal entries should be present in the file; we cannot create
      // them for file-based data distributions.  
      // Do an error check to ensure they are all there.
      size_t localss = diagset.size();
      size_t globalss;
      Teuchos::reduceAll<int, size_t>(*comm, Teuchos::REDUCE_SUM, 1,
                                      &localss, &globalss);
      TEUCHOS_TEST_FOR_EXCEPTION(globalss != nRow, std::logic_error,
        "File-based nonzero distribution requires all diagonal "
         << "entries to be present in the file.  \n" 
         << "Number of diagonal entries in file = " << globalss << "\n"
         << "Number of matrix rows = " << nRow << "\n");
    }
    else {
      for (global_ordinal_type i = 0;
                               i < static_cast<global_ordinal_type>(nRow); i++)
      {
        if (dist->Mine(i,i)) {
          if (diagset.find(i) == diagset.end()) {
            nzindex_t idx = std::make_pair(i,i);
            localNZ[idx] = 0;   
          }
        }
      }
    }
  }
  // Done with diagset; free its memory
  std::set<global_ordinal_type>().swap(diagset);

}

//   (This format is used by the upcoming readPerProcessBinary function)    //
// 
// 
// FILE FORMAT:
//     globalNumRows globalNumCols localNumNnzs row_1 col_1 ... row_n col_n, 
//  where n = localNumNnzs
// 
//
// ASSUMPTIONS:
//  - The nonzeros should be sorted by their row indices within each file.
//  - Nonzeros have global row and column indices.
//  - The user specifies the base <filename>, where rank i reads <filename>.i.cooBin.  
// 
// 
// TYPES:
//  localNumNnzs:      unsigned long long
//  everything else:   unsigned int
//
//
// BASE OF INDEXING: 1
//
// 
static
void 
readPerProcessBinary(
  const std::string &filename,
  const Teuchos::RCP<const Teuchos::Comm<int> > &comm,  
  const Teuchos::ParameterList &params,
  size_t &nRow,
  size_t &nCol,
  typename Distribution<global_ordinal_type,scalar_type>::LocalNZmap_t &localNZ,
  Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > &dist,
  unsigned int* &buffer,
  size_t &nNz
)
{

  int me = comm->getRank();

  bool verbose = false;   // Print status as reading
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("verbose");
  if (pe != NULL) 
    verbose = pe->getValue<bool>(&verbose);
  }

  std::string distribution = "1D";  // Default distribution is 1D row-based
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("distribution");
  if (pe != NULL) 
    distribution = pe->getValue<std::string>(&distribution);
  }

  if (verbose && me == 0)
    std::cout << "Reading per-process binary files... " << filename << std::endl;


  std::string rankFileName = filename + "." + std::to_string(me) + ".cooBin";
  FILE *fp = NULL;

  fp = fopen(rankFileName.c_str(), "rb");
  if (fp == NULL) {
    std::cout << "Error:  cannot open file " << filename << std::endl;
    throw std::runtime_error("Error:  non-existing input file: " + rankFileName);
  }
 
  // The header in each per-process file:  globalNumRows globalNumCols localNumNonzeros
  unsigned int globalNumRows = 0, globalNumCols = 0;
  unsigned long long localNumNonzeros = 0;
  if (fread(&globalNumRows, sizeof(unsigned int), 1, fp) != 1)
    throw std::runtime_error("Error: bad number of read values.");
  if (fread(&globalNumCols, sizeof(unsigned int), 1, fp) != 1)
    throw std::runtime_error("Error: bad number of read values.");
  if (fread(&localNumNonzeros, sizeof(unsigned long long), 1, fp) != 1)
    throw std::runtime_error("Error: bad number of read values.");
 
  nRow = static_cast<size_t>(globalNumRows);
  nCol = static_cast<size_t>(globalNumCols);
  nNz = static_cast<size_t>(localNumNonzeros);

  // Fill the simple buffer array instead of a std::map
  // S. Acer: With large graphs, we can't afford std::map
  buffer = new unsigned int[nNz*2];

  if(nNz > 0) {
    size_t ret = fread(buffer, sizeof(unsigned int), 2*nNz, fp);
    if (ret == 0) {
      std::cout << "Unexpected end of matrix file: " << rankFileName << std::endl;
      std::cout.flush();
      delete [] buffer;
      exit(-1);
    }
  }
  if (fp != NULL) fclose(fp);

  // This barrier is not necessary but useful for debugging
  comm->barrier();
  if(verbose && me == 0)
    std::cout << "All ranks finished reading their nonzeros from their individual files\n";

  // Create distribution based on nRow, nCol, npRow, npCol
  dist = buildDistribution<global_ordinal_type,scalar_type>(distribution,
                  nRow, nCol, params,
                  comm);
}

public:

// This is the default interface.
static Teuchos::RCP<sparse_matrix_type>
readSparseFile(
  const std::string &filename,    // MatrixMarket file to read
  const Teuchos::RCP<const Teuchos::Comm<int> > &comm,  
  const Teuchos::ParameterList &params
)
{
  Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > dist;
  return readSparseFile(filename, comm, params, dist);
}

// This version has the Distribution object as an output parameter.
// S. Acer needs the distribution object to get the chunk cuts from
// LowerTriangularBlock distribution.  
static Teuchos::RCP<sparse_matrix_type>
readSparseFile(
  const std::string &filename,    // MatrixMarket file to read
  const Teuchos::RCP<const Teuchos::Comm<int> > &comm,  
  const Teuchos::ParameterList &params,
  Teuchos::RCP<Distribution<global_ordinal_type,scalar_type> > &dist 
)
{
  bool useTimers = false;   // should we time various parts of the reader?
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("useTimers");
  if (pe != NULL) 
    useTimers = pe->getValue<bool>(&useTimers);
  }

  Teuchos::RCP<Teuchos::TimeMonitor> timer = Teuchos::null;
  if (useTimers) 
    timer = rcp(new Teuchos::TimeMonitor(
                    *Teuchos::TimeMonitor::getNewTimer("RSF parameterSetup")));

  int me = comm->getRank();

  // Check parameters to determine how to process the matrix while reading
  // TODO:  Add validators for the parameters

  bool verbose = false;   // Print status as reading
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("verbose");
  if (pe != NULL) 
    verbose = pe->getValue<bool>(&verbose);
  }
  
  bool callFillComplete = true;   // should we fillComplete the new CrsMatrix?
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("callFillComplete");
  if (pe != NULL) 
    callFillComplete = pe->getValue<bool>(&callFillComplete);
  }

  // Don't want to use MatrixMarket's coordinate reader, because don't want
  // entire matrix on one processor.
  // Instead, Proc 0 reads nonzeros in chunks and broadcasts chunks to all 
  // processors.
  // All processors insert nonzeros they own into a std::map

  // Storage for this processor's nonzeros.
  using localNZmap_t = 
        typename Distribution<global_ordinal_type,scalar_type>::LocalNZmap_t;
  localNZmap_t localNZ;

  bool binary = false;   // should we read a binary file?
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("binary");
  if (pe != NULL) 
    binary = pe->getValue<bool>(&binary);
  }

  bool readPerProcess = false;   // should we read a separate file per process?
  {
  const Teuchos::ParameterEntry *pe = params.getEntryPtr("readPerProcess");
  if (pe != NULL) 
    readPerProcess = pe->getValue<bool>(&readPerProcess);
  }

  if (useTimers) {
    const char *timername = (binary?"RSF readBinary":"RSF readMatrixMarket");
    timer = Teuchos::null;
    timer = rcp(new Teuchos::TimeMonitor(
                   *Teuchos::TimeMonitor::getNewTimer(timername)));
  }
  
  // Read nonzeros from the given file(s)
  size_t nRow = 0, nCol = 0;
  unsigned int *buffer=0; size_t nNz = 0;
  if(binary){
    if(readPerProcess)
      readPerProcessBinary(filename, comm, params, nRow, nCol, localNZ, dist, buffer, nNz);
    else 
      readBinary(filename, comm, params, nRow, nCol, localNZ, dist);
  }
  else
    readMatrixMarket(filename, comm, params, nRow, nCol, localNZ, dist);

  if(readPerProcess == false){

    // Redistribute nonzeros as needed to satisfy the Distribution
    // For most Distributions, this is a no-op    
    if (useTimers) {
      timer = Teuchos::null;
      timer = rcp(new Teuchos::TimeMonitor(
                   *Teuchos::TimeMonitor::getNewTimer("RSF redistribute")));
    }
  
    dist->Redistribute(localNZ);
  }

  if (useTimers) {
    timer = Teuchos::null;
    timer = rcp(new Teuchos::TimeMonitor(
                   *Teuchos::TimeMonitor::getNewTimer("RSF nonzerosConstruction")));
  }

  //  Now construct the matrix.
  //  Count number of entries in each row for best use of StaticProfile
  if (verbose && me == 0)
    std::cout << std::endl << "Constructing the matrix" << std::endl;

  Teuchos::Array<global_ordinal_type> rowIdx;
  Teuchos::Array<size_t> nnzPerRow;
  Teuchos::Array<global_ordinal_type> colIdx;
  Teuchos::Array<scalar_type> val;
  Teuchos::Array<size_t> offsets;

  if(readPerProcess) {
    global_ordinal_type prevI = std::numeric_limits<global_ordinal_type>::max();
    for (size_t it = 0; it < nNz; it++){
      global_ordinal_type I = buffer[2*it]-1;
      global_ordinal_type J = buffer[2*it+1]-1;
      if (prevI != I) {
  prevI = I;
  rowIdx.push_back(I);
  nnzPerRow.push_back(0);
      }
      nnzPerRow.back()++;
      colIdx.push_back(J);
    }
    delete [] buffer;      
  }
  else {
    // Exploit fact that map has entries sorted by I, then J
    global_ordinal_type prevI = std::numeric_limits<global_ordinal_type>::max();
    for (auto it = localNZ.begin(); it != localNZ.end(); it++) {
      global_ordinal_type I = it->first.first;
      global_ordinal_type J = it->first.second;
      scalar_type V = it->second;
      if (prevI != I) {
  prevI = I;
  rowIdx.push_back(I);
  nnzPerRow.push_back(0);
      }
      nnzPerRow.back()++;
      colIdx.push_back(J);
      val.push_back(V);
    }

    // Done with localNZ map; free it.
    localNZmap_t().swap(localNZ);
  }

  // Compute prefix sum in offsets array
  offsets.resize(rowIdx.size() + 1);
  offsets[0] = 0;
  for (size_t row = 0; row < static_cast<size_t>(rowIdx.size()); row++)
    offsets[row+1] = offsets[row] + nnzPerRow[row];

  if (useTimers) {
    timer = Teuchos::null;
    timer = rcp(new Teuchos::TimeMonitor(
                   *Teuchos::TimeMonitor::getNewTimer("RSF insertNonzeros")));
  }

  // Create a new RowMap with only rows having non-zero entries.
  size_t dummy = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
  Teuchos::RCP<const Tpetra::Map<> > rowMap = 
       Teuchos::rcp(new Tpetra::Map<>(dummy, rowIdx(), 0, comm));

  Teuchos::RCP<sparse_matrix_type> A = 
           rcp(new sparse_matrix_type(rowMap, nnzPerRow()));

  //  Insert the global values into the matrix row-by-row.
  if (verbose && me == 0) 
    std::cout << "Inserting global values" << std::endl;

  if(readPerProcess){
    const scalar_type ONE = Teuchos::ScalarTraits<scalar_type>::one();
    for (int i = 0; i < rowIdx.size(); i++) {
      size_t nnz = nnzPerRow[i];
      size_t off = offsets[i];
      val.assign(nnz, ONE);
      // ReadPerProcess routine does not read any numeric values from the file,
      // So we insert dummy values here. 
      A->insertGlobalValues(rowIdx[i], colIdx(off, nnz), val());
    }
  }
  else {
    for (int i = 0; i < rowIdx.size(); i++) {
      size_t nnz = nnzPerRow[i];
      size_t off = offsets[i];
      A->insertGlobalValues(rowIdx[i], colIdx(off, nnz), val(off, nnz));
    }
  }

  // free memory that is no longer needed
  Teuchos::Array<size_t>().swap(nnzPerRow);
  Teuchos::Array<size_t>().swap(offsets);
  Teuchos::Array<global_ordinal_type>().swap(rowIdx);
  Teuchos::Array<global_ordinal_type>().swap(colIdx);
  Teuchos::Array<scalar_type>().swap(val);

  if (useTimers)
    timer = Teuchos::null;

  if (callFillComplete) {

    if (verbose && me == 0) 
      std::cout << "Building vector maps" << std::endl;

    if (useTimers) {
      timer = Teuchos::null;
      timer = rcp(new Teuchos::TimeMonitor(
                     *Teuchos::TimeMonitor::getNewTimer("RSF buildVectorMaps")));
    }

    // Build domain map that corresponds to distribution 
    Teuchos::Array<global_ordinal_type> vectorSet;
    for (global_ordinal_type i = 0;
                             i < static_cast<global_ordinal_type>(nCol); i++) 
      if (dist->VecMine(i)) vectorSet.push_back(i);

    dummy = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
    Teuchos::RCP<const Tpetra::Map<> > domainMap = 
           Teuchos::rcp(new Tpetra::Map<>(dummy, vectorSet(), 0, comm));
  
    Teuchos::Array<global_ordinal_type>().swap(vectorSet);

    // Build range map that corresponds to distribution
    for (global_ordinal_type i = 0; 
                             i < static_cast<global_ordinal_type>(nRow); i++) 
      if (dist->VecMine(i)) vectorSet.push_back(i);
  
    dummy = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
    Teuchos::RCP<const Tpetra::Map<> > rangeMap = 
           Teuchos::rcp(new Tpetra::Map<>(dummy, vectorSet(), 0, comm));

    Teuchos::Array<global_ordinal_type>().swap(vectorSet);

    // FillComplete the matrix
    if (useTimers) {
      timer = Teuchos::null;
      timer = rcp(new Teuchos::TimeMonitor(
                     *Teuchos::TimeMonitor::getNewTimer("RSF fillComplete")));
    }

    if (verbose && me == 0) 
      std::cout << "Calling fillComplete" << std::endl;

    A->fillComplete(domainMap, rangeMap);

    if (useTimers) 
      timer = Teuchos::null;

    if (verbose) {
      std::cout << "\nRank " << A->getComm()->getRank() 
                << "  nRows " << A->getLocalNumRows()
                << "  minRow " << A->getRowMap()->getMinGlobalIndex()
                << "  maxRow " << A->getRowMap()->getMaxGlobalIndex()
                << "\nRank " << A->getComm()->getRank() 
                << "  nCols " << A->getLocalNumCols()
                << "  minCol " << A->getColMap()->getMinGlobalIndex()
                << "  maxCol " << A->getColMap()->getMaxGlobalIndex()
                << "\nRank " << A->getComm()->getRank() 
                << "  nDomain " << A->getDomainMap()->getLocalNumElements()
                << "  minDomain " << A->getDomainMap()->getMinGlobalIndex()
                << "  maxDomain " << A->getDomainMap()->getMaxGlobalIndex()
                << "\nRank " << A->getComm()->getRank() 
                << "  nRange " << A->getRangeMap()->getLocalNumElements()
                << "  minRange " << A->getRangeMap()->getMinGlobalIndex()
                << "  maxRange " << A->getRangeMap()->getMaxGlobalIndex()
                << "\nRank " << A->getComm()->getRank() 
                << "  nEntries " << A->getLocalNumEntries()
                << std::endl;
    }
  }

  return A;
}



#endif