matrix.C

#include <iostream>
#include <vector>
#include <chrono>
#include <stdlib.h> // RAND_MAX

typedef std::vector<std::vector<double>> Matrix;

// Function prototypes
void print(const Matrix & matrix);
void resize(Matrix & matrix, unsigned int N);
void fill(Matrix & matrix);

// main program
int main(int argc, const char * argv[])
{
  // Number of rows/columns in the matrix.
  const unsigned int N = 1000;

  // Fill A with random values in [0,1)
  Matrix A;
  resize(A, N);
  fill(A);
  // print(A);

  // Fill B with random values in [0,1)
  Matrix B;
  resize(B, N);
  fill(B);
  // print(B);

  // Matrix multiplication: C = A*B
  Matrix C;
  resize(C, N);

  typedef std::chrono::system_clock clock_type;
  std::chrono::time_point<clock_type> start, end;

  start = clock_type::now();

  // C_ij = A_ik * B_kj
  for (unsigned int i=0; i<A.size(); ++i)
    for (unsigned int j=0; j<B.size(); ++j)
      for (unsigned int k=0; k<C.size(); ++k)
        C[i][j] += A[i][k]*B[k][j];

  end = clock_type::now();

  clock_type::duration dur = end - start;

  // Print elapsed time in microseconds, but multiplied by 1.e-6,
  // since the duration_cast produces and integral number of
  // microseconds.
  std::cout << "Elapsed time: "
            << std::chrono::duration_cast<std::chrono::microseconds>(dur).count() * 1.e-6
            << " s\n";


  // print(C);
  // Print something to guarantee that the compiler did not optimize
  // away the entire calculation?
  std::cout << C[N/2][N/2] << std::endl;

  return 0;
}

// Results with -O3
// N    s
// 500  0.269248
// 550  0.282318
// 1000 6.848208
// 2000 60.8494

void print(const Matrix & matrix)
{
  for (unsigned int i=0; i<matrix.size(); ++i)
    {
      for (unsigned int j=0; j<matrix.size(); ++j)
        std::cout << matrix[i][j] << ' ';
      std::cout << std::endl;
    }
}

void resize(Matrix & matrix, unsigned int N)
{
  matrix.resize(N);

  for (unsigned int row=0; row<matrix.size(); ++row)
    matrix[row].resize(N);
}

void fill(Matrix & matrix)
{
  for (unsigned int i=0; i<matrix.size(); ++i)
    for (unsigned int j=0; j<matrix.size(); ++j)
      matrix[i][j] = random()/static_cast<double>(RAND_MAX);
}