Julia vs C++

cartesian-grid.cpp

// compile with either:
// clang++ -std=c++11 -O2 cartesian-grid.cpp -o cartesian-grid
// g++ -std=c++11 -O2 cartesian-grid.cpp -o cartesian-grid
#include <chrono>
#include <iostream>
#include <vector>

typedef double scalar_t;
typedef int64_t idx_t;

struct CartesianGrid
{
  idx_t nx;
  idx_t ny;
  scalar_t h;

  idx_t nnodes() const { return nx*ny; }

  std::pair<idx_t,idx_t> xyindices(const idx_t i) const
  {
    const idx_t ix = i % nx;
    return std::make_pair(ix, (i-ix) / nx);
  }
};

template<typename VecT>
void fill_array(VecT& A, const CartesianGrid& g)
{
  auto start = std::chrono::steady_clock::now();
  const idx_t n_elems = A.size();
  for(idx_t i = 0; i != n_elems; ++i)
  {
    const auto p = g.xyindices(i);
    const idx_t ix = p.first;
    const idx_t iy = p.second;
    A[i] = ix + iy;
  }
  auto end = std::chrono::steady_clock::now();
  std::cout << std::chrono::duration<double>(end-start).count() << " s" << std::endl;
}

int main()
{
  CartesianGrid grid = {1001,1001,1.0};

  std::vector<idx_t> A(grid.nnodes());

  fill_array(A,grid);
  fill_array(A,grid);
  fill_array(A,grid);

  return 0;
}

cartesian-grid.jl

immutable CartesianGrid
  nx::Int # Number of points in the X direction
  ny::Int # Number of points in the Y direction
  h::Float64 # Spacing
end

nnodes(g::CartesianGrid) = g.nx*g.ny

function xyindices(g::CartesianGrid, i)
  ix = i % g.nx
  iy = (i-ix) ÷ g.nx
  return (ix,iy)
end

function fill_array!(A, g::CartesianGrid)
  for i in linearindices(A)
    (ix,iy) = xyindices(g,i)
    A[i] = ix + iy
  end
end

grid = CartesianGrid(1001,1001,1)
A = zeros(Int,nnodes(grid))

@time fill_array!(A,grid)
@time fill_array!(A,grid)
@time fill_array!(A,grid)