(C++) Fastest way possible for reading in matrix files (arbitrary size)

Question

I'm developing a bioinformatic tool, which requires reading in millions of matrix files (average dimension = (20k, 20k)). They are tab-delimited text files, and they look something like:

0.53  0.11

0.24  0.33

Because the software reads the matrix files one at a time, memory is not an issue, but it's very slow. The following is my current function for reading in a matrix file. I first make a matrix object using a double pointer, then fill in the matrix by looping through an input file .

float** make_matrix(int nrow, int ncol, float val){
    float** M = new float *[nrow];
    for(int i = 0; i < nrow; i++) {
        M[i] = new float[ncol];
        for(int j = 0; j < ncol; j++) {
           M[i][j] = val;
        }
    }
    return M;
}


float** read_matrix(string fname, int dim_1, int dim_2){

    float** K = make_matrix(dim_1, dim_2, 0);

    ifstream ifile(fname);
    for (int i = 0; i < dim_1; ++i) {
        for (int j = 0; j < dim_2; ++j) {
            ifile >> K[i][j];
        }
    }

    ifile.clear();
    ifile.seekg(0, ios::beg);
    return K;
}

Is there a much faster way to do this? From my experience with python, reading in a matrix file using pandas is so much faster than using python for-loops. Is there a trick like that in c++?

---

(added)

Thanks so much everyone for all your suggestions and comments!

Answer 1

Just for fun, I measured the program posted above (using a 20,000x20,000 ASCII input file, as described) on my Mac Mini (3.2GHz i7 with SSD drive) and found that it took about 102 seconds to parse in the file using the posted code.

Then I wrote a version of the same function that uses the C stdio API (fopen()/fread()/fclose()) and does character-by-character parsing into a 1D float array. This implementation takes about 13 seconds to parse in the file on the same hardware, so it's about 7 times faster.

Both programs were compiled with g++ -O3 test_read_matrix.cpp.

float* faster_read_matrix(string fname, int numRows, int numCols)
{
    FILE * fpIn = fopen(fname.c_str(), "r");
    if (fpIn == NULL)
    {
       printf("Couldn't open file [%s] for input!\n", fname.c_str());
       return NULL;
    }

    float* K = new float[numRows*numCols];

    // We'll hold the current number in (numberBuf) until we're ready to parse it
    char numberBuf[128] = {'\0'};
    int numCharsInBuffer = 0;

    int curRow = 0, curCol = 0;
    while(curRow < numRows)
    {
       char tempBuf[4*1024];  // an arbitrary size
       const size_t bytesRead = fread(tempBuf, 1, sizeof(tempBuf), fpIn);
       if (bytesRead <= 0)
       {
          if (bytesRead < 0) perror("fread");
          break;
       }

       for (size_t i=0; i<bytesRead; i++)
       {
          const char c = tempBuf[i];
          if ((c=='.')||(c=='+')||(c=='-')||(isdigit(c)))
          {
             if ((numCharsInBuffer+1) < sizeof(numberBuf)) numberBuf[numCharsInBuffer++] = c;
             else
             {
                printf("Error, number string was too long for numberBuf!\n");
             }
          }
          else
          {
             if (numCharsInBuffer > 0)
             {
                // Parse the current number-chars we have assembled into (numberBuf) and reset (numberBuf) to empty
                numberBuf[numCharsInBuffer] = '\0';
                if (curCol < numCols) K[curRow*numCols+curCol] = strtod(numberBuf, NULL);
                else
                {
                    printf("Error, too many values in row %i!  (Expected %i, found at least %i)\n", curRow, numCols, curCol);
                }
                curCol++;
             }
             numCharsInBuffer = 0;

             if (c == '\n') 
             {
                curRow++;
                curCol = 0;
                if (curRow >= numRows) break;
             }
          }
       }
    }
    fclose(fpIn);

    if (curRow != numRows) printf("Warning:  I read %i lines in the file, but I expected there would be %i!\n", curRow, numRows);

    return K;
}

