Multiplication of 3×3 and N×N square matrices I recently wrote a matrix module in C++.

During the development process I quoted some source code and found a problem.

For example, matrix multiplication:

This way is suitable for all N×N matrices：

void multiplyMatrix(const float32* a, const float32* b, float32* dst,
                    int32 aColumns, int32 bColumns, int32 dstColumns, int32 dstRows) {

    for (int32 i = 0; i &lt; dstRows; i++) {

        for (int32 j = 0; j &lt; dstColumns; j++) 
            dst[i * dstColumns + j] = dotMatrix(a, b, aColumns, bColumns, j, i);

    }
}

float32 dotMatrix(const float32* a, const float32* b, 
                         int32 aColumns, int32 bColumns, 
                         int32 column, int32 row) {

    float32 result = 0.0f;

    int32 index = aColumns * row;
    for (int32 i = 0; i &lt; aColumns; i++) {
        result += a[index++] * b[column];
        column += bColumns;
    }

    return result;
}


Next, I wrote a 3x3 matrix class.

class Matrix3x3
{
public:
float32 m11, m12, m13, 
        m21, m22, m23, 
        m31, m32, m33;

float32 element[9];

void multiply(float32 ma11, float32 ma12, float32 ma13,
              float32 ma21, float32 ma22, float32 ma23,
              float32 ma31, float32 ma32, float32 ma33) {

    float32 temp1 = m11 * ma11 + m21 * ma12 + m31 * ma13;
    float32 temp2 = m12 * ma11 + m22 * ma12 + m32 * ma13;

    m13 = m13 * ma11 + m23 * ma12 + m33 * ma13;
    m11 = temp1;
    m12 = temp2;

    temp1 = m11 * ma21 + m21 * ma22 + m31 * ma23;
    temp2 = m12 * ma21 + m22 * ma22 + m32 * ma23;
    m23   = m13 * ma21 + m23 * ma22 + m33 * ma23;
    m21   = temp1;
    m22   = temp2;

    temp1 = m11 * ma31 + m21 * ma32 + m31 * ma33;
    temp2 = m12 * ma31 + m22 * ma32 + m32 * ma33;
    m31   = m13 * ma31 + m23 * ma32 + m33 * ma33;
    m32   = temp1;
    m33   = temp2;
}
}


Obviously the first one is very convenient.

Next, I tested the time it took to calculate：

        float32 e1[9];
        e1[0] = 2.1018f;   e1[1] = -1.81754f; e1[2] = 1.2541f;
        e1[3] = 0.54194f;  e1[4] = 2.75391f;  e1[5] = -0.1167f;
        e1[6] = -5.81652f; e1[7] = -7.9381f;  e1[8] = 4.2816f;

        float32 e2[9];
        e2[0] = 2.1018f;   e2[1] = -1.81754f; e2[2] = 1.2541f;
        e2[3] = 0.54194f;  e2[4] = 2.75391f;  e2[5] = -0.1167f;
        e2[6] = -5.81652f; e2[7] = -7.9381f;  e2[8] = 4.2816f;

        Matrix3x3 a;
        a.m11 = 2.1018f;   a.m12 = -1.81754f; a.m13 = 1.2541f; 
        a.m21 = 0.54194f;  a.m22 = 2.75391f;  a.m23 = -0.1167f;
        a.m31 = -5.81652f; a.m32 = -7.9381f;  a.m33 = 4.2816f;

        Matrix3x3 b = a;

        float64 timeSpent = 0;
        LARGE_INTEGER nFreq;
        LARGE_INTEGER nBeginTime;
        LARGE_INTEGER nEndTime;

        QueryPerformanceFrequency(&amp;nFreq); // statistical frequency
        QueryPerformanceCounter(&amp;nBeginTime);// start timer

        for (int32 i = 0; i &lt; 100000; i++) {
            multiplyMatrix(e1, e2, dst, 3, 3, 3, 3);
        }

    QueryPerformanceCounter(&amp;nEndTime); //end timer
            timeSpent = (float64)(nEndTime.QuadPart - nBeginTime.QuadPart) / (nFreq.QuadPart);

   printf("timeSpent1:%f\n", timeSpent);

   QueryPerformanceCounter(&amp;nBeginTime);
   for (int32 i = 0; i &lt; 100000; i++) {
            b.multiply(a.m11, a.m12, a.m13, 
                       a.m21, a.m22, a.m23, 
                       a.m31, a.m32, a.m33);
        }
   QueryPerformanceCounter(&amp;nEndTime);
        timeSpent = (float64)(nEndTime.QuadPart - nBeginTime.QuadPart) / (nFreq.QuadPart);

