matmul.cpp

/*
  Code to perform L1, L2 cache locality experiments
  https://kaushikghose.wordpress.com/2020/01/27/profiling-and-improving-cache-performance/
*/
#include <iostream>
#include <vector>

const size_t N = 1200;

class Matrix {

public:
    Matrix(size_t n, size_t m)
        : n(n)
        , m(m)
    {
        data.resize(n * m);
    }

    size_t rows() const { return n; }
    size_t cols() const { return m; }

    double&      get(size_t i, size_t j) { return data[_index(i, j)]; };
    const double get(size_t i, size_t j) const { return data[_index(i, j)]; };

    void fill_with_test_data()
    {
        double x = 1.0;
        for (size_t i = 0; i < n; i++) {
            for (size_t j = 0; j < m; j++) {
                get(i, j) = x++;
            }
        }
    }

    double sum()
    {
        double s = 0;
        for (auto d : data) {
            s += d;
        }
        return s;
    }

protected:
    virtual size_t _index(size_t i, size_t j) const = 0;

    size_t n, m;

    std::vector<double> data;
};

class MatrixRowMajor : public Matrix {

public:
    MatrixRowMajor(size_t n, size_t m)
        : Matrix(n, m)
    {
    }

private:
    size_t _index(size_t i, size_t j) const { return i * m + j; };
};

class MatrixColMajor : public Matrix {

public:
    MatrixColMajor(size_t n, size_t m)
        : Matrix(n, m)
    {
    }

private:
    size_t _index(size_t i, size_t j) const { return i + j * n; };
};

void mul(const Matrix& lhs, const Matrix& rhs, Matrix& ans)
{
    for (size_t i = 0; i < lhs.rows(); i++) {
        for (size_t j = 0; j < rhs.cols(); j++) {
            double v = 0;
            for (size_t k = 0; k < lhs.cols(); k++) {
                v += lhs.get(i, k) * rhs.get(k, j);
            }
            ans.get(i, j) = v;
        }
    }
}

MatrixRowMajor operator*(const Matrix& lhs, const Matrix& rhs)
{
    MatrixRowMajor ans(lhs.rows(), rhs.cols());
    mul(lhs, rhs, ans);
    return ans;
}

std::ostream& operator<<(std::ostream& out, const Matrix& mat)
{
    for (size_t i = 0; i < mat.rows(); i++) {
        for (size_t j = 0; j < mat.cols(); j++) {
            out << mat.get(i, j) << " ";
        }
        out << std::endl;
    }

    return out;
}

void row_major_mul()
{
    MatrixRowMajor m1(N, N);
    m1.fill_with_test_data();

    MatrixRowMajor m2(N, N);
    m2.fill_with_test_data();

    auto m3 = m1 * m2;
    std::cout << m3.sum();
}

void row_col_major_mul()
{
    MatrixRowMajor m1(N, N);
    m1.fill_with_test_data();

    MatrixColMajor m2(N, N);
    m2.fill_with_test_data();

    auto m3 = m1 * m2;
    std::cout << m3.sum();
}

void col_row_major_mul()
{
    MatrixRowMajor m1(N, N);
    m1.fill_with_test_data();

    MatrixColMajor m2(N, N);
    m2.fill_with_test_data();

    auto m3 = m2 * m1;
    std::cout << m3.sum();
}

int main(int argc, char* argv[])
{
    row_major_mul();
    row_col_major_mul();
    col_row_major_mul();
    return 0;
}