Answer 2

The fastest way, by far, is to change the way you write those files: write in binary format, two int first (width, height) then just dump your values.

You will be able to load it in just three read calls.

Answer 3

I am dissatisfied with Jeremy Friesner’s otherwise excellent answer because it:

• blames the problem to be with C++'s I/O system (which it is not)
• fixes the problem by circumventing the actual I/O problem without being explicit about how it is a significant contributor to speed
• modifies memory accesses which (may or may not) contribute to speed, and does so in a way that very large matrices may not be supported

The reason his code runs so much faster is because he removes the single most important bottleneck: unoptimized disk access. JWO’s original code can be brought to match with three extra lines of code:

float** read_matrix(std::string fname, int dim_1, int dim_2){

    float** K = make_matrix(dim_1, dim_2, 0);
    
    std::size_t buffer_size = 4*1024;               // 1
    char buffer[buffer_size];                       // 2

    std::ifstream ifile(fname);
    ifile.rdbuf()->pubsetbuf(buffer, buffer_size);  // 3
    
    for (int i = 0; i < dim_1; ++i) {
        for (int j = 0; j < dim_2; ++j) {
            ss >> K[i][j];
        }
    }

//    ifile.clear();
//    ifile.seekg(0, std::ios::beg);
    return K;
}

The addition exactly replicates Friesner’s design, but using the C++ library capabilities without all the extra programming grief on our end.

You’ll notice I also removed a couple lines at the bottom that should be inconsequential to program function and correctness, but which may cause a minor cumulative time issue as well. (If they are not inconsequential, that is a bug and should be fixed!)

How much difference this all makes depends entirely on the quality of the C++ Standard Library implementation. AFAIK the big three modern C++ compilers (MSVC, GCC, and Clang) all have sufficiently-optimized I/O handling to make the issue moot.

locale

One other thing that may also make a difference is to .imbue() the stream with the default "C" locale, which avoids a lot of special handling for numbers in locale-dependent formats other than what your files use. You only need to bother to do this if you have changed your global locale, though.

ifile.imbue(std::locale(""));

redundant initialization

Another thing that is killing your time is the effort to zero-initialize the array when you create it. Don’t do that if you don’t need it! (You don’t need it here because you know the total extents and will fill them properly. C++17 and later is nice enough to give you a zero value if the input stream goes bad, too. So you get zeros for unread values either way.)

dynamic memory block size

Finally, keeping memory accesses to an array of array should not significantly affect speed, but it still might be worth testing if you can change it. This is assuming that the resulting matrix will never be too large for the memory manager to return as a single block (and consequently crash your program).

A common design is to allocate the entire array as a single block, with the requested size plus size for the array of pointers to the rest of the block. This allows you to delete the array in a single delete[] statement. Again, I don’t believe this should be an optimization issue you need to care about until your profiler says so.

Answer 4

At the risk of the answer being considered incomplete (no code examples), I would like to add to the other answers additional options how to tackle the problem:

1. Use a binary format (width,height, values...) as file format and then use file mapping (MapViewOfFile() on Windows, mmap() or so on posix/unix systems).
   Then, you can simply point your "matrix structure" pointer to the mapped address space and you are done. And in case, you do something like sparse access to the matrix, it can even save some real IO. If you always do full access to all elements of the matrix (no sparse matrices etc.), it is still quite elegant and probably faster than malloc/read.

2. Replacements for c++ iostream, which is known to be quite slow and should not be used for performance critical stuff:
   Have a look at the {fmt} library, which has become quite popular in recent years and claims to be quite fast.

Back in the days, when I did a lot of numerics on large data sets, I always opted for binary files for storage. (It was back in the days, when the fastest CPU you get your hands on were the Pentium 1 (with the floating point bug :)). Back then, all was slower, memory was much more limited (we had MB not GB as units for RAM in our systems) and all in all, nearly 20 years have passed since.