printf("timeSpent2:%f\n", timeSpent);


Output：

timeSpent1:0.014277
timeSpent2:0.004649
timeSpent1:0.012684
timeSpent2:0.004522
.......
.......
timeSpent1:0.003414
timeSpent2:0.001166
timeSpent1:0.003407
timeSpent2:0.001242


Is this difference in efficiency significant or negligible?
A factor of 3 is large, but in my opinion not unexpected or abnormal. The functions that can handle a variable size matrix in their natural form (ie as they would be compiled without knowledge of the size, for example if the functions are defined in a different compilation unit than they are used in and LTO is not applied) have a lot of overhead: non-linear control flow (3 nested loops), more complicated address computation (involving multiplication by a variable).

Basically, that is the cost of generality .. but there is more to it.

From your use of QueryPerformanceCounter I assume you use MSVC (other compilers aren't much different for the following considerations). MSVC likes to unroll loops such as the one in dotMatrix by 4. It does not like to unroll such loops by 3, though it can be persuaded to do so anyway, for example by giving it a loop that makes exactly 3 iterations. So the cost of generality would work out much differently if the relevant matrix was of size 4x4 or 8x8, as in those cases only the faster unrolled codepath would be used (this still comes with overhead, but less). 3 is a bad case, only ever using the fallback codepath.

Additionally, the general matrix multiply implemented by multiplyMatrix is not scalable: it does not implement cache blocking, so for any matrix that does not fit in L1 cache it will perform badly (and even more badly when going beyond the L2 and L3 sizes). That is normal for code in general, but matrix multiplication is special in that it does not have to suffer significantly from that common effect thanks to its "O(n2) data in O(n3) time" property.

Both the general matrix multiply and the special 3x3 one could use SIMD intrinsics for extra efficiency. 3x3 is an awkward size that would cause some "wasted lanes", but it would still help. For example, it could be done like this (not tested):

#include &lt;xmmintrin.h&gt;

class Matrix3x3
{
public:
float32 m11, m21, m31, 
        m12, m22, m32, 
        m13, m23, m33, padding;

void multiply(float32 ma11, float32 ma12, float32 ma13,
              float32 ma21, float32 ma22, float32 ma23,
              float32 ma31, float32 ma32, float32 ma33) {

    __m128 col1 = _mm_loadu_ps(&amp;m11);
    __m128 col2 = _mm_loadu_ps(&amp;m12);
    __m128 col3 = _mm_loadu_ps(&amp;m13);
    __m128 t1 = _mm_add_ps(_mm_add_ps(
        _mm_mul_ps(col1, _mm_set1_ps(ma11)),
        _mm_mul_ps(col2, _mm_set1_ps(ma21))),
        _mm_mul_ps(col3, _mm_set1_ps(ma31)));
    __m128 t2 = _mm_add_ps(_mm_add_ps(
        _mm_mul_ps(col1, _mm_set1_ps(ma12)),
        _mm_mul_ps(col2, _mm_set1_ps(ma22))),
        _mm_mul_ps(col3, _mm_set1_ps(ma32)));
    __m128 t3 = _mm_add_ps(_mm_add_ps(
        _mm_mul_ps(col1, _mm_set1_ps(ma13)),
        _mm_mul_ps(col2, _mm_set1_ps(ma23))),
        _mm_mul_ps(col3, _mm_set1_ps(ma33)));

    _mm_storeu_ps(&amp;m11, t1);
    _mm_storeu_ps(&amp;m12, t2);
    _mm_storeu_ps(&amp;m13, t3);
}
};


The padding is a bit unfortunate (and shouldn't be private, because that makes its positioning relative to the actual matrix elements undefined), but simplifies the SIMD logic, chunks of 16 bytes are easier to deal with. It is possible to avoid the padding if required. Anyway, this results in a significant reduction in code and should be more efficient (without AVX the set1s cost more, that shouldn't be enough to undo the improvement but I didn't try it). The dllexport in the code on godbolt is not really part of the code, I just put that there to force code to be generated for an otherwise unused method.

Column-major order is used here because the columns of the result are a linear combination of the columns of the left hand matrix, which we have access to in packed memory. Similarly, the rows of the output are a linear combination of the rows of the right hand side, but we have no packed access to the rows of the right hand side, so they would be inefficient to gather. A row-oriented version of the above could be arranged for example if the right hand side was passed in as a reference to a Matrix3x3.

Passing the right hand side as matrix is probably a nicer interface anyway, with 9 separate arguments there is no choice but to write them all out separately even if the RHS is available as a matrix object, as you already experienced in your benchmark code.