array_heap.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;

typedef double ** Matrix;

// Function prototypes
void print(const Matrix & matrix);
void resize(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;
  resize(A);
  fill(A);
  // print(A);

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

  // Matrix multiplication: C = A*B
  Matrix C;
  resize(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.276017
// 550  0.284725
// 1000 6.74756
// 2000 56.1554

// 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;
    }
}

void resize(Matrix & matrix)
{
  // Allocate N rows
  matrix = static_cast<double **>(malloc(N*sizeof(double *)));

  // Allocate N cols for each row
  for (unsigned int row=0; row<N; ++row)
    matrix[row] = static_cast<double *>(malloc(N*sizeof(double)));
}

// 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);
}