Multiplication Of A Matrix In C++

Matrix multiplication stands as a cornerstone operation in linear algebra, serving as the computational backbone for fields ranging from computer graphics and machine learning to physics simulations and cryptography. In real terms, when implementing this operation in C++, developers gain precise control over memory management and execution speed, making it the language of choice for high-performance numerical computing. Understanding how to translate the mathematical definition into efficient, readable C++ code requires a grasp of both the algorithmic logic and the language features that optimize data access patterns.

The Mathematical Foundation

Before diving into syntax, Visualize the rule governing the operation — this one isn't optional. Which means given two matrices, A with dimensions m x n (m rows, n columns) and B with dimensions n x p, the multiplication C = A * B is only valid if the number of columns in A equals the number of rows in B. The resulting matrix C will possess dimensions m x p.

Each element c[i][j] in the result matrix is computed as the dot product of the i-th row of matrix A and the j-th column of matrix B. Mathematically, this is expressed as:

$C_{ij} = \sum_{k=1}^{n} A_{ik} \times B_{kj}$

This triple-nested summation structure—iterating over rows of A, columns of B, and the shared dimension k—forms the basis of the standard naive implementation.

Setting Up the C++ Environment

A clean implementation begins with the necessary headers. For matrix operations, <iostream> handles input and output, <vector> provides dynamic, bounds-checked arrays (preferred over raw pointers for safety), and <iomanip> assists with formatted output alignment.

#include 
#include 
#include 

// Using namespace std for brevity in examples
using namespace std;

Using std::vector<std::vector<double>> (or int, float depending on precision needs) represents a jagged array or vector of vectors. While convenient, this structure stores rows in potentially non-contiguous memory blocks. For high-performance scenarios, a flattened 1D vector (std::vector<double> of size rows * cols) is superior due to cache locality, but the 2D vector approach remains the standard pedagogical starting point.

And yeah — that's actually more nuanced than it sounds.

The Naive Implementation: Triple Nested Loops

The most direct translation of the mathematical formula involves three nested loops. The outer loop iterates through rows of A (index i), the middle loop iterates through columns of B (index j), and the inner loop performs the summation over the shared dimension k.

vector> multiplyMatrices(const vector>& A, const vector>& B) {
    int rowsA = A.size();
    int colsA = A[0].size();
    int rowsB = B.size();
    int colsB = B[0].size();

    // Validation: Columns of A must equal Rows of B
    if (colsA != rowsB) {
        cerr << "Error: Incompatible matrix dimensions for multiplication." << endl;
        return {}; // Return empty matrix
    }

    // Initialize Result Matrix C with dimensions rowsA x colsB, filled with 0.0
    vector> C(rowsA, vector(colsB, 0.0));

    // The Triple Loop
    for (int i = 0; i < rowsA; ++i) {
        for (int j = 0; j < colsB; ++j) {
            double sum = 0.0;
            for (int k = 0; k < colsA; ++k) {
                sum += A[i][k] * B[k][j];
            }
            C[i][j] = sum;
        }
    }
    return C;
}

Key Observations:

1. Const Correctness: The input matrices are passed as const vector<vector<double>>& to avoid expensive deep copies.
2. Initialization: C is initialized with zeros immediately. Accumulating into a local variable sum inside the j loop (rather than C[i][j] += ...) reduces repeated memory writes to the result matrix, a minor but measurable optimization known as scalar replacement.
3. Validation: Checking colsA != rowsB prevents undefined behavior or crashes.

Optimizing for Cache Performance: Loop Interchange

The naive implementation above accesses B[k][j] in the innermost loop. Still, since C++ stores vectors in row-major order, elements in the same row are contiguous in memory. Accessing B[k][j] jumps across different rows (strided access), causing frequent cache misses.

A classic optimization is Loop Interchanging (specifically swapping the j and k loops). This transforms the access pattern to be sequential for both A and B.

vector> multiplyOptimized(const vector>& A, const vector>& B) {
    int rowsA = A.size();
    int colsA = A[0].size();
    int colsB = B[0].size();

    if (colsA != (int)B.size()) return {};

    vector> C(rowsA, vector(colsB, 0.0));

    // Loop Order: i -> k -> j
    for (int i = 0; i < rowsA; ++i) {
        for (int k = 0; k < colsA; ++k) {
            double a_ik = A[i][k]; // Load once per k iteration
            for (int j = 0; j < colsB; ++j) {
                // Sequential access: C[i][j] and B[k][j] are row-contiguous
                C[i][j] += a_ik * B[k][j];
            }
        }
    }
    return C;
}

Why this is faster:

• A[i][k] is invariant in the inner loop; hoisting it to a register (a_ik) eliminates repeated loads.
• C[i][j] writes sequentially across memory.
• B[k][j] reads sequentially across memory. This arrangement exploits spatial locality, allowing the CPU prefetcher to pull entire cache lines, drastically reducing latency for large matrices.

The Flattened 1D Array Approach

