small matrix test https://godbolt.org/z/tctouH

matmul.cpp

//https://godbolt.org/z/tctouH

#include <iostream>
#include <random>
#include <functional>
#include <exception>
#include <cstring>
#include <chrono>

using namespace std;

#define MUL(other, i, j, k) m[i][k] * other.m[k][j]

#define SUM(other, i, j, dst) dst.m[i][j] = MUL(other, i, j, 0) + MUL(other, i, j, 1) + MUL(other, i, j, 2) + MUL(other, i, j, 3) ;   

#define MUL_ROW(other, r, dst) SUM(other, r, 0, dst) SUM(other, r, 1, dst) SUM(other, r, 2, dst) SUM(other, r, 3, dst)


#define SUM_ROW(r) m[r][0] + m[r][1] + m[r][2] + m[r][3] 

#define RANDOM_ROW(r, gen) m[r][0] = gen(); m[r][1] = gen(); m[r][2] = gen(); m[r][3] = gen();


struct Matrix4x4
{
    float m[4][4];

    Matrix4x4() {memset(m, 0, sizeof(m));}    

    void rand() {
        auto gen = std::bind(uniform_real_distribution<float>(0.f,1.f), mt19937{random_device{}()});
        
        // for (int i = 0; i < 4; ++i)
        //     for(int j = 0; j < 4; ++j)
        //         m[i][j] = gen();
        
        // return;
        
        RANDOM_ROW(0, gen)
        RANDOM_ROW(1, gen)
        RANDOM_ROW(2, gen)
        RANDOM_ROW(3, gen)
    }

    //1.5x slower
    Matrix4x4 mul_v1 (const Matrix4x4& other) const
    {
        Matrix4x4 res;
        for (int i = 0; i < 4; ++i)
            for(int j = 0; j < 4; ++j){
                for(int k = 0; k < 4; ++k){
                    res.m[i][j] += m[i][k] * other.m[k][j];                                    
                }
        }
        return res;
    }
    
   Matrix4x4 operator * (const Matrix4x4& other) const
    {
        //return mul_v1(other);
        
        Matrix4x4 res;
        MUL_ROW(other, 0, res)
        MUL_ROW(other, 1, res)        
        MUL_ROW(other, 2, res)        
        MUL_ROW(other, 3, res)
        return res;
        
    }

    inline float& at(int i, int j){
        return m[i][j];
    }    


    float sum() const{
        return SUM_ROW(0) + SUM_ROW(1) + SUM_ROW(2) + SUM_ROW(3);
    }

};

#undef MUL
#undef SUM
#undef SUM_ROW
#undef MUL_ROW
#undef RANDOM_ROW


#define DECLARE_ROW(r) float m ## r ## 0 = 0, m ## r ## 1 = 0, m ## r ## 2 = 0, m ## r ## 3 = 0;  

#define RANDOM_ROW(r, gen) m ## r ## 0 = gen(); m ## r ## 1 = gen(); m ## r ## 2 = gen(); m ## r ## 3 = gen();


#define MUL(other, i, j, k) m ## i ## k * other.m ## k ## j

#define SUM(other, i, j, dst) dst.m ## i ## j = MUL(other, i, j, 0) + MUL(other, i, j, 1) + MUL(other, i, j, 2) + MUL(other, i, j, 3) ;   

#define MUL_ROW(other, r, dst) SUM(other, r, 0, dst) SUM(other, r, 1, dst) SUM(other, r, 2, dst) SUM(other, r, 3, dst)

#define SUM_ROW(r) m ## r ## 0 + m ## r ## 1 + m ## r ## 2 + m ## r ## 3 

#define CASE_IN_ROW(r, selector) { switch(selector) { \
case 0: return m ## r ## 0; \
case 1: return m ## r ## 1; \
case 2: return m ## r ## 2; \
case 3: return m ## r ## 3; \
} }

#define CASE(i, j) switch(i) {\
case 0: CASE_IN_ROW(0, j) \
case 1: CASE_IN_ROW(1, j) \
case 2: CASE_IN_ROW(2, j) \
case 3: CASE_IN_ROW(3, j) \
}

struct Matrix4x4m
{
    DECLARE_ROW(0)
    DECLARE_ROW(1)  
    DECLARE_ROW(2)
    DECLARE_ROW(3)
    

    void rand()
    {
        auto gen = std::bind(uniform_real_distribution<float>(0.f,1.f), mt19937{random_device{}()});
        RANDOM_ROW(0, gen)
        RANDOM_ROW(1, gen)
        RANDOM_ROW(2, gen)
        RANDOM_ROW(3, gen)
    }

    Matrix4x4m operator * (const Matrix4x4m& other) const
    {
        Matrix4x4m res;
        MUL_ROW(other, 0, res)
        MUL_ROW(other, 1, res)        
        MUL_ROW(other, 2, res)        
        MUL_ROW(other, 3, res)
        return res;
    }

    float& at(int i, int j)
    {
        CASE(i, j)
        
        throw range_error("range error");    
    };

    float sum() const{
        return SUM_ROW(0) + SUM_ROW(1) + SUM_ROW(2) + SUM_ROW(3);
    }

};


template <typename MatT>
void test(volatile int rand_step, volatile int rand_pos)
{
    MatT first;
    
    first.rand();

    MatT second;
    second.rand();

    float sum = 0;
    

    auto start = chrono::high_resolution_clock::now();

    for (size_t i = 1; i < 10000000u; ++i){

        if (i % rand_step == rand_pos){
            first.rand();
            second.rand();
        }
        
        auto mul = first * second;
        sum += mul.sum();
    }

    auto elapsed = chrono::high_resolution_clock::now() - start;
   
    cout << "sum: " << sum << endl;

    cout << "elapsed " << chrono::duration_cast<chrono::milliseconds>(elapsed).count() << endl;
}

int main()
{
    
    int rand_step = 500 + rand() % 1000;

    int rand_pos = rand() % rand_step;

    cout << "array:\n";
    test<Matrix4x4>(rand_step, rand_pos);

    
    cout << "members:\n";
    test<Matrix4x4m>(rand_step, rand_pos);

    
    return 0;
}