So, as a refresher, I did write some code to show, how much faster than iostream and text files you can do if you do not have extra constraints (such as endianess of different cpus etc.).

So far, my little test only has an iostream and a binary file version with a) stdio fread() kind of loading and b) mmap(). Since I sit in front of a debian bullseye computer, my code uses linux specific stuff for the mmap() approach. To run it on Windows, you have to change a few lines of code and some includes.

Edit: I added a save function using {fmt} now as well.
Edit: I added a load function with stdio now as well.
Edit: To reduce memory workload, I reordered the code somewhat and now only keep 2 matrix instances in memory at any given time.

The program does the following:

1. create a 20k x 20k matrix in ram (in a struct named Matrix_t). With random values, slowly generated by std::random.
2. Write the matrix with iostream to a text file.
3. Write the matrix with stdio to a binary file.
4. Create a new matrix textMatrix by loading its data from the text file.
5. Create a new matrix inMemoryMatrix by loading its data from the binary file with a few fread() calls.
6. mmap() the binary file and use it under the name mappedMatrix.
7. Compare each of the loaded matrices to the original randomMatrix to see if the round-trip worked.

Here the results I got on my machine after compiling this work of wonder with clang++ -O3 -o fmatio fast-matrix-io.cpp -lfmt:

./fmatio
creating random matrix (20k x 20k) (27.0775seconds)
the first 10 floating values in randomMatrix are:
57970.2 -365700 -986079 44657.8 826968 -506928 668277 398241 -828176 394645
saveMatrixAsText_IOSTREAM()
saving matrix with iostream. (192.749seconds)
saveMatrixAsText_FMT(mat0_fmt.txt)
saving matrix with {fmt}. (34.4932seconds)
saveMatrixAsBinary()
saving matrix into a binary file. (30.7591seconds)
loadMatrixFromText_IOSTREAM()
loading matrix from text file with iostream. (102.074seconds)
randomMatrix == textMatrix
comparing randomMatrix with textMatrix. (0.125328seconds)
loadMatrixFromText_STDIO(mat0_fmt.txt)
loading matrix from text file with stdio. (71.2746seconds)
randomMatrix == textMatrix
comparing randomMatrix with textMatrix (stdio). (0.124684seconds)
loadMatrixFromBinary(mat0.bin)
loading matrix from binary file into memory. (0.495685seconds)
randomMatrix == inMemoryMatrix
comparing randomMatrix with inMemoryMatrix. (0.124206seconds)
mapMatrixFromBinaryFile(mat0.bin)
mapping a view to a matrix in a binary file. (4.5883e-05seconds)
randomMatrix == mappedMatrix
comparing randomMatrix with mappedMatrix. (0.158459seconds)

And here is the code:

#include <cinttypes>
#include <memory>
#include <random>
#include <iostream>
#include <fstream>
#include <cstring>
#include <string>
#include <chrono>
#include <limits>
#include <iomanip>

// includes for mmap()...
#include <sys/mman.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <cstdio>
#include <cstdlib>
#include <unistd.h>

// includes for {fmt}...
#include <fmt/core.h>
#include <fmt/os.h>

struct StopWatch {
  using Clock = std::chrono::high_resolution_clock;
  using TimePoint =
    std::chrono::time_point<Clock>;
  using Duration =
    std::chrono::duration<double>;
  
  void start(const char* description) {
    this->description = std::string(description);
    tstart = Clock::now();
  }
  void stop() {
    TimePoint tend = Clock::now();
    Duration elapsed = tend - tstart;
    std::cout << description << " (" << elapsed.count()
          << "seconds)" << std::endl;
  }
  TimePoint tstart;
  std::string description;
};

struct Matrix_t {
  uint32_t ncol;
  uint32_t nrow;
  float values[];
  inline uint32_t to_index(uint32_t col, uint32_t row) const {
    return ncol * row + col;
  }
};

