G06: Scalable architecture and measured performance

Prerequisites and observable goals#

The central question: which resource bounds this workload?#

A program uses instruction throughput, memory bandwidth, cache capacity, synchronization, allocation, and I/O. Improving a nonlimiting resource may not improve latency. Start with the specified output and input distribution, then count data movement and arithmetic. Treat a cost model as a prediction to test, not as a replacement for measurement.

For dense $C=AB$ with $A$ of shape $m\times k$ and $B$ of shape $k\times n$, the conventional algorithm computes

There are $mkn$ multiplications and $mn(k-1)$ additions when each dot product is accumulated from its first product, often summarized as approximately $2mkn$ floating operations. A fused multiply-add may count as two arithmetic operations in a performance convention while using one instruction; state the convention. The arithmetic expression does not tell you how often a cache line is transferred.

Three related bounds#

For fixed total work with serial fraction $s$ and perfectly scaling remainder, ideal speedup with $p$ workers is

This is Amdahl’s model. It ignores overhead, contention, and workload changes. For a computation DAG, total work $W$ and longest dependency-chain work $D$ (span) give $T_p\ge\max(W/p,D)$. Even unlimited processors cannot execute dependent operations simultaneously.

For a chosen memory level, define arithmetic intensity $I$ as floating operations divided by bytes moved through that level. If compute throughput is bounded by $P_{\max}$ and bandwidth by $B$, then a simple roofline-style bound is

The units check: $(\mathrm{bytes/s})(\mathrm{operations/byte})=\mathrm{operations/s}$. This is an upper envelope, not a promise that a kernel achieves it. Cache misses, dependencies, issue limits, and instruction mix can reduce achieved throughput further. Choose the memory level explicitly: DRAM traffic and L1 traffic are different quantities.

Why loop order changes locality#

Assume row-major contiguous arrays. In the $i,j,p$ order, consecutive $p$ reads from $A$ are contiguous but reads from $B_{pj}$ stride by an entire row. In $i,p,j$ order, hold one $A_{ip}$ value while sweeping contiguous entries of $B$ and $C$. The formulas describe the same real-number result, but floating-point rounding and vectorization choices must still be checked.

Blocking processes tiles so parts of $A$, $B$, and $C$ can be reused before eviction. A rough working-set estimate for three $b\times b$ tiles of doubles is $3b^2(8)$ bytes. For $b=24$, that is 13,824 bytes, ignoring metadata, associativity, prefetching, other active data, and exact loop structure. A cache of that nominal size does not guarantee the tiles fit effectively. Benchmark multiple justified candidates; do not select the tile solely because one capacity division looks neat.

This is a schematic tile picture, not a diagram of actual cache contents. The selected $C$ block accumulates contributions from successive matching inner-dimension blocks of $A$ and $B$. Edge blocks can be smaller; treating every dimension as a multiple of the tile size is a correctness bug.

Parallelism without ownership confusion#

For parallel matrix output, assign disjoint output tiles to workers and let each worker accumulate the complete inner sum for its tile. Sharing a single output element across workers requires a reduction or synchronization; an unsynchronized += is not made safe by using doubles. Partitioning by the inner dimension produces partial matrices that need a defined merge. Partitioning by output rows simplifies ownership but may interact differently with scheduling and locality.

False sharing occurs when workers update distinct data that share a coherence unit, usually a cache line. The objects can be logically independent yet still cause hardware ownership traffic. Padding is not an automatic fix: it consumes memory and can harm locality. Confirm with measurements and the actual hardware’s line size and access pattern before asserting a diagnosis.

SIMD executes an operation across lanes; threads schedule independent instruction streams; GPU kernels expose another parallel execution and memory model. These are related tools, not interchangeable meanings of “parallel.” A small host-to-device transfer plus launch overhead can exceed the time saved by a GPU kernel. No GPU execution is claimed by this CPU lab.

A benchmark is an experiment with an acceptance test#

The companion benchmarks the two C++ implementations at $n=192$, tile 24, with a seeded input, a warm-up pair, and seven measured pairs whose method order alternates. Both functions include input validation and output allocation in the timed region. Correctness comparison happens after each timed call. The reference result is computed before the measurement loop; a checksum is consumed so the output is not simply ignored.

This is one local environment and workload, not a claim about a consumer GPU such as an RTX 5090, about Apple hardware, or about all matrices. Do not divide one printed timing by another and call that a universal acceleration factor. Examine paired ratios/differences and spread, record build settings, and repeat on your target. The CSV is included with the actual environment record. Unchanged operation count does not contradict a measured improvement from reduced memory traffic.

Production judgment beyond a fast kernel#

A production library must also consider overflow in sizes, invalid dimensions, allocation failure, thread lifecycle, ABI stability, numerical tolerances, cancellation policy, portability, and performance across many shapes. The lab validates small bounded dimensions so shape products remain representable, but it is not a replacement for a tuned BLAS library. A correct educational implementation can be useful precisely because it reveals the mechanics before a library is introduced.

