EpetraExt_CrsMatrixIn.cpp

Go to the documentation of this file.

//@HEADER
// ***********************************************************************
//
//     EpetraExt: Epetra Extended - Linear Algebra Services Package
//                 Copyright (2011) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
//@HEADER
#include "Epetra_ConfigDefs.h"
#include "EpetraExt_CrsMatrixIn.h"
#include "Epetra_Comm.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_IntVector.h"
#include "Epetra_IntSerialDenseVector.h"
#include "Epetra_Import.h"
#include "Epetra_Time.h"
#include "Epetra_Util.h"
 
#include <fstream>
 
#if defined(__PGI)
#include <sstream>
#endif
 
// 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 "chunk_read" can be changed to modify the size of the chunks
// read from the file.  Larger values of chunk_read lead to faster execution
// and greater memory use.
 
using namespace EpetraExt;
namespace EpetraExt {
 
template<typename int_type>
static void sort_three(int_type* list, int_type *parlista, double *parlistb,
                       int start, int end);
 
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
int MatrixMarketFileToCrsMatrix(const char *filename,
    const Epetra_Comm & comm, Epetra_CrsMatrix * & A)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix(const char *filename,
    const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A,
                                                   0, 0, 0, 0, transpose));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix(const char *filename,
    const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose,
    const bool verbose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename, comm, A,
                                                   0, 0, 0, 0,
                                                   transpose, verbose));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix(const char *filename,
    const Epetra_Map & rowMap, const Epetra_Map& rangeMap,
    const Epetra_Map& domainMap, Epetra_CrsMatrix * & A, const bool transpose,
    const bool verbose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename,
                                                   rowMap.Comm(), A,
                                                   &rowMap, 0,
                                                   &rangeMap, &domainMap,
                                                   transpose, verbose));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix(const char *filename,
    const Epetra_Map & rowMap, Epetra_CrsMatrix * & A, const bool transpose,
    const bool verbose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename,
                                                   rowMap.Comm(), A,
                                                   &rowMap, 0, 0, 0,
                                                   transpose, verbose));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix(const char *filename,
    const Epetra_Map & rowMap, const Epetra_Map & colMap,
    Epetra_CrsMatrix * & A, const bool transpose,
    const bool verbose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename,
                                                   rowMap.Comm(), A,
                                                   &rowMap, &colMap, 0, 0,
                                                   transpose, verbose));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix(const char *filename,
    const Epetra_Map & rowMap, const Epetra_Map & colMap,
    const Epetra_Map& rangeMap, const Epetra_Map& domainMap,
    Epetra_CrsMatrix * & A, const bool transpose,
    const bool verbose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle(filename,
                                                   rowMap.Comm(), A,
                                                   &rowMap, &colMap,
                                                   &rangeMap, &domainMap,
                                                   transpose, verbose));
  return(0);
}
#endif
 
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
int MatrixMarketFileToCrsMatrix64(const char *filename,
    const Epetra_Comm & comm, Epetra_CrsMatrix * & A)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, comm, A));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix64(const char *filename,
    const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, comm, A,
                                                   0, 0, 0, 0, transpose));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix64(const char *filename,
    const Epetra_Comm & comm, Epetra_CrsMatrix * & A, const bool transpose,
    const bool verbose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename, comm, A,
                                                   0, 0, 0, 0,
                                                   transpose, verbose));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix64(const char *filename,
    const Epetra_Map & rowMap, const Epetra_Map& rangeMap,
    const Epetra_Map& domainMap, Epetra_CrsMatrix * & A, const bool transpose,
    const bool verbose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename,
                                                   rowMap.Comm(), A,
                                                   &rowMap, 0,
                                                   &rangeMap, &domainMap,
                                                   transpose, verbose));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix64(const char *filename,
    const Epetra_Map & rowMap, Epetra_CrsMatrix * & A, const bool transpose,
    const bool verbose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename,
                                                   rowMap.Comm(), A,
                                                   &rowMap, 0, 0, 0,
                                                   transpose, verbose));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix64(const char *filename,
    const Epetra_Map & rowMap, const Epetra_Map & colMap,
    Epetra_CrsMatrix * & A, const bool transpose,
    const bool verbose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename,
                                                   rowMap.Comm(), A,
                                                   &rowMap, &colMap, 0, 0,
                                                   transpose, verbose));
  return(0);
}
 