template <class Initializer>
Matrix_t *createMatrix
( uint32_t ncol,
  uint32_t nrow,
  Initializer initFn
) {
  size_t nfloats = ncol*nrow;
  size_t nbytes = UINTMAX_C(8) + nfloats * sizeof(float);
  Matrix_t * result =
    reinterpret_cast<Matrix_t*>(operator new(nbytes));
  if (nullptr != result) {
    result->ncol = ncol;
    result->nrow = nrow;
    for (uint32_t row = 0; row < nrow; row++) {
      for (uint32_t col = 0; col < ncol; col++) {
    result->values[result->to_index(col,row)] =
      initFn(ncol,nrow,col,row);
      }
    }
  }
  return result;
}

void saveMatrixAsText_IOSTREAM(const char* filePath,
                   const Matrix_t* matrix) {
  std::cout << "saveMatrixAsText_IOSTREAM()" << std::endl;
  if (nullptr == matrix) {
    std::cout << "cannot save matrix - no matrix!" << std::endl;
  }
  std::ofstream outFile(filePath);
  if (outFile) {
    outFile << matrix->ncol << " " << matrix->nrow << std::endl;
    const auto defaultPrecision = outFile.precision();
    outFile.precision
      (std::numeric_limits<float>::max_digits10); 
    for (uint32_t row = 0; row < matrix->nrow; row++) {
      for (uint32_t col = 0; col < matrix->ncol; col++) {
    outFile << matrix->values[matrix->to_index(col,row)]
        << " ";
      }
      outFile << std::endl;
    }
  } else {
    std::cout << "could not open " << filePath << " for writing."
          << std::endl;
  }
}

void saveMatrixAsText_FMT(const char* filePath,
              const Matrix_t* matrix) {
  std::cout << "saveMatrixAsText_FMT(" << filePath << ")"
        << std::endl;
  if (nullptr == matrix) {
    std::cout << "cannot save matrix - no matrix!" << std::endl;
  }
  auto outFile = fmt::output_file(filePath);
  outFile.print("{} {}\n", matrix->ncol, matrix->nrow);
  for (uint32_t row = 0; row < matrix->nrow; row++) {
    outFile.print("{}", matrix->values[matrix->to_index(0,row)]);
    for (uint32_t col = 1; col < matrix->ncol; col++) {
      outFile.print(" {}",
            matrix->values[matrix->to_index(col,row)]);
    }
    outFile.print("\n");
  }
  
}

void saveMatrixAsBinary(const char* filePath,
            const Matrix_t* matrix) {
  std::cout << "saveMatrixAsBinary()" << std::endl;
  FILE * outFile = fopen(filePath, "wb");
  if (nullptr != outFile) {
    fwrite( &matrix->ncol, 4, 1, outFile);
    fwrite( &matrix->nrow, 4, 1, outFile);
    size_t nfloats = matrix->ncol * matrix->nrow;
    fwrite( &matrix->values, sizeof(float), nfloats, outFile);
    fclose(outFile);
  } else {
    std::cout << "could not open " << filePath << " for writing."
          << std::endl;
  }
}

Matrix_t* loadMatrixFromText_IOSTREAM(const char* filePath) {
  std::cout << "loadMatrixFromText_IOSTREAM()" << std::endl;
  std::ifstream inFile(filePath);
  if (inFile) {
    uint32_t ncol;
    uint32_t nrow;
    inFile >> ncol;
    inFile >> nrow;
    uint32_t nfloats = ncol * nrow;
    auto loader =
      [&inFile]
      (uint32_t , uint32_t , uint32_t , uint32_t )
      -> float
      {
    float value;
    inFile >> value;
    return value;
      };

    Matrix_t * matrix = createMatrix( ncol, nrow, loader);
    return matrix;
  } else {
    std::cout << "could not open " << filePath << "for reading."
          << std::endl;
  }
  return nullptr;
}

