main.cpp

#include <iostream>
#include <vector>
#include <immintrin.h> // For SIMD intrinsics (optional for AVX/AVX2)

#include <string>
#include <random>
#include <assert.h>
#include <algorithm>
#include <sys/stat.h>
#include <unistd.h>
#include <sys/time.h>
#include <Eigen/Dense>
// #include <sys/resource.h>

void display_timing(double wtime, std::string message="")
{
    int wh;      //, ch;
    int wmin;    //, cpmin;
    double wsec; //, csec;
    wh = (int)wtime / 3600;
    wmin = ((int)wtime % 3600) / 60;
    wsec = wtime - (3600. * wh + 60. * wmin);
    std::cout << message << " ";
    printf("WTime : %d hours and %d minutes and %.4f seconds.\n", wh, wmin, wsec);
}

double get_wall_time()
{
    /*!
    measure real passed time
    \return wall time in second
    */

    struct timeval time;
    if (gettimeofday(&time, NULL))
    {
        //  Handle error
        return 0;
    }
    return (double)time.tv_sec + (double)time.tv_usec * .000001;
}

void compareVec(const std::vector<double> &vec1, const std::vector<double> &vec2)
{
    if (vec1.size() != vec2.size())
    {
        std::cerr << "Vectors have different sizes!" << std::endl;
        return;
    }

    // compare the sum of differences
    double sum_diff = 0.0;
    for (size_t i = 0; i < vec1.size(); ++i)
    {
        sum_diff += std::abs(vec1[i] - vec2[i]);
    }
    std::cout << "Sum of differences: " << sum_diff << std::endl;
}



void matVec(const std::vector<std::vector<double>> &matrix,
            const std::vector<double> &vec,
            std::vector<double> &result)
{
    int n = matrix.size();

    for (int i = 0; i < n; ++i)
    {
        double sum = 0.0;

        for (int j = 0; j < n; ++j)
        {
            sum += matrix[i][j] * vec[j];
        }

        result[i] = sum;
    }
}

void matVecOpt(const std::vector<std::vector<double>> &matrix,
               const std::vector<double> &vec,
               std::vector<double> &result)
{
    int n = matrix.size();

    for (int i = 0; i < n; ++i)
    {
        double sum = 0.0;

        // Use pointer arithmetic for cache-friendly access
        const double *row = matrix[i].data();
        const double *vectorData = vec.data();

        // Process data in chunks (loop unrolling)
        int j = 0;
        for (; j + 4 <= n; j += 4)
        {
            sum += row[j] * vectorData[j] + row[j + 1] * vectorData[j + 1] + row[j + 2] * vectorData[j + 2] + row[j + 3] * vectorData[j + 3];
        }

        // Handle remaining elements
        for (; j < n; ++j)
        {
            sum += row[j] * vectorData[j];
        }

        result[i] = sum;
    }
}




int main()
{
    int n = 100; // Example size (you can increase this for benchmarking)

    // Initialize matrix and vector
    std::vector<std::vector<double>> matrix(n, std::vector<double>(n, 1.0));
    std::vector<double> vec(n, 2.0);
    std::vector<double> result(n, 0.0);

    // Perform optimized matrix-vector multiplication
    double wtime = get_wall_time();
    for (int i = 0; i < 100; ++i)
        matVec(matrix, vec, result);
    wtime = get_wall_time() - wtime;
    display_timing(wtime, "Regular C++  : ");

    // Perform regular matrix-vector multiplication
    std::fill(result.begin(), result.end(), 0.0);
    wtime = get_wall_time();
    for (int i = 0; i < 100; ++i)
        matVecOpt(matrix, vec, result);
    wtime = get_wall_time() - wtime;
    display_timing(wtime, "Optimized C++: ");

    
    Eigen::VectorXd vecEigen = Eigen::VectorXd::Constant(n, 2.0);
    Eigen::MatrixXd matrixEigen = Eigen::MatrixXd::Constant(n, n, 1.0);
    Eigen::VectorXd resultEigen = Eigen::VectorXd::Zero(n);

    // Perform Eigen matrix-vector multiplication
    wtime = get_wall_time();
    for (int i = 0; i < 100; ++i)
        resultEigen = matrixEigen * vecEigen;
    wtime = get_wall_time() - wtime;
    display_timing(wtime, "Eigen C++    : ");

    // Convert Eigen result to std::vector
    std::vector<double> resultEigenVec(resultEigen.data(), resultEigen.data() + resultEigen.size());
    compareVec(result, resultEigenVec);
    compareVec(result, result);


    return 0;
}


// g++ -O3 main.cpp -o prog && ./prog