Array Benchmark

/*TO COMPILE: g++ -I /usr/include/eigen3 test.cpp -Wall -O3 -DNDEBUG -DBOOST_DISABLE_ASSERTS -fopenmp -ffast-math -funroll-loops -march=native -mtune=native -o array_test */

// Please note that Eigen version 3 is required


#include <Eigen/Dense>
#include <boost/multi_array.hpp>

#include <vector>
#include <iostream>

#include <omp.h>

using namespace std;

const size_t
  rows = 25000,
  cols = 25000,
  cols_subset = 10000,
  iterations = 1e1;

/*** SEQUENTIAL ACCESS ***/

template<typename MatrixType>
void sequential_access_test_eigen(const MatrixType &eigen_matrix)
{
    double
            tmp_sum = 0,
            tmp_time = omp_get_wtime();

    for(size_t i=0; i<iterations; i++)
    {
        for(size_t j=0; j<rows; j++)
        {
            for(size_t k=0; k<cols_subset; k++)
            {
                tmp_sum += eigen_matrix(j, k);
            }
        }
    }
    double sequential_access_time = omp_get_wtime() - tmp_time;

    std::cout << "SEQUENTIAL ACCESS (EIGEN)          :" << sequential_access_time << "s" << std::endl;
    std::cout << "tmp_sum == " << tmp_sum << std::endl;
}

template<typename ArrayType>
void sequential_access_test_boost(ArrayType &boost_matrix)
{
    double
            tmp_sum = 0,
            tmp_time = omp_get_wtime();

    for(size_t i=0; i<iterations; i++)
    {
        for(size_t j=0; j<rows; j++)
        {
            for(size_t k=0; k<cols_subset; k++)
            {
                tmp_sum += boost_matrix[j][k];
            }
        }
    }

    double sequential_access_time = omp_get_wtime() - tmp_time;
    std::cout << "SEQUENTIAL ACCESS (BOOST)          :" << sequential_access_time << "s" << std::endl;
    std::cout << "tmp_sum == " << tmp_sum << std::endl;
}

template<typename ArrayType>
void sequential_access_test_boost_mapped(ArrayType &boost_matrix)
{
    double
            tmp_sum = 0,
            tmp_time = omp_get_wtime();

    Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > boost_matrix_mapped(boost_matrix.data(), rows, cols);

    for(size_t i=0; i<iterations; i++)
    {
        for(size_t j=0; j<rows; j++)
        {
            for(size_t k=0; k<cols_subset; k++)
            {
                tmp_sum += boost_matrix_mapped(j,k);
            }
        }
    }

    double sequential_access_time = omp_get_wtime() - tmp_time;
    std::cout << "SEQUENTIAL ACCESS (BOOST_MAPPED)   :" << sequential_access_time << "s" << std::endl;
    std::cout << "tmp_sum == " << tmp_sum << std::endl;
}

/*** CONDITIONAL ACCESS ***/

template<typename MatrixType>
void conditional_access_test_eigen(const MatrixType &eigen_matrix)
{
    double
            tmp_sum = 0,
            tmp_time = omp_get_wtime();

    for(size_t i=0; i<iterations; i++)
    {
        for(size_t j=0; j<rows; j++)
        {
            for(size_t k=0; k<cols_subset; k++)
            {
                if(k%2==0)
                {
                    tmp_sum += eigen_matrix(j, k);
                }
            }
        }
    }
    double sequential_access_time = omp_get_wtime() - tmp_time;

    std::cout << "CONDITIONAL ACCESS (EIGEN)         :" << sequential_access_time << "s" << std::endl;
    std::cout << "tmp_sum == " << tmp_sum << std::endl;
}

template<typename ArrayType>
void conditional_access_test_boost(ArrayType &boost_matrix)
{
    double
            tmp_sum = 0,
            tmp_time = omp_get_wtime();

    for(size_t i=0; i<iterations; i++)
    {
        for(size_t j=0; j<rows; j++)
        {
            for(size_t k=0; k<cols_subset; k++)
            {
                if(k%2==0)
                {
                    tmp_sum += boost_matrix[j][k];
                }
            }
        }
    }

    double sequential_access_time = omp_get_wtime() - tmp_time;
    std::cout << "CONDITIONAL ACCESS (BOOST)         :" << sequential_access_time << "s" << std::endl;
    std::cout << "tmp_sum == " << tmp_sum << std::endl;
}

template<typename ArrayType>
void conditional_access_test_boost_mapped(ArrayType &boost_matrix)
{
    double
            tmp_sum = 0,
            tmp_time = omp_get_wtime();

    Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > boost_matrix_mapped(boost_matrix.data(), rows, cols);

    for(size_t i=0; i<iterations; i++)
    {
        for(size_t j=0; j<rows; j++)
        {
            for(size_t k=0; k<cols_subset; k++)
            {
                if(k%2==0)
                {
                    tmp_sum += boost_matrix_mapped(j,k);
                }
            }
        }
    }

    double sequential_access_time = omp_get_wtime() - tmp_time;
    std::cout << "CONDITIONAL ACCESS (BOOST_MAPPED)  :" << sequential_access_time << "s" << std::endl;
    std::cout << "tmp_sum == " << tmp_sum << std::endl;
}



int main()
{
    //// INITIALISE RANDOM TEST DATA ////

    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> eigen_matrix;
    eigen_matrix = Eigen::MatrixXd::Random(rows, cols);

    boost::multi_array<double, 2, Eigen::aligned_allocator<double> > boost_matrix(boost::extents[rows][cols]);

    for(size_t j=0; j<rows; j++)
    {
        for(size_t k=0; k<cols; k++)
        {
            boost_matrix[j][k] = eigen_matrix(j,k);
        }
    }

    //// RUN TESTS ////

    sequential_access_test_boost_mapped(boost_matrix);
    sequential_access_test_eigen(eigen_matrix);
    sequential_access_test_boost(boost_matrix);

    conditional_access_test_boost_mapped(boost_matrix);
    conditional_access_test_eigen(eigen_matrix);
    conditional_access_test_boost(boost_matrix);

    return 0;

}