Matrix_t* loadMatrixFromText_STDIO(const char* filePath) {
  std::cout << "loadMatrixFromText_STDIO(" << filePath << ")"
        << std::endl;
  Matrix_t* matrix = nullptr;
  FILE * inFile = fopen(filePath, "rt");
  if (nullptr != inFile) {
    uint32_t ncol;
    uint32_t nrow;
    fscanf(inFile, "%d %d", &ncol, &nrow);
    auto loader =
      [&inFile]
      (uint32_t , uint32_t , uint32_t , uint32_t )
      -> float
      {
    float value;
    fscanf(inFile, "%f", &value);
    return value;
      };
    matrix = createMatrix( ncol, nrow, loader);
    fclose(inFile);
  } else {
    std::cout << "could not open " << filePath << "for reading."
          << std::endl;
  }
  return matrix;
}

Matrix_t* loadMatrixFromBinary(const char* filePath) {
  std::cout << "loadMatrixFromBinary(" << filePath << ")"
        << std::endl;
  FILE * inFile = fopen(filePath, "rb");
  if (nullptr != inFile) {
    uint32_t ncol;
    uint32_t nrow;
    fread( &ncol, 4, 1, inFile);
    fread( &nrow, 4, 1, inFile);
    uint32_t nfloats = ncol * nrow;
    uint32_t nbytes = nfloats * sizeof(float) + UINT32_C(8);
    Matrix_t* matrix =
      reinterpret_cast<Matrix_t*>
      (operator new (nbytes));
    if (nullptr != matrix) {
      matrix->ncol = ncol;
      matrix->nrow = nrow;
      fread( &matrix->values[0], sizeof(float), nfloats, inFile);
      return matrix;
    } else {
      std::cout << "could not find memory for the matrix."
        << std::endl;
    }
    fclose(inFile);
  } else {
    std::cout << "could not open file "
          << filePath << " for reading." << std::endl;
  }
  return nullptr;
}

void freeMatrix(Matrix_t* matrix) {
  operator delete(matrix);
}

Matrix_t* mapMatrixFromBinaryFile(const char* filePath) {
  std::cout << "mapMatrixFromBinaryFile(" << filePath << ")"
        << std::endl;
  Matrix_t * matrix = nullptr;
  int fd = open( filePath, O_RDONLY);
  if (-1 != fd) {
    struct stat sb;
    if (-1 != fstat(fd, &sb)) {
      auto fileSize = sb.st_size;
      matrix =
    reinterpret_cast<Matrix_t*>
    (mmap(nullptr, fileSize, PROT_READ, MAP_PRIVATE, fd, 0));
      if (nullptr == matrix) {
    std::cout << "mmap() failed!" << std::endl;
      }
    } else {
      std::cout << "fstat() failed!" << std::endl;
    }
    close(fd);
  } else {
    std::cout << "open() failed!" << std::endl;
  }
  return matrix;
}

void unmapMatrix(Matrix_t* matrix) {
  if (nullptr == matrix)
    return;
  size_t nbytes =
    UINTMAX_C(8) +
    sizeof(float) * matrix->ncol * matrix->nrow;
  munmap(matrix, nbytes);
}

bool areMatricesEqual( const Matrix_t* m1, const Matrix_t* m2) {
  if (nullptr == m1) return false;
  if (nullptr == m2) return false;
  if (m1->ncol != m2->ncol) return false;
  if (m1->nrow != m2->nrow) return false;
  // both exist and have same size...
  size_t nfloats = m1->ncol * m1->nrow;
  size_t nbytes = nfloats * sizeof(float);
  return 0 == memcmp( m1->values, m2->values, nbytes);
}

