利用模板偏特化为特定尺寸矩阵实现循环展开

Hey there! I totally get the frustration here—subclassing your matrix template for specific dimensions like 2x2 leads to messy, redundant code, which defeats the purpose of using templates in the first place. Let’s walk through two clean, efficient solutions to implement loop-unrolled operator*= for your 2x2 matrices without all the extra baggage.

Option 1: Full Template Specialization (Works for All C++ Versions) #

If you need compatibility with older C++ standards (pre-C++17), full template specialization is the way to go. We’ll explicitly define the 2x2 matrix class with its unrolled multiplication logic, while keeping the generic template for all other sizes.

Here’s the updated code:

#include <iostream>

// Generic matrix template for all sizes
template<class T, size_t M, size_t N>
class matrix {
public:
    T data[M][N];

    // Default constructor
    matrix() = default;
    // Convenience constructor for easy initialization
    template<typename... Args>
    matrix(Args&&... args) : data{std::forward<Args>(args)...} {}

    matrix<T,M,N>& operator*= (const matrix<T,M,N> &B);
};

// Generic matrix multiplication implementation
template<class T, size_t M, size_t N>
matrix<T,M,N>& matrix<T,M,N>::operator*= (const matrix<T,M,N> &B) {
    // Standard nested loop for general cases
    for (size_t i = 0; i < M; ++i) {
        for (size_t j = 0; j < N; ++j) {
            T sum = T{};
            for (size_t k = 0; k < N; ++k) {
                sum += data[i][k] * B.data[k][j];
            }
            data[i][j] = sum;
        }
    }
    return *this;
}

// Full specialization for 2x2 matrices
template<class T>
class matrix<T, 2, 2> {
public:
    T data[2][2];

    matrix() = default;
    matrix(T a, T b, T c, T d) : data{{a, b}, {c, d}} {}

    matrix<T,2,2>& operator*= (const matrix<T,2,2> &B);
};

// Unrolled multiplication for 2x2 matrices
template<class T>
matrix<T,2,2>& matrix<T,2,2>::operator*= (const matrix<T,2,2> &B) {
    // No loops here—fully unrolled for maximum performance
    const T a = data[0][0] * B.data[0][0] + data[0][1] * B.data[1][0];
    const T b = data[0][0] * B.data[0][1] + data[0][1] * B.data[1][1];
    const T c = data[1][0] * B.data[0][0] + data[1][1] * B.data[1][0];
    const T d = data[1][0] * B.data[0][1] + data[1][1] * B.data[1][1];

    data[0][0] = a;
    data[0][1] = b;
    data[1][0] = c;
    data[1][1] = d;

    return *this;
}

int main() {
    // Test the 2x2 unrolled multiplication
    matrix<double,2,2> m1(1.0, 2.0, 3.0, 4.0);
    matrix<double,2,2> m2(5.0, 6.0, 7.0, 8.0);
    m1 *= m2;

    std::cout << "2x2 Result:\n";
    std::cout << m1.data[0][0] << " " << m1.data[0][1] << "\n";
    std::cout << m1.data[1][0] << " " << m1.data[1][1] << "\n";

    // Test a generic size (3x3) to confirm it uses the generic implementation
    matrix<int,3,3> m3{{1,0,0}, {0,1,0}, {0,0,1}};
    m3 *= m3; // Identity matrix squared stays identity
    return 0;
}

Pros of This Approach: #

• Backward compatible: Works with C++98 and later.
• Clean separation: The 2x2 logic is isolated, making it easy to modify or debug.
• Max optimization: The compiler can fully optimize the unrolled code with no loop overhead.

Option 2: if constexpr for Inline Specialization (C++17+) #

If you’re using C++17 or newer, if constexpr is a far more concise solution. It lets you embed the unrolled 2x2 logic directly in the generic operator*= method—compiler will eliminate the unused branch at compile time, so there’s no runtime cost.

Here’s how it looks:

#include <iostream>

template<class T, size_t M, size_t N>
class matrix {
public:
    T data[M][N];

    matrix() = default;
    template<typename... Args>
    matrix(Args&&... args) : data{std::forward<Args>(args)...} {}

    matrix<T,M,N>& operator*= (const matrix<T,M,N> &B);
};

template<class T, size_t M, size_t N>
matrix<T,M,N>& matrix<T,M,N>::operator*= (const matrix<T,M,N> &B) {
    if constexpr (M == 2 && N == 2) {
        // Unrolled 2x2 multiplication—only compiled for 2x2 matrices
        const T a = data[0][0] * B.data[0][0] + data[0][1] * B.data[1][0];
        const T b = data[0][0] * B.data[0][1] + data[0][1] * B.data[1][1];
        const T c = data[1][0] * B.data[0][0] + data[1][1] * B.data[1][0];
        const T d = data[1][0] * B.data[0][1] + data[1][1] * B.data[1][1];

        data[0][0] = a;
        data[0][1] = b;
        data[1][0] = c;
        data[1][1] = d;
    } else {
        // Generic loop-based multiplication for all other sizes
        for (size_t i = 0; i < M; ++i) {
            for (size_t j = 0; j < N; ++j) {
                T sum = T{};
                for (size_t k = 0; k < N; ++k) {
                    sum += data[i][k] * B.data[k][j];
                }
                data[i][j] = sum;
            }
        }
    }
    return *this;
}

int main() {
    matrix<int,2,2> m1{{1,2}, {3,4}};
    matrix<int,2,2> m2{{5,6}, {7,8}};
    m1 *= m2;

    std::cout << "2x2 Result:\n";
    std::cout << m1.data[0][0] << " " << m1.data[0][1] << "\n";
    std::cout << m1.data[1][0] << " " << m1.data[1][1] << "\n";

    // Test 3x3 matrix to verify generic logic works
    matrix<float,3,3> m4;
    // Initialize and test as needed
    return 0;
}

Pros of This Approach: #

• No redundant code: All logic lives in one class—no need to redefine constructors or members for specializations.
• Compile-time efficiency: The unused branch is completely stripped out during compilation, so no runtime overhead.
• Scalable: Easy to add more special cases (like 3x3 or 4x4) by adding more if constexpr blocks.

Quick Optimization Tips #

• Use contiguous storage: If you switch to a 1D array (T data[M*N]), the compiler can optimize memory access patterns better, pairing nicely with loop unrolling.
• Mark methods as inline: Putting the unrolled operator*= directly in the class body (or using inline) helps the compiler inline the code, further boosting performance.
• Test with your target types: For floating-point types, some compilers might auto-unroll small loops, but explicit unrolling guarantees the optimization regardless of compiler settings.

内容的提问来源于stack exchange，提问作者Tom de Geus