array.C

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

// In this case, N is a compile time global constant.
const unsigned int N = 1000;

// Note the weird syntax for typedef'ing a fixed size array...
typedef double Matrix[N][N];

// Function prototypes
void print(const Matrix & matrix);
void fill(Matrix & matrix);

// main program
int main(int argc, const char * argv[])
{
  // Fill A with random values in [0,1)
  Matrix A;
  fill(A);
  // print(A);

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

  // Matrix multiplication: C = A*B
  Matrix C;

  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<N; ++i)
    for (unsigned int j=0; j<N; ++j)
      for (unsigned int k=0; k<N; ++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.059728
// 550  0.092765
// 1000 --
// 2000 --

// This function relies on knowing the size N since we can't ask the
// array.
void print(const Matrix & matrix)
{
  for (unsigned int i=0; i<N; ++i)
    {
      for (unsigned int j=0; j<N; ++j)
        std::cout << matrix[i][j] << ' ';
      std::cout << std::endl;
    }
}

// This function relies on knowing the size N since we can't ask the
// array.
void fill(Matrix & matrix)
{
  for (unsigned int i=0; i<N; ++i)
    for (unsigned int j=0; j<N; ++j)
      matrix[i][j] = random()/static_cast<double>(RAND_MAX);
}

// Back of the envelope stack size estimation:
// echo "500*500*8*3/1024/1024" | bc -l
// 5.72204589843750000000
// echo "550*550*8*3/1024/1024" | bc -l
// 6.92367553710937500000
// echo "600*600*8*3/1024/1024" | bc -l
// 8.23974609375000000000
// 500 is OK
// 550 is OK
// 600 segfaults
// 750 segfaults
// 1000 segfaults