Independent practice#

F foundation C standard course work T transfer S optional stretch L runnable lab or project. Attempt the complete set before opening its solutions.

Smallest useful hints#

For .02, convert all rate prefixes before comparing. For .05, identifying a unique writer for each output is often easier than trying to lock every addition. For .06, cap each block end with the actual dimension.

Complete solutions: G06#

Compare the first step where your work differs. A corrected answer is not yet an independent success.

Open every solution

#pragma once
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <vector>
namespace graduate {
inline void validate(const std::vector<double>& a,
                     const std::vector<double>& b,
                     std::size_t m, std::size_t k, std::size_t n) {
    if (m==0 || k==0 || n==0 || m>512 || k>512 || n>512)
        throw std::invalid_argument("dimensions must be in 1..512");
    if (a.size()!=m*k || b.size()!=k*n)
        throw std::invalid_argument("shape mismatch");
    for (double x:a) if (!std::isfinite(x))
        throw std::invalid_argument("nonfinite A");
    for (double x:b) if (!std::isfinite(x))
        throw std::invalid_argument("nonfinite B");
}
inline void check_result(const std::vector<double>& c) {
    for (double x:c) if (!std::isfinite(x))
        throw std::overflow_error("nonfinite result");
}
inline std::vector<double> naive(const std::vector<double>& a,
                                 const std::vector<double>& b,
                                 std::size_t m,std::size_t k,std::size_t n) {
    validate(a,b,m,k,n);
    std::vector<double> c(m*n,0.);
    for (std::size_t i=0;i<m;++i)
        for (std::size_t j=0;j<n;++j)
            for (std::size_t p=0;p<k;++p)
                c[i*n+j]+=a[i*k+p]*b[p*n+j];
    check_result(c);
    return c;
}
inline std::vector<double> blocked(const std::vector<double>& a,
                                   const std::vector<double>& b,
                                   std::size_t m,std::size_t k,std::size_t n,
                                   std::size_t tile=24) {
    validate(a,b,m,k,n);
    if (tile==0 || tile>1024)
        throw std::invalid_argument("tile must be in 1..1024");
    std::vector<double> c(m*n,0.);
    for (std::size_t ii=0;ii<m;ii+=tile)
        for (std::size_t pp=0;pp<k;pp+=tile)
            for (std::size_t jj=0;jj<n;jj+=tile)
                for (std::size_t i=ii;i<std::min(ii+tile,m);++i)
                    for (std::size_t p=pp;p<std::min(pp+tile,k);++p) {
                        const double value=a[i*k+p];
                        for (std::size_t j=jj;j<std::min(jj+tile,n);++j)
                            c[i*n+j]+=value*b[p*n+j];
                    }
    check_result(c);
    return c;
}
} // namespace graduate

Benchmark driver and exact timing boundary#

#include <matrix.hpp>
#include <chrono>
#include <iomanip>
#include <iostream>
#include <random>
int main() {
    constexpr std::size_t n=192,tile=24;
    std::mt19937 rng(21);
    std::uniform_real_distribution<double> sample(-1.,1.);
    std::vector<double> a(n*n),b(n*n);
    for(auto& x:a) x=sample(rng);
    for(auto& x:b) x=sample(rng);
    const auto reference=graduate::naive(a,b,n,n,n);
    auto run=[&](bool blocked) {
        const auto start=std::chrono::steady_clock::now();
        const auto c=blocked ? graduate::blocked(a,b,n,n,n,tile)
                             : graduate::naive(a,b,n,n,n);
        const auto end=std::chrono::steady_clock::now();
        double checksum=0.;
        for(std::size_t i=0;i<c.size();++i) {
            if(std::abs(c[i]-reference[i])>1e-10*(1+std::abs(reference[i])))
                throw std::runtime_error("result mismatch");
            checksum+=c[i];
        }
        return std::pair{std::chrono::duration<double,std::milli>(end-start).count(),checksum};
    };
    (void)run(false); (void)run(true);
    std::cout<<"pair,first,naive_ms,blocked_ms,checksum\n"<<std::setprecision(12);
    for(int i=0;i<7;++i) {
        std::pair<double,double> old,newer;
        if(i%2==0) {old=run(false);newer=run(true);}
        else {newer=run(true);old=run(false);}
        std::cout<<i<<','<<(i%2==0?"naive":"blocked")<<','<<old.first<<','
                 <<newer.first<<','<<newer.second<<'\n';
    }
}

Readiness and continuation#

Advance when you can predict which loop strided access harms, explain the output ownership of a parallel version, and defend why a benchmark compares equivalent tasks. Further work includes hardware counters, vectorization reports, NUMA, GPU kernels, asynchronous pipelines, and target-specific libraries. The current compiler documentation was consulted for warnings/sanitizer controls, but sanitizer success is not proof of race freedom or all memory safety. [G-GCC]