For production-grade numerical libraries (like Eigen, Blaze, or custom high-performance kernels), the vector<vector<T>> abstraction is abandoned in favor of a single contiguous block of memory. This eliminates the double indirection (pointer chasing) inherent in vectors of vectors and guarantees perfect cache locality for row access Not complicated — just consistent..

// Flattened representation: index = row * num_cols + col
using Matrix = vector;

Matrix multiplyFlattened(const Matrix& A, const Matrix& B, int rowsA, int colsA, int colsB) {
    if (colsA != (int)B.size() / colsB) return {}; // Simplified check

    Matrix C(rowsA * colsB, 0.0);

    for (int i = 0; i < rowsA; ++i) {
        for (int k = 0; k < colsA; ++k) {
            double a_ik = A[i * colsA + k];
            int offsetB = k * colsB;
            int offsetC = i * colsB;
            for (int j = 0; j < colsB; ++j) {
                C[offsetC + j] += a_ik * B[offsetB + j];
            }
        }
    }
    return C;
}

This manual index calculation (row * stride + col) is verbose but unlocks the ability to use **SIMD (Single

Instruction Multiple Data) vectorization. With contiguous memory and predictable strides, the compiler can auto-vectorize the inner j loop (using #pragma omp simd or -O3 flags), or developers can write explicit intrinsics (AVX2, AVX-512, NEON) to process 4, 8, or 16 double elements per instruction. The inner loop effectively becomes a fused multiply-add (FMA) stream:

// Conceptual AVX2 intrinsic style (processing 4 doubles at a time)
for (int j = 0; j < colsB; j += 4) {
    __m256d a_vec = _mm256_set1_pd(a_ik);           // Broadcast scalar
    __m256d b_vec = _mm256_loadu_pd(&B[offsetB + j]); // Load 4 contiguous doubles
    __m256d c_vec = _mm256_loadu_pd(&C[offsetC + j]); // Load 4 contiguous doubles
    c_vec = _mm256_fmadd_pd(a_vec, b_vec, c_vec);     // C += A * B (FMA)
    _mm256_storeu_pd(&C[offsetC + j], c_vec);         // Store results
}

Handling remainder elements (when colsB % 4 != 0) requires a scalar cleanup loop or masked loads (AVX-512).

Cache Blocking (Tiling): Conquering the Memory Wall

Loop interchange and SIMD optimize the core compute, but for large matrices ($N > \text{L1/L2 cache size}$), the working set exceeds cache capacity. On top of that, data loaded for the first k-iteration is evicted before it can be reused for subsequent i or j iterations. Cache Blocking (Tiling) solves this by restructuring the computation to operate on small tiles (blocks) that fit entirely within the L1/L2 cache Still holds up..

We partition matrices $A$, $B$, and $C$ into sub-matrices (tiles) of size $T \times T$ (e.Which means g. So , $64 \times 64$ or $128 \times 128$). The algorithm computes $C_{tile} = \sum A_{tile} \times B_{tile}$, keeping the active tiles resident in fast cache.

void multiplyBlocked(const Matrix& A, const Matrix& B, Matrix& C, 
                     int rowsA, int colsA, int colsB, int blockSize) {
    // Initialize C to zero (assumed done prior or here)
    // std::fill(C.begin(), C.end(), 0.0);

    for (int i0 = 0; i0 < rowsA; i0 += blockSize) {
        for (int k0 = 0; k0 < colsA; k0 += blockSize) {
            // Prefetch/Load Block A[i0:i0+bs, k0:k0+bs] into L1/L2
            for (int j0 = 0; j0 < colsB; j0 += blockSize) {
                
                // Micro-kernel: Compute C_block += A_block * B_block
                // Limits handle edge cases where matrix dim % blockSize != 0
                int iMax = std::min(i0 + blockSize, rowsA);
                int kMax = std::min(k0 + blockSize, colsA);
                int jMax = std::min(j0 + blockSize, colsB);

                for (int i = i0; i < iMax; ++i) {
                    for (int k = k0; k < kMax; ++k) {
                        double a_ik = A[i * colsA + k];
                        int offsetB = k * colsB + j0;
                        int offsetC = i * colsB + j0;
                        // Inner loop (j) is now tiny (blockSize), perfect for SIMD unrolling
                        for (int j = j0; j < jMax; ++j) {
                            C[offsetC + (j - j0)] += a_ik * B[offsetB + (j - j0)];
                        }
                    }
                }
            }
        }
    }
}

Why Tiling Wins:

1. Data Reuse: A tile of $A$ (size $T \times T$) is loaded once and reused for $T$ columns of $B$. A tile of $B$ is loaded once and reused for $T$ rows of $A$. The $C$ tile stays in registers/L1 for the entire $k_0$ iteration.
2. Cache Fit: Choosing blockSize such that $3 \times T^2 \times \text{sizeof(double)} \le \text{L1 Cache Size}$ (typically 32KB–48KB) ensures zero capacity misses inside the micro-kernel.
3. TLB Efficiency: Accessing