int MatrixMarketFileToCrsMatrix64(const char *filename,
    const Epetra_Map & rowMap, const Epetra_Map & colMap,
    const Epetra_Map& rangeMap, const Epetra_Map& domainMap,
    Epetra_CrsMatrix * & A, const bool transpose,
    const bool verbose)
{
  EPETRA_CHK_ERR(MatrixMarketFileToCrsMatrixHandle64(filename,
                                                   rowMap.Comm(), A,
                                                   &rowMap, &colMap,
                                                   &rangeMap, &domainMap,
                                                   transpose, verbose));
  return(0);
}
#endif
 
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
int MatrixMarketFileToCrsMatrixHandle(const char *filename,
  const Epetra_Comm & comm,
  Epetra_CrsMatrix * & A,
  const Epetra_Map * rowMap,
  const Epetra_Map * colMap,
  const Epetra_Map * rangeMap,
  const Epetra_Map * domainMap,
  const bool transpose,
  const bool verbose
)
{
  const int chunk_read = 500000;  //  Modify this variable to change the
                                  //  size of the chunks read from the file.
  const int headerlineLength = 257;
  const int lineLength = 81;
  const int tokenLength = 35;
  char line[headerlineLength];
  char token1[tokenLength];
  char token2[tokenLength];
  char token3[tokenLength];
  char token4[tokenLength];
  char token5[tokenLength];
  int M, N, NZ;      // Matrix dimensions
  int me = comm.MyPID();
 
  Epetra_Time timer(comm);
 
  // Make sure domain and range maps are either both defined or undefined
  if ((domainMap!=0 && rangeMap==0) || (domainMap==0 && rangeMap!=0)) {
    EPETRA_CHK_ERR(-3);
  }
 
  // check maps to see if domain and range are 1-to-1
 
  if (domainMap!=0) {
    if (!domainMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
    if (!rangeMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
  }
  else {
    // If domain and range are not specified,
    // rowMap becomes both and must be 1-to-1 if specified
    if (rowMap!=0) {
      if (!rowMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
    }
  }
 
  FILE * handle = 0;
  int rv=0;
  if (me == 0) {
    if (verbose) std::cout << "Reading MatrixMarket file " << filename << std::endl;
    handle = fopen(filename,"r");  // Open file
    if (handle == 0) {
      rv = -1; // file not found
    }
 
    // Check first line, which should be
    // %%MatrixMarket matrix coordinate real general
    if( (rv != -1) && fgets(line, headerlineLength, handle)==0 ) {
      if (handle!=0) fclose(handle);
      rv = -1;
    }
    if( (rv != -1) && sscanf(line, "%s %s %s %s %s", token1,token2,token3,token4,token5)==0) {
      if (handle!=0) fclose(handle);
      rv = -1;
    }
    if ( (rv != -1) && (strcmp(token1, "%%MatrixMarket") ||
        strcmp(token2, "matrix") ||
        strcmp(token3, "coordinate") ||
        strcmp(token4, "real") ||
        strcmp(token5, "general")) ) {
      if (handle!=0) fclose(handle);
      rv = -1;
    }
 
    // Next, strip off header lines (which start with "%")
    if (rv != -1) {
      do {
        if(fgets(line, headerlineLength, handle)==0) {
          if (handle!=0) fclose(handle);
          rv = -1;
          break;
        }
      } while (line[0] == '%');
    }
 
    // Next get problem dimensions: M, N, NZ
    if((rv != -1) && sscanf(line, "%d %d %d", &M, &N, &NZ)==0) {
      if (handle!=0) fclose(handle);
      rv = -1;
    }
  }
  // If there's an error reading the file, broadcast to all processes in communicator.
  comm.Broadcast(&rv, 1, 0);
  if (rv==-1)
    EPETRA_CHK_ERR(-1);
  comm.Broadcast(&M, 1, 0);
  comm.Broadcast(&N, 1, 0);
  comm.Broadcast(&NZ, 1, 0);
 
  // Now create matrix using user maps if provided.
 
 
  // Now read in chunks of triplets and broadcast them to other processors.
  char *buffer = new char[chunk_read*lineLength];
  int nchunk;
  int nmillion = 0;
  int nread = 0;
  int rlen;
 
  // Storage for this processor's nonzeros.
  const int localblock = 100000;
  int localsize = NZ / comm.NumProc() + localblock;
  int *iv = (int *) malloc(localsize * sizeof(int));
  int *jv = (int *) malloc(localsize * sizeof(int));
  double *vv = (double *) malloc(localsize * sizeof(double));
  int lnz = 0;   //  Number of non-zeros on this processor.
 
  if (!iv || !jv || !vv)
    EPETRA_CHK_ERR(-1);
 
  Epetra_Map *rowMap1;
  bool allocatedHere=false;
  if (rowMap != 0)
    rowMap1 = const_cast<Epetra_Map *>(rowMap);
  else {
    rowMap1 = new Epetra_Map(M, 0, comm);
    allocatedHere = true;
  }
  int ioffset = rowMap1->IndexBase()-1;
  int joffset = (colMap != 0 ? colMap->IndexBase()-1 : ioffset);
 
  int rowmajor = 1;  // Assume non-zeros are listed in row-major order,
  int prevrow = -1;  // but test to detect otherwise.  If non-zeros
                     // are row-major, we can avoid the sort.
 
  while (nread < NZ) {
    if (NZ-nread > chunk_read) nchunk = chunk_read;
    else nchunk = NZ - nread;
 
    if (me == 0) {
      char *eof;
      rlen = 0;
      for (int i = 0; i < nchunk; i++) {
        eof = fgets(&buffer[rlen],lineLength,handle);
        if (eof == NULL) {
          fprintf(stderr, "%s", "Unexpected end of matrix file.");
          EPETRA_CHK_ERR(-1);
        }
        rlen += strlen(&buffer[rlen]);
      }
      buffer[rlen++]='\n';
    }
    comm.Broadcast(&rlen, 1, 0);
    comm.Broadcast(buffer, rlen, 0);
 
    buffer[rlen++] = '\0';
    nread += nchunk;
 
    // Process the received data, saving non-zeros belonging on this proc.
    char *lineptr = buffer;
    for (rlen = 0; rlen < nchunk; rlen++) {
      char *next = strchr(lineptr,'\n');
      int I = atoi(strtok(lineptr," \t\n")) + ioffset;
      int J = atoi(strtok(NULL," \t\n")) + joffset;
      double V = atof(strtok(NULL," \t\n"));
      lineptr = next + 1;
      if (transpose) {
        // swap I and J indices.
        int tmp = I;
        I = J;
        J = tmp;
      }
      if (rowMap1->MyGID(I)) {
        //  This processor keeps this non-zero.
        if (lnz >= localsize) {
          // Need more memory to store nonzeros.
          localsize += localblock;
          iv = (int *) realloc(iv, localsize * sizeof(int));
          jv = (int *) realloc(jv, localsize * sizeof(int));
          vv = (double *) realloc(vv, localsize * sizeof(double));
        }
        iv[lnz] = I;
        jv[lnz] = J;
        vv[lnz] = V;
        lnz++;
        if (I < prevrow) rowmajor = 0;
        prevrow = I;
      }
    }
 
    // Status check
    if (nread / 1000000 > nmillion) {
      nmillion++;
      if (verbose && me == 0) std::cout << nmillion << "M ";
    }
  }
 
  delete [] buffer;
 
  if (!rowmajor) {
    // Sort into row-major order (by iv) so can insert entire rows at once.
    // Reorder jv and vv to parallel iv.
    if (verbose && me == 0) std::cout << std::endl << "   Sorting local nonzeros" << std::endl;
    sort_three(iv, jv, vv, 0, lnz-1);
  }
 
  //  Compute number of entries per local row for use in constructor;
  //  saves reallocs in FillComplete.
 
  //  Now construct the matrix.
  //  Count number of entries in each row so can use StaticProfile=2.
  if (verbose && me == 0) std::cout << std::endl << "   Constructing the matrix" << std::endl;
  int numRows = rowMap1->NumMyElements();
  int *numNonzerosPerRow = new int[numRows];
  for (int i = 0; i < numRows; i++) numNonzerosPerRow[i] = 0;
  for (int i = 0; i < lnz; i++)
    numNonzerosPerRow[rowMap1->LID(iv[i])]++;
 
  if (rowMap!=0 && colMap !=0)
    A = new Epetra_CrsMatrix(Copy, *rowMap, *colMap, numNonzerosPerRow);
  else if (rowMap!=0) {
    // Construct with StaticProfile=true since we know numNonzerosPerRow.
    // Less memory will be needed in FillComplete.
    A = new Epetra_CrsMatrix(Copy, *rowMap, numNonzerosPerRow, true);
  }
  else {
    // Construct with StaticProfile=true since we know numNonzerosPerRow.
    // Less memory will be needed in FillComplete.
    A = new Epetra_CrsMatrix(Copy, *rowMap1, numNonzerosPerRow, true);
  }
  A->SetTracebackMode(1);
 
  // Rows are inserted in ascending global number, and the insertion uses numNonzerosPerRow.
  // Therefore numNonzerosPerRow must be permuted, as it was created in ascending local order.
  int *gidList = new int[numRows];
  for (int i = 0; i < numRows; i++) gidList[i] = rowMap1->GID(i);
  Epetra_Util::Sort(true,numRows,gidList,0,0,1,&numNonzerosPerRow);
  delete [] gidList;
 
  //  Insert the global values into the matrix row-by-row.
  if (verbose && me == 0) std::cout << "   Inserting global values" << std::endl;
  {
    int i = 0;
    for (int sum = 0; i < numRows; i++) {
      if (numNonzerosPerRow[i]) {
        int ierr = A->InsertGlobalValues(iv[sum], numNonzerosPerRow[i],
                                         &vv[sum], &jv[sum]);
        if (ierr<0) EPETRA_CHK_ERR(ierr);
        sum += numNonzerosPerRow[i];
      }
    }
  }
 
  delete [] numNonzerosPerRow;
  free(iv);
  free(jv);
  free(vv);
 
  if (verbose && me == 0) std::cout << "   Completing matrix fill" << std::endl;
  if (rangeMap != 0 && domainMap != 0) {
    EPETRA_CHK_ERR(A->FillComplete(*domainMap, *rangeMap));
  }
  else if (M!=N) {
    Epetra_Map newDomainMap(N,rowMap1->IndexBase(), comm);
    EPETRA_CHK_ERR(A->FillComplete(newDomainMap, *rowMap1));
  }
  else {
    EPETRA_CHK_ERR(A->FillComplete());
  }
 
  if (allocatedHere) delete rowMap1;
 
  if (handle!=0) fclose(handle);
  double dt = timer.ElapsedTime();
  if (verbose && me == 0) std::cout << "File Read time (secs):  " << dt << std::endl;
  return(0);
}
#endif
 
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
int MatrixMarketFileToCrsMatrixHandle64(const char *filename,
  const Epetra_Comm & comm,
  Epetra_CrsMatrix * & A,
  const Epetra_Map * rowMap,
  const Epetra_Map * colMap,
  const Epetra_Map * rangeMap,
  const Epetra_Map * domainMap,
  const bool transpose,
  const bool verbose
)
{
  const int chunk_read = 500000;  //  Modify this variable to change the
                                  //  size of the chunks read from the file.
  const int headerlineLength = 257;
  const int lineLength = 81;
  const int tokenLength = 35;
  char line[headerlineLength];
  char token1[tokenLength];
  char token2[tokenLength];
  char token3[tokenLength];
  char token4[tokenLength];
  char token5[tokenLength];
  long long M, N, NZ;      // Matrix dimensions
  int me = comm.MyPID();
 
  Epetra_Time timer(comm);
 
  // Make sure domain and range maps are either both defined or undefined
  if ((domainMap!=0 && rangeMap==0) || (domainMap==0 && rangeMap!=0)) {
    EPETRA_CHK_ERR(-3);
  }
 
  // check maps to see if domain and range are 1-to-1
 
  if (domainMap!=0) {
    if (!domainMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
    if (!rangeMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
  }
  else {
    // If domain and range are not specified,
    // rowMap becomes both and must be 1-to-1 if specified
    if (rowMap!=0) {
      if (!rowMap->UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
    }
  }
 
  FILE * handle = 0;
  if (me == 0) {
    if (verbose) std::cout << "Reading MatrixMarket file " << filename << std::endl;
    handle = fopen(filename,"r");  // Open file
    if (handle == 0)
      EPETRA_CHK_ERR(-1); // file not found
 
    // Check first line, which should be
    // %%MatrixMarket matrix coordinate real general
    if(fgets(line, headerlineLength, handle)==0) {
      if (handle!=0) fclose(handle);
      EPETRA_CHK_ERR(-1);
    }
    if(sscanf(line, "%s %s %s %s %s", token1,token2,token3,token4,token5)==0) {
      if (handle!=0) fclose(handle);
      EPETRA_CHK_ERR(-1);
    }
    if (strcmp(token1, "%%MatrixMarket") ||
        strcmp(token2, "matrix") ||
        strcmp(token3, "coordinate") ||
        strcmp(token4, "real") ||
        strcmp(token5, "general")) {
      if (handle!=0) fclose(handle);
      EPETRA_CHK_ERR(-1);
    }
 
    // Next, strip off header lines (which start with "%")
    do {
      if(fgets(line, headerlineLength, handle)==0) {
        if (handle!=0) fclose(handle);
        EPETRA_CHK_ERR(-1);
      }
    } while (line[0] == '%');
 
    // Next get problem dimensions: M, N, NZ
    if(sscanf(line, "%lld %lld %lld", &M, &N, &NZ)==0) {
      if (handle!=0) fclose(handle);
      EPETRA_CHK_ERR(-1);
    }
  }
  comm.Broadcast(&M, 1, 0);
  comm.Broadcast(&N, 1, 0);
  comm.Broadcast(&NZ, 1, 0);
 
  // Now create matrix using user maps if provided.
 
 
  // Now read in chunks of triplets and broadcast them to other processors.
  char *buffer = new char[chunk_read*lineLength];
  long long nchunk;
  int nmillion = 0;
  long long nread = 0;
  int rlen;
 
  // Storage for this processor's nonzeros.
  const int localblock = 100000;
  int localsize = (int) (NZ / comm.NumProc()) + localblock;
  long long *iv = (long long *) malloc(localsize * sizeof(long long));
  long long *jv = (long long *) malloc(localsize * sizeof(long long));
  double *vv = (double *) malloc(localsize * sizeof(double));
  int lnz = 0;   //  Number of non-zeros on this processor.
 
  if (!iv || !jv || !vv)
    EPETRA_CHK_ERR(-1);
 
  Epetra_Map *rowMap1;
  bool allocatedHere=false;
  if (rowMap != 0)
    rowMap1 = const_cast<Epetra_Map *>(rowMap);
  else {
    rowMap1 = new Epetra_Map(M, 0, comm);
    allocatedHere = true;
  }
  long long ioffset = rowMap1->IndexBase64()-1;
  long long joffset = (colMap != 0 ? colMap->IndexBase64()-1 : ioffset);
 
  int rowmajor = 1;  // Assume non-zeros are listed in row-major order,
  long long prevrow = -1;  // but test to detect otherwise.  If non-zeros
                     // are row-major, we can avoid the sort.
 
  while (nread < NZ) {
    if (NZ-nread > chunk_read) nchunk = chunk_read;
    else nchunk = NZ - nread;
 
    if (me == 0) {
      char *eof;
      rlen = 0;
      for (int i = 0; i < nchunk; i++) {
        eof = fgets(&buffer[rlen],lineLength,handle);
        if (eof == NULL) {
          fprintf(stderr, "%s", "Unexpected end of matrix file.");
          EPETRA_CHK_ERR(-1);
        }
        rlen += strlen(&buffer[rlen]);
      }
      buffer[rlen++]='\n';
    }
    comm.Broadcast(&rlen, 1, 0);
    comm.Broadcast(buffer, rlen, 0);
 
    buffer[rlen++] = '\0';
    nread += nchunk;
 
    // Process the received data, saving non-zeros belonging on this proc.
    char *lineptr = buffer;
    for (rlen = 0; rlen < nchunk; rlen++) {
      char *next = strchr(lineptr,'\n');
      char *endp;
      const int base = 10;
#if defined(_MSC_VER)
      long long I = _strtoi64(strtok(lineptr," \t\n"), &endp, base) + ioffset;
      long long J = _strtoi64(strtok(NULL," \t\n"), &endp, base) + joffset;
#else
#if defined(__PGI)
      long long I, J;
      std::istringstream ssI(strtok(lineptr," \t\n"));
      ssI >> I; I += ioffset;
      std::istringstream ssJ(strtok(NULL," \t\n"));
      ssJ >> J; J += joffset;
#else
      long long I = strtoll(strtok(lineptr," \t\n"), &endp, base) + ioffset;
      long long J = strtoll(strtok(NULL," \t\n"), &endp, base) + joffset;
#endif
#endif
      double V = atof(strtok(NULL," \t\n"));
      lineptr = next + 1;
      if (transpose) {
        // swap I and J indices.
        long long tmp = I;
        I = J;
        J = tmp;
      }
      if (rowMap1->MyGID(I)) {
        //  This processor keeps this non-zero.
        if (lnz >= localsize) {
          // Need more memory to store nonzeros.
          localsize += localblock;
          iv = (long long *) realloc(iv, localsize * sizeof(long long));
          jv = (long long *) realloc(jv, localsize * sizeof(long long));
          vv = (double *) realloc(vv, localsize * sizeof(double));
        }
        iv[lnz] = I;
        jv[lnz] = J;
        vv[lnz] = V;
        lnz++;
        if (I < prevrow) rowmajor = 0;
        prevrow = I;
      }
    }
 
    // Status check
    if (nread / 1000000 > nmillion) {
      nmillion++;
      if (verbose && me == 0) std::cout << nmillion << "M ";
    }
  }
 
  delete [] buffer;
 
  if (!rowmajor) {
    // Sort into row-major order (by iv) so can insert entire rows at once.
    // Reorder jv and vv to parallel iv.
    if (verbose && me == 0) std::cout << std::endl << "   Sorting local nonzeros" << std::endl;
    sort_three(iv, jv, vv, 0, lnz-1);
  }
 
  //  Compute number of entries per local row for use in constructor;
  //  saves reallocs in FillComplete.
 
  //  Now construct the matrix.
  //  Count number of entries in each row so can use StaticProfile=2.
  if (verbose && me == 0) std::cout << std::endl << "   Constructing the matrix" << std::endl;
  int numRows = rowMap1->NumMyElements();
  int *numNonzerosPerRow = new int[numRows];
  for (int i = 0; i < numRows; i++) numNonzerosPerRow[i] = 0;
  for (int i = 0; i < lnz; i++)
    numNonzerosPerRow[rowMap1->LID(iv[i])]++;
 
  if (rowMap!=0 && colMap !=0)
    A = new Epetra_CrsMatrix(Copy, *rowMap, *colMap, numNonzerosPerRow);
  else if (rowMap!=0) {
    // Construct with StaticProfile=true since we know numNonzerosPerRow.
    // Less memory will be needed in FillComplete.
    A = new Epetra_CrsMatrix(Copy, *rowMap, numNonzerosPerRow, true);
  }
  else {
    // Construct with StaticProfile=true since we know numNonzerosPerRow.
    // Less memory will be needed in FillComplete.
    A = new Epetra_CrsMatrix(Copy, *rowMap1, numNonzerosPerRow, true);
  }
  A->SetTracebackMode(1);
 
  // Rows are inserted in ascending global number, and the insertion uses numNonzerosPerRow.
  // Therefore numNonzerosPerRow must be permuted, as it was created in ascending local order.
  long long *gidList = new long long[numRows];
  for (int i = 0; i < numRows; i++) gidList[i] = rowMap1->GID64(i);
  Epetra_Util::Sort(true,numRows,gidList,0,0,1,&numNonzerosPerRow,0,0);
  delete [] gidList;
 
  //  Insert the global values into the matrix row-by-row.
  if (verbose && me == 0) std::cout << "   Inserting global values" << std::endl;
  {
    int i = 0;
    for (int sum = 0; i < numRows; i++) {
      if (numNonzerosPerRow[i]) {
        int ierr = A->InsertGlobalValues(iv[sum], numNonzerosPerRow[i],
                                         &vv[sum], &jv[sum]);
        if (ierr<0) EPETRA_CHK_ERR(ierr);
        sum += numNonzerosPerRow[i];
      }
    }
  }
 
  delete [] numNonzerosPerRow;
  free(iv);
  free(jv);
  free(vv);
 
  if (verbose && me == 0) std::cout << "   Completing matrix fill" << std::endl;
  if (rangeMap != 0 && domainMap != 0) {
    EPETRA_CHK_ERR(A->FillComplete(*domainMap, *rangeMap));
  }
  else if (M!=N) {
    Epetra_Map newDomainMap(N,rowMap1->IndexBase64(), comm);
    EPETRA_CHK_ERR(A->FillComplete(newDomainMap, *rowMap1));
  }
  else {
    EPETRA_CHK_ERR(A->FillComplete());
  }
 
  if (allocatedHere) delete rowMap1;
 
  if (handle!=0) fclose(handle);
  double dt = timer.ElapsedTime();
  if (verbose && me == 0) std::cout << "File Read time (secs):  " << dt << std::endl;
  return(0);
}
#endif
 
// Sorting values in array list in increasing order. Criteria is int.
// Also rearrange values in arrays parlista and partlistb
// to match the new order of list.
 
template<typename int_type>
static void quickpart_list_inc_int (
  int_type *list, int_type *parlista, double *parlistb,
  int start, int end, int *equal, int *larger)
{
int i;
int_type key, itmp;
double dtmp;
 
  key = list ? list[(end+start)/2] : 1;
 
  *equal = *larger = start;
  for (i = start; i <= end; i++)
    if (list[i] < key) {
      itmp                = parlista[i];
      parlista[i]         = parlista[*larger];
      parlista[(*larger)] = parlista[*equal];
      parlista[(*equal)]  = itmp;
      dtmp                = parlistb[i];
      parlistb[i]         = parlistb[*larger];
      parlistb[(*larger)] = parlistb[*equal];
      parlistb[(*equal)]  = dtmp;
      itmp                = list[i];
      list[i]             = list[*larger];
      list[(*larger)++]   = list[*equal];
      list[(*equal)++]    = itmp;
    }
    else if (list[i] == key) {
      itmp                = parlista[i];
      parlista[i]         = parlista[*larger];
      parlista[(*larger)] = itmp;
      dtmp                = parlistb[i];
      parlistb[i]         = parlistb[*larger];
      parlistb[(*larger)] = dtmp;
      list[i]             = list[*larger];
      list[(*larger)++]   = key;
    }
}
 
template<typename int_type>
static void sort_three(int_type* list, int_type *parlista, double *parlistb,
                       int start, int end)
{
int  equal, larger;
 
  if (start < end) {
    quickpart_list_inc_int(list, parlista, parlistb, start, end,
                           &equal, &larger);
    sort_three(list, parlista, parlistb, start,  equal-1);
    sort_three(list, parlista, parlistb, larger, end);
  }
}
 
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
int MatlabFileToCrsMatrix(const char *filename,
                                const Epetra_Comm & comm,
                                Epetra_CrsMatrix * & A)
{
  const int lineLength = 1025;
  char line[lineLength];
  int I, J;
  double V;
 
  FILE * handle = 0;
 
  handle = fopen(filename,"r");  // Open file
  if (handle == 0)
    EPETRA_CHK_ERR(-1); // file not found
 
  int numGlobalRows = 0;
  int numGlobalCols = 0;
  while(fgets(line, lineLength, handle)!=0) {
    if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
    if (I>numGlobalRows) numGlobalRows = I;
    if (J>numGlobalCols) numGlobalCols = J;
  }
 
  if (handle!=0) fclose(handle);
  Epetra_Map rangeMap(numGlobalRows, 0, comm);
  Epetra_Map domainMap(numGlobalCols, 0, comm);
  A = new Epetra_CrsMatrix(Copy, rangeMap, 0);
 
  // Now read in each triplet and store to the local portion of the matrix if the row is owned.
  const Epetra_Map & rowMap1 = A->RowMap();
 
  handle = 0;
 
  handle = fopen(filename,"r");  // Open file
  if (handle == 0)
    EPETRA_CHK_ERR(-1); // file not found
 
  while (fgets(line, lineLength, handle)!=0) {
    if(sscanf(line, "%d %d %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
    I--; J--; // Convert to Zero based
    if (rowMap1.MyGID(I)) {
      int ierr = A->InsertGlobalValues(I, 1, &V, &J);
      if (ierr<0) EPETRA_CHK_ERR(ierr);
    }
  }
 
  EPETRA_CHK_ERR(A->FillComplete(domainMap, rangeMap));
 
  if (handle!=0) fclose(handle);
  return(0);
}
#endif
 
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
int MatlabFileToCrsMatrix64(const char *filename,
                                const Epetra_Comm & comm,
                                Epetra_CrsMatrix * & A)
{
  const int lineLength = 1025;
  char line[lineLength];
  long long I, J;
  double V;
 
  FILE * handle = 0;
 
  handle = fopen(filename,"r");  // Open file
  if (handle == 0)
    EPETRA_CHK_ERR(-1); // file not found
 
  long long numGlobalRows = 0;
  long long numGlobalCols = 0;
  while(fgets(line, lineLength, handle)!=0) {
    if(sscanf(line, "%lld %lld %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
    if (I>numGlobalRows) numGlobalRows = I;
    if (J>numGlobalCols) numGlobalCols = J;
  }
 
  if (handle!=0) fclose(handle);
  Epetra_Map rangeMap(numGlobalRows, 0, comm);
  Epetra_Map domainMap(numGlobalCols, 0, comm);
  A = new Epetra_CrsMatrix(Copy, rangeMap, 0);
 
  // Now read in each triplet and store to the local portion of the matrix if the row is owned.
  const Epetra_Map & rowMap1 = A->RowMap();
 
  handle = 0;
 
  handle = fopen(filename,"r");  // Open file
  if (handle == 0)
    EPETRA_CHK_ERR(-1); // file not found
 
  while (fgets(line, lineLength, handle)!=0) {
    if(sscanf(line, "%lld %lld %lg\n", &I, &J, &V)==0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
    I--; J--; // Convert to Zero based
    if (rowMap1.MyGID(I)) {
      int ierr = A->InsertGlobalValues(I, 1, &V, &J);
      if (ierr<0) EPETRA_CHK_ERR(ierr);
    }
  }
 
  EPETRA_CHK_ERR(A->FillComplete(domainMap, rangeMap));
 
  if (handle!=0) fclose(handle);
  return(0);
}
#endif
 
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
int HypreFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix){
  int MyPID = comm.MyPID();
  // This double will be in the format we want for the extension besides the leading zero
  double filePID = (double)MyPID/(double)100000;
  std::ostringstream stream;
  // Using setprecision() puts it in the std::string
  stream << std::setiosflags(std::ios::fixed) << std::setprecision(5) << filePID;
  // Then just ignore the first character
  std::string fileName(filename);
  fileName += stream.str().substr(1,7);
  // Open the file
  std::ifstream file(fileName.c_str());
  std::string line;
  if(file.is_open()){
    std::getline(file, line);
    int ilower, iupper;
    std::istringstream istream(line);
    // The first line of the file has the beginning and ending rows
    istream >> ilower;
    istream >> iupper;
    // Using those we can create a row map
    Epetra_Map RowMap(-1, iupper-ilower+1, 0, comm);
    Matrix = new Epetra_CrsMatrix(Copy, RowMap, 0);
    int currRow = -1;
    int counter = 0;
    std::vector<int> indices;
    std::vector<double> values;
    while(!file.eof()){
      std::getline(file, line);
      std::istringstream lineStr(line);
      int row, col;
      double val;
      lineStr >> row;
      lineStr >> col;
      lineStr >> val;
      if(currRow == -1) currRow = row; // First line
      if(row == currRow){
        // add to the vector
        counter = counter + 1;
        indices.push_back(col);
        values.push_back(val);
      } else {
        Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
        indices.clear();
        values.clear();
        counter = 0;
        currRow = row;
        // make a new vector
        indices.push_back(col);
        values.push_back(val);
        counter = counter + 1;
      }
    }
    Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
    Matrix->Comm().Barrier();
    Matrix->FillComplete();
    file.close();
    return 0;
  } else {
    std::cout << "\nERROR:\nCouldn't open " << fileName << ".\n";
    return -1;
  }
}
#endif
 
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
int HypreFileToCrsMatrix64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix){
  int MyPID = comm.MyPID();
  // This double will be in the format we want for the extension besides the leading zero
  double filePID = (double)MyPID/(double)100000;
  std::ostringstream stream;
  // Using setprecision() puts it in the std::string
  stream << std::setiosflags(std::ios::fixed) << std::setprecision(5) << filePID;
  // Then just ignore the first character
  std::string fileName(filename);
  fileName += stream.str().substr(1,7);
  // Open the file
  std::ifstream file(fileName.c_str());
  std::string line;
  if(file.is_open()){
    std::getline(file, line);
    int ilower, iupper;
    std::istringstream istream(line);
    // The first line of the file has the beginning and ending rows
    istream >> ilower;
    istream >> iupper;
    // Using those we can create a row map
    Epetra_Map RowMap(-1LL, iupper-ilower+1, 0LL, comm);
    Matrix = new Epetra_CrsMatrix(Copy, RowMap, 0);
    long long currRow = -1;
    int counter = 0;
    std::vector<long long> indices;
    std::vector<double> values;
    while(!file.eof()){
      std::getline(file, line);
      std::istringstream lineStr(line);
      long long row, col;
      double val;
      lineStr >> row;
      lineStr >> col;
      lineStr >> val;
      if(currRow == -1) currRow = row; // First line
      if(row == currRow){
        // add to the vector
        counter = counter + 1;
        indices.push_back(col);
        values.push_back(val);
      } else {
        Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
        indices.clear();
        values.clear();
        counter = 0;
        currRow = row;
        // make a new vector
        indices.push_back(col);
        values.push_back(val);
        counter = counter + 1;
      }
    }
    Matrix->InsertGlobalValues(currRow, counter, &values[0], &indices[0]);
    Matrix->Comm().Barrier();
    Matrix->FillComplete();
    file.close();
    return 0;
  } else {
    std::cout << "\nERROR:\nCouldn't open " << fileName << ".\n";
    return -1;
  }
}
#endif
 
} // namespace EpetraExt