int main(int argc, const char* argv[]) {
    std::random_device rdev;
    std::default_random_engine reng(rdev());
    std::uniform_real_distribution<> rdist(-1.0E6F, 1.0E6F);
    StopWatch sw;
    
    auto randomInitFunction =
      [&reng,&rdist]
      (uint32_t ncol, uint32_t nrow, uint32_t col, uint32_t row)
      -> float
      {
    return rdist(reng);
      };
    sw.start("creating random matrix (20k x 20k)");
    Matrix_t * randomMatrix =
      createMatrix(UINT32_C(20000),
           UINT32_C(20000),
           randomInitFunction);
    sw.stop();
    if (nullptr != randomMatrix) {
      std::cout
    << "the first 10 floating values in randomMatrix are: "
    << std::endl;
      std::cout << randomMatrix->values[0];
      for (size_t i = 1; i < 10; i++) {
    std::cout << " " << randomMatrix->values[i];
      }
      std::cout << std::endl;
    
      sw.start("saving matrix with iostream.");
      saveMatrixAsText_IOSTREAM("mat0_iostream.txt", randomMatrix);
      sw.stop();
      sw.start("saving matrix with {fmt}.");
      saveMatrixAsText_FMT("mat0_fmt.txt", randomMatrix);
      sw.stop();
      sw.start("saving matrix into a binary file.");
      saveMatrixAsBinary("mat0.bin", randomMatrix);
      sw.stop();
      
      sw.start("loading matrix from text file with iostream.");
      Matrix_t* textMatrix =
    loadMatrixFromText_IOSTREAM("mat0_iostream.txt");
      sw.stop();
      sw.start("comparing randomMatrix with textMatrix.");
      if (!areMatricesEqual(randomMatrix, textMatrix)) {
    std::cout << "randomMatrix != textMatrix!" << std::endl;
      } else {
    std::cout << "randomMatrix == textMatrix" << std::endl;
      }
      sw.stop();
      freeMatrix(textMatrix);
      textMatrix = nullptr;
      
      sw.start("loading matrix from text file with stdio.");
      textMatrix =
    loadMatrixFromText_STDIO("mat0_fmt.txt");
      sw.stop();
      sw.start("comparing randomMatrix with textMatrix (stdio).");
      if (!areMatricesEqual(randomMatrix, textMatrix)) {
    std::cout << "randomMatrix != textMatrix!" << std::endl;
      } else {
    std::cout << "randomMatrix == textMatrix" << std::endl;
      }
      sw.stop();
      freeMatrix(textMatrix);
      textMatrix = nullptr;
      
      sw.start("loading matrix from binary file into memory.");
      Matrix_t* inMemoryMatrix =
    loadMatrixFromBinary("mat0.bin");
      sw.stop();
      sw.start("comparing randomMatrix with inMemoryMatrix.");
      if (!areMatricesEqual(randomMatrix, inMemoryMatrix)) {
    std::cout << "randomMatrix != inMemoryMatrix!"
          << std::endl;
      } else {
    std::cout << "randomMatrix == inMemoryMatrix" << std::endl;
      }
      sw.stop();
      freeMatrix(inMemoryMatrix);
      inMemoryMatrix = nullptr;
      
      sw.start("mapping a view to a matrix in a binary file.");
      Matrix_t* mappedMatrix =
    mapMatrixFromBinaryFile("mat0.bin");
      sw.stop();
      sw.start("comparing randomMatrix with mappedMatrix.");
      if (!areMatricesEqual(randomMatrix, mappedMatrix)) {
    std::cout << "randomMatrix != mappedMatrix!"
          << std::endl;
      } else {
    std::cout << "randomMatrix == mappedMatrix" << std::endl;
      }
      sw.stop();
      unmapMatrix(mappedMatrix);
      mappedMatrix = nullptr;

      freeMatrix(randomMatrix);
    } else {
      std::cout << "could not create random matrix!" << std::endl;
    }

    return 0;
}

Please note, that binary formats where you simply cast to a struct pointer also depend on how the compiler does alignment and padding within structures. In my case, I was lucky and it worked. On other systems, you might have to tweak a little (#pragma pack(4) or something along that line) to make it work.