Parallel Matrix Multiplication

Why Parallel Matrix Multiplication Matters

A company has an order matrix A whose rows are branches and columns are products. A second matrix B maps products to several price scenarios. Multiplying them produces every branch's total under every scenario. Each result combines one row of A with one column of B; many different output entries can be computed independently.

Parallel matrix multiplication appears in linear transformations, dense numerical methods, neural-network layers, batched simulations, and many data-processing pipelines. Its arithmetic is regular, but a fast implementation must also reuse data, expose enough independent work, and prevent workers from corrupting shared sums. A hundred available workers help little if one large task owns all the output, or if they repeatedly fetch the same inputs from slow memory.

This lesson implements classical dense multiplication with exclusive output-tile ownership. Each tile job computes all shared-index contributions to its own output cells. Jobs can run independently because they read the same immutable inputs and write disjoint output regions. We group the shared dimension into bands and carefully handle partial edge tiles.

The three programs are deterministic worker simulators, not threaded implementations. They print every scalar product, every updated partial sum, tile snapshots after each band, and a clearly defined arithmetic-only timing model. This makes correctness and scheduling visible before introducing actual concurrency or hardware-specific optimization.

Core Concepts and Terminology

Let $A\in\mathbb{Z}^{M\times K}$ and $B\in\mathbb{Z}^{K\times N}$. The inner dimensions must match. Their product $C=AB\in\mathbb{Z}^{M\times N}$ is defined by

$$C_{ij}=\sum_{t=0}^{K-1}A_{it}B_{tj},\qquad 0\le i\lt M,\quad 0\le j\lt N.$$

Indices start at zero. This is row times column, not entry-by-entry multiplication and not row times row. In general $AB\ne BA$, and one of these products may not even have compatible dimensions. If $K=0$, the empty sum is zero: C can be a nonempty zero matrix even though both flat input arrays are empty. If $M=0$ or $N=0$, C has no entries.

The implementation stores matrices in flat row-major arrays. A matrix with R rows and S columns stores row i, column j at offset $iS+j$. Therefore A uses $iK+t$, B uses $tN+j$, and C uses $iN+j$. Explicit M, K, N preserve zero-dimension shapes that cannot be recovered reliably from a bare empty nested list. A flat input must have exactly $MK$ or $KN$ entries respectively; there are no hidden row lengths or strides in this API.

With tile side T, output tile row r and column c covers rows $[rT,\min((r+1)T,M))$ and columns $[cT,\min((c+1)T,N))$. The half-open bounds include their start and exclude their end. Its owner processes all t values, grouped into bands $[q,\min(q+T,K))$. Clamping all three dimensions handles matrices whose sizes are not multiples of T.

A scalar MAC step in our model performs one multiplication and one addition: $C_{ij}\leftarrow C_{ij}+A_{it}B_{tj}$. We charge that pair one logical tick. This is an analysis unit, not a claim that integer multiplication/addition execute as one instruction or that a floating-point fused multiply-add is used. The code counts multiplications and additions separately, including addition to the initial zero.

A job executes its owned tile serially, one scalar MAC step at a time. Independent jobs can run on different workers; jobs within a wave start together, and the next wave waits for the slowest job. Work counts all MAC steps; tile span means the longest independent tile-job chain before scheduling barriers. The actual wave schedule can take longer. The diagram illustrates this distinction on the worked example: job costs are 12, 6, 6, and 3, so MAC work is 27, tile span is 12, and two-worker wave time is 18 ticks. None of these counts includes memory traffic or setup.

How Parallel Matrix Multiplication Works, Step by Step

Use immutable flat A and B, a fresh zero-filled C, tile side T, and P logical workers. Every tile owns a disjoint rectangular region of C for its entire computation, not just one shared-index band.

1. Validate before indexing. Require worker count $1\le P\le16$, tile side $1\le T\le8$, dimensions $0\le M,K,N\le8$, exact flat-array lengths, and integer values in $[-1000,1000]$. Validate in that order. The limits keep complete traces readable and all arithmetic exact in every language.
2. Allocate C once. Fill its $MN$ entries with zero. Do not zero it again when the next K band begins. For $K=0$ this initialization already produces the required zero result. With no output entries there are no tile jobs.
3. Enumerate output owners. Visit row-tile starts in increasing order and, within each, column-tile starts in increasing order. Clamp the row and column ends. Job j uses worker $j\bmod P$ and wave $\lfloor j/P\rfloor$. The tile rectangles partition all output coordinates exactly once.
4. Start a wave. All its jobs have the same modeled start tick. Within one job, loop over shared-index bands, then t within a band, then owned rows i, then owned columns j. Read $A_{it}$ and $B_{tj}$, compute their product, and add it to this owner's $C_{ij}$. Increment the job's MAC step count and record its operands, product, old/new accumulator, and logical completion tick.
5. Keep and record the partial tile. After each K band, save the owned tile's current values in row-major order. That snapshot is a partial result, not a reset or a new independent output. Continue through every t value in increasing order.
6. Finish the wave and continue. A job with d MAC steps finishes at start plus d. The wave ends at its longest job's finish tick; shorter jobs idle until then. The next wave starts at that end tick. Return C, all trace records, operation counts, job/wave counts, tile span, and schedule ticks.

The code executes one job at a time on the host while recording the arithmetic schedule they could follow concurrently. Its stdout groups updates by job, so a later printed job may contain earlier tick values. These are logical event records, not a globally time-sorted execution log. Within any one job, steps and shared indices increase normally.

For real threads, distinct tile owners need no lock for their C entries as long as A and B stay read-only and result memory does not alias either input. There is no required barrier between K bands of different tiles in this ownership design. The wave barrier is a scheduling choice; production task queues can let idle workers take the next tile instead. A cooperative GPU implementation that shares loaded A/B tiles among threads introduces its own memory-scope barriers and ownership rules, which are not implemented here.

flowchart TD
    A["Validate workers, tile, dimensions, sizes, values"] --> B["Allocate zero C; enumerate disjoint output tiles"]
    B --> C{"Any wave left?"}
    C -->|No| D["Return output, counts, and trace"]
    C -->|Yes| E["Start all wave jobs at the same logical tick"]
    E --> F["For each owner, keep its accumulator tile"]
    F --> G{"More shared-index bands?"}
    G -->|Yes| H["For each valid t, row, column: multiply and add"]
    H --> I["Record every update and the partial owned tile"]
    I --> G
    G -->|No| J{"More jobs in this wave?"}
    J -->|Yes| F
    J -->|No| K["Wave barrier: wait for maximum job finish"]
    K --> C

Worked Example with a Full Trace

Take $M=K=N=3$, $T=2$, and $P=2$, with rows

• A: [1,2,3], [4,5,6], [7,8,9].
• B: [9,8,7], [6,5,4], [3,2,1].

The expected C rows are [30,24,18], [84,69,54], and [138,114,90]. For example, $C_{12}=4\cdot7+5\cdot4+6\cdot1=54$. Each cell has three products. Starting from zero means three accumulation additions per cell, even though a hand-written sum of three already formed products can use only two additions.

There are four output jobs in row-major tile order. Job 0 owns rows $[0,2)$ and columns $[0,2)$; job 1 owns rows $[0,2)$ and columns $[2,3)$; job 2 owns rows $[2,3)$ and columns $[0,2)$; job 3 owns rows $[2,3)$ and columns $[2,3)$. Their cell counts are 4, 2, 2, and 1, so their MAC costs are 12, 6, 6, and 3. Every job uses K bands $[0,2)$ and $[2,3)$.

Job 0, worker 0, wave 0. In its first band, t=0 updates the owned tile from zero to rows [9,8] and [36,32]. Then t=1 adds rows [12,10] and [30,25], producing partial rows [21,18] and [66,57], flattened as [21,18,66,57]. In the final band, t=2 adds rows [9,6] and [18,12], yielding [30,24,84,69]. Logical ticks 1–4 perform t=0, 5–8 perform t=1, and 9–12 perform t=2. The job finishes at tick 12.

Job 1, worker 1, wave 0. For its right-edge tile, t=0 contributes [7,28]; t=1 adds [8,20], so the first band finishes with [15,48]. The final band adds [3,6] and gives [18,54]. It finishes at tick 6 and idles until wave 0 ends at tick 12. At that barrier C is [[30,24,18],[84,69,54],[0,0,0]].

Job 2, worker 0, wave 1. It starts at tick 12. For the lower-left tile, t=0 gives [63,56] and t=1 adds [48,40], yielding [111,96]. The final band adds [27,18] and gives [138,114]. Its six MAC steps complete at tick 18.

Job 3, worker 1, wave 1. Also starting at tick 12, its lone cell follows $0\to49\to81\to90$ from products 49, 32, and 9. It finishes at tick 15 and idles until wave 1 ends at tick 18. Now every cell is final. The diagram below lists all 27 products and all 27 resulting partial sums, with the same completion ticks as the programs.

Total arithmetic is 27 multiplications and 27 additions, or 27 MAC steps under the paired-step model. Two waves take $\max(12,6)+\max(6,3)=18$ modeled ticks, while independent tile span is 12. One worker takes 27 ticks; three workers take 15 under this wave policy; four or more take 12. Those differences reflect task grouping and idle time, not a change in C or in the total number of products.

Implementation in Python, C++, and Java

Each panel is a complete standalone program with identical inputs, loop order, edge bounds, trace records, and output. Python needs 3.9 or later. C++ uses C++17; compile with g++ -std=c++17 -O2 main.cpp -o matmul and run ./matmul. Save Java as Main.java, then use javac Main.java and java Main. Run Python with python3 main.py. The programs simulate workers sequentially; no threading or numeric libraries are required.

The shared API is multiply(a,b,m,k,n,tile,processors), with both input matrices supplied as flat row-major arrays. Explicit dimensions support zero rows, zero columns, and a zero shared dimension without ambiguity. C++ uses long long entries, Java uses long, and Python uses exact integers. Dimensions, tile sizes, and worker counts are bounded small integers. Validate errors in order: PROCESSORS, TILE, DIMENSION, SIZE, VALUE. Python additionally rejects booleans, floats, strings, and nulls where exact integers are required. The input containers must be non-null and flat; C++/Java primitive numeric arrays/vectors enforce their element types. Inputs are never modified, and C is freshly allocated.

With $K\le8$ and input magnitudes at most 1000, every product has magnitude at most $10^6$ and every partial sum at most $8\cdot10^6$. The programs therefore agree exactly without floating-point rounding or integer overflow. Every MAC increments both multiplication and addition counters, including $0+\text{first product}$. A printed tick is the logical completion time of that modeled MAC, not a timestamp or an observed CPU cycle.

Update records operands, product, old/new sum, coordinates, shared index, and completion tick. Block stores the K-band bounds, its updates, and the owned tile's snapshot. Job stores its exclusive output rectangle, worker, start/finish ticks, and blocks. Wave stores its jobs and barrier time; Result stores the output and totals. The printed tile_span is the largest independent serial tile cost, before this P-worker scheduling policy adds wave barriers. Retained snapshots are of owned tiles, not copies of the whole C after every arithmetic update.

For zero shared dimension, jobs still describe output ownership, but have no blocks or arithmetic steps. All output entries are zero, and model_ticks and tile_span are zero even if several zero-duration waves are listed. Allocation, validation, job bookkeeping, snapshots, input loads, output stores, and printing still take real time; they are deliberately excluded from the arithmetic-only model.

Verification: each language passed 6,567 input/configuration cases, including 6,555 accepted and 12 rejected. The suite covers every shape with M, K, N from 0 through 8 at six tile/worker combinations, exhaustive small integer/binary matrices, identity multiplication on both sides, extreme values, edge tiles, and seeded random cases. Independent dot-product and prefix-sum oracles checked 323,317 MAC updates, all intermediate products/partial sums and tile snapshots, exact per-cell ownership, logical event ticks, wave durations, counters, and unchanged inputs. Full formatted traces agree across languages; sample stdout matches byte for byte. Eleven additional Python checks reject invalid runtime types. C++ passed address/undefined-behavior sanitizer runs; Java compiled with all lint warnings enabled.

from dataclasses import dataclass

@dataclass
class Update:
    i: int
    j: int
    t: int
    left: int
    right: int
    product: int
    old: int
    new: int
    tick: int

@dataclass
class Block:
    lo: int
    hi: int
    updates: list[Update]
    tile_after: list[int]

@dataclass
class Job:
    number: int
    worker: int
    r0: int
    r1: int
    c0: int
    c1: int
    start: int
    finish: int
    blocks: list[Block]

@dataclass
class Wave:
    number: int
    start: int
    finish: int
    jobs: list[Job]

@dataclass
class Result:
    output: list[int]
    waves: list[Wave]
    multiplications: int
    additions: int
    ticks: int
    span: int
    jobs: int

def multiply(a: list[int], b: list[int], m: int, k: int, n: int,
             tile: int, processors: int) -> Result:
    if type(processors) is not int or not 1 <= processors <= 16:
        raise ValueError("PROCESSORS")
    if type(tile) is not int or not 1 <= tile <= 8:
        raise ValueError("TILE")
    if any(type(d) is not int or not 0 <= d <= 8 for d in (m, k, n)):
        raise ValueError("DIMENSION")
    if len(a) != m * k or len(b) != k * n:
        raise ValueError("SIZE")
    if any(type(x) is not int or not -1000 <= x <= 1000 for values in (a, b) for x in values):
        raise ValueError("VALUE")
    c = [0] * (m * n)
    rectangles = [(r0, min(r0 + tile, m), c0, min(c0 + tile, n))
                  for r0 in range(0, m, tile) for c0 in range(0, n, tile)]
    waves: list[Wave] = []
    ticks = count = span = 0
    for base in range(0, len(rectangles), processors):
        start = ticks
        jobs: list[Job] = []
        end = start
        for number in range(base, min(base + processors, len(rectangles))):
            r0, r1, c0, c1 = rectangles[number]
            step = 0
            blocks: list[Block] = []
            for lo in range(0, k, tile):
                hi = min(lo + tile, k)
                updates: list[Update] = []
                for t in range(lo, hi):
                    for i in range(r0, r1):
                        for j in range(c0, c1):
                            left, right = a[i * k + t], b[t * n + j]
                            product = left * right
                            index = i * n + j
                            old = c[index]
                            c[index] = old + product
                            step += 1
                            updates.append(Update(i, j, t, left, right, product,
                                                  old, c[index], start + step))
                snapshot = [c[i * n + j] for i in range(r0, r1)
                            for j in range(c0, c1)]
                blocks.append(Block(lo, hi, updates, snapshot))
            finish = start + step
            jobs.append(Job(number, number % processors, r0, r1, c0, c1,
                            start, finish, blocks))
            count += step
            span = max(span, step)
            end = max(end, finish)
        waves.append(Wave(base // processors, start, end, jobs))
        # The modeled wave barrier waits for its longest job.
        ticks = end
    return Result(c, waves, count, count, ticks, span, len(rectangles))

def array_text(values: list[int]) -> str:
    return "[" + ",".join(map(str, values)) + "]"

def matrix_text(values: list[int], rows: int, cols: int) -> str:
    return "[" + ",".join(array_text(values[i * cols:(i + 1) * cols])
                           for i in range(rows)) + "]"

def format_result(result: Result, m: int, n: int) -> str:
    lines: list[str] = []
    for wave in result.waves:
        lines.append(f"wave={wave.number} ticks=[{wave.start},{wave.finish})")
        for job in wave.jobs:
            lines.append(f"  job={job.number} worker={job.worker} "
                         f"rows=[{job.r0},{job.r1}) cols=[{job.c0},{job.c1}) "
                         f"ticks=[{job.start},{job.finish}) macs={job.finish - job.start}")
            for block in job.blocks:
                lines.append(f"    k=[{block.lo},{block.hi})")
                for u in block.updates:
                    lines.append(f"      C[{u.i},{u.j}] t={u.t} tick={u.tick} "
                                 f"a={u.left} b={u.right} product={u.product} old={u.old} new={u.new}")
                lines.append("    tile_after=" + array_text(block.tile_after))
    lines.append("C=" + matrix_text(result.output, m, n))
    lines.append(f"multiplications={result.multiplications} additions={result.additions} "
                 f"model_ticks={result.ticks} tile_span={result.span} jobs={result.jobs} "
                 f"waves={len(result.waves)}")
    return "\n".join(lines)

if __name__ == "__main__":
    m = k = n = 3
    a = [1, 2, 3, 4, 5, 6, 7, 8, 9]
    b = [9, 8, 7, 6, 5, 4, 3, 2, 1]
    tile, processors = 2, 2
    print(f"m={m} k={k} n={n} tile={tile} processors={processors}")
    print("A=" + matrix_text(a, m, k))
    print("B=" + matrix_text(b, k, n))
    print(format_result(multiply(a, b, m, k, n, tile, processors), m, n))

#include <algorithm>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
using Array = std::vector<long long>;
struct Update {
    int i, j, t;
    long long left, right, product, old, next;
    int tick;
};
struct Block { int lo, hi; std::vector<Update> updates; Array tile_after; };
struct Job {
    int number, worker, r0, r1, c0, c1, start, finish;
    std::vector<Block> blocks;
};
struct Wave { int number, start, finish; std::vector<Job> jobs; };
struct Result {
    Array output;
    std::vector<Wave> waves;
    int multiplications, additions, ticks, span, jobs;
};
struct Rectangle { int r0, r1, c0, c1; };
Result multiply(const Array& a, const Array& b, int m, int k, int n,
                int tile, int processors) {
    if (processors < 1 || processors > 16) throw std::invalid_argument("PROCESSORS");
    if (tile < 1 || tile > 8) throw std::invalid_argument("TILE");
    if (m < 0 || m > 8 || k < 0 || k > 8 || n < 0 || n > 8)
        throw std::invalid_argument("DIMENSION");
    if (a.size() != static_cast<std::size_t>(m*k) || b.size() != static_cast<std::size_t>(k*n))
        throw std::invalid_argument("SIZE");
    for (long long x : a) if (x < -1000 || x > 1000) throw std::invalid_argument("VALUE");
    for (long long x : b) if (x < -1000 || x > 1000) throw std::invalid_argument("VALUE");
    Array c(m*n,0); std::vector<Rectangle> rectangles;
    for (int r0 = 0; r0 < m; r0 += tile) {
        for (int c0 = 0; c0 < n; c0 += tile)
            rectangles.push_back(Rectangle{r0,std::min(r0+tile,m),c0,std::min(c0+tile,n)});
    }
    Result result{ {}, {}, 0, 0, 0, 0, static_cast<int>(rectangles.size())};
    for (int base = 0; base < static_cast<int>(rectangles.size()); base += processors) {
        int start = result.ticks, end = start;
        Wave wave{base/processors,start,start,{} };
        for (int number = base; number < std::min(base+processors,static_cast<int>(rectangles.size())); ++number) {
            const auto& rect = rectangles[number];
            Job job{number,number%processors,rect.r0,rect.r1,rect.c0,rect.c1,start,start,{} };
            int step = 0;
            for (int lo = 0; lo < k; lo += tile) {
                int hi = std::min(lo+tile,k); Block block{lo,hi,{}, {} };
                for (int t = lo; t < hi; ++t) {
                    for (int i = rect.r0; i < rect.r1; ++i) {
                        for (int j = rect.c0; j < rect.c1; ++j) {
                            long long left = a[i*k+t], right = b[t*n+j];
                            long long product = left*right, old = c[i*n+j];
                            c[i*n+j] = old+product; ++step;
                            block.updates.push_back(Update{i,j,t,left,right,product,old,c[i*n+j],start+step});
                        }
                    }
                }
                for (int i = rect.r0; i < rect.r1; ++i)
                    for (int j = rect.c0; j < rect.c1; ++j) block.tile_after.push_back(c[i*n+j]);
                job.blocks.push_back(block);
            }
            job.finish = start+step; wave.jobs.push_back(job);
            result.multiplications += step; result.additions += step;
            result.span = std::max(result.span,step); end = std::max(end,job.finish);
        }
        wave.finish = end; result.waves.push_back(wave);
        // The modeled wave barrier waits for its longest job.
        result.ticks = end;
    }
    result.output = c; return result;
}
std::string array_text(const Array& a) {
    std::ostringstream out; out << '[';
    for (std::size_t i = 0; i < a.size(); ++i) { if (i) out << ','; out << a[i]; }
    out << ']'; return out.str();
}
std::string matrix_text(const Array& a, int rows, int cols) {
    std::ostringstream out; out << '[';
    for (int i = 0; i < rows; ++i) {
        if (i) out << ',';
        out << '[';
        for (int j = 0; j < cols; ++j) { if (j) out << ','; out << a[i*cols+j]; }
        out << ']';
    }
    out << ']'; return out.str();
}
std::string format_result(const Result& r, int m, int n) {
    std::ostringstream out;
    for (const auto& wave : r.waves) {
        out << "wave=" << wave.number << " ticks=[" << wave.start << ',' << wave.finish << ")\n";
        for (const auto& job : wave.jobs) {
            out << "  job=" << job.number << " worker=" << job.worker
                << " rows=[" << job.r0 << ',' << job.r1 << ") cols=[" << job.c0 << ',' << job.c1
                << ") ticks=[" << job.start << ',' << job.finish << ") macs=" << job.finish-job.start << '\n';
            for (const auto& block : job.blocks) {
                out << "    k=[" << block.lo << ',' << block.hi << ")\n";
                for (const auto& u : block.updates) {
                    out << "      C[" << u.i << ',' << u.j << "] t=" << u.t << " tick=" << u.tick
                        << " a=" << u.left << " b=" << u.right << " product=" << u.product
                        << " old=" << u.old << " new=" << u.next << '\n';
                }
                out << "    tile_after=" << array_text(block.tile_after) << '\n';
            }
        }
    }
    out << "C=" << matrix_text(r.output,m,n) << '\n';
    out << "multiplications=" << r.multiplications << " additions=" << r.additions
        << " model_ticks=" << r.ticks << " tile_span=" << r.span << " jobs=" << r.jobs << " waves=" << r.waves.size();
    return out.str();
}
int main() {
    int m = 3, k = 3, n = 3, tile = 2, processors = 2;
    Array a{1,2,3,4,5,6,7,8,9}, b{9,8,7,6,5,4,3,2,1};
    std::cout << "m=" << m << " k=" << k << " n=" << n << " tile=" << tile << " processors=" << processors << '\n';
    std::cout << "A=" << matrix_text(a,m,k) << '\n' << "B=" << matrix_text(b,k,n) << '\n';
    std::cout << format_result(multiply(a,b,m,k,n,tile,processors),m,n) << '\n';
    return 0;
}

import java.util.ArrayList;
import java.util.List;

public class Main {
    static final class Update {
        int i, j, t, tick;
        long left, right, product, old, next;
        Update(int i, int j, int t, long left, long right, long product,
               long old, long next, int tick) {
            this.i=i; this.j=j; this.t=t; this.left=left; this.right=right;
            this.product=product; this.old=old; this.next=next; this.tick=tick;
        }
    }
    static final class Block {
        int lo, hi;
        List<Update> updates = new ArrayList<>();
        List<Long> tileAfter = new ArrayList<>();
        Block(int lo, int hi) { this.lo=lo; this.hi=hi; }
    }
    static final class Job {
        int number, worker, r0, r1, c0, c1, start, finish;
        List<Block> blocks = new ArrayList<>();
        Job(int number, int worker, int r0, int r1, int c0, int c1, int start) {
            this.number=number; this.worker=worker; this.r0=r0; this.r1=r1;
            this.c0=c0; this.c1=c1; this.start=start;
        }
    }
    static final class Wave {
        int number, start, finish;
        List<Job> jobs = new ArrayList<>();
        Wave(int number, int start) { this.number=number; this.start=start; }
    }
    static final class Result {
        long[] output;
        List<Wave> waves = new ArrayList<>();
        int multiplications, additions, ticks, span, jobs;
    }
    static Result multiply(long[] a, long[] b, int m, int k, int n,
                           int tile, int processors) {
        if (processors<1 || processors>16) throw new IllegalArgumentException("PROCESSORS");
        if (tile<1 || tile>8) throw new IllegalArgumentException("TILE");
        if (m<0 || m>8 || k<0 || k>8 || n<0 || n>8)
            throw new IllegalArgumentException("DIMENSION");
        if (a.length!=m*k || b.length!=k*n) throw new IllegalArgumentException("SIZE");
        for (long x : a) if (x < -1000 || x > 1000) throw new IllegalArgumentException("VALUE");
        for (long x : b) if (x < -1000 || x > 1000) throw new IllegalArgumentException("VALUE");
        long[] c = new long[m*n]; List<int[]> rectangles = new ArrayList<>();
        for (int r0=0; r0<m; r0+=tile)
            for (int c0=0; c0<n; c0+=tile)
                rectangles.add(new int[]{r0,Math.min(r0+tile,m),c0,Math.min(c0+tile,n)});
        Result result = new Result(); result.jobs=rectangles.size();
        for (int base=0; base<rectangles.size(); base+=processors) {
            int start=result.ticks, end=start; Wave wave=new Wave(base/processors,start);
            for (int number=base; number<Math.min(base+processors,rectangles.size()); ++number) {
                int[] rect=rectangles.get(number);
                Job job=new Job(number,number%processors,rect[0],rect[1],rect[2],rect[3],start);
                int step=0;
                for (int lo=0; lo<k; lo+=tile) {
                    int hi=Math.min(lo+tile,k); Block block=new Block(lo,hi);
                    for (int t=lo; t<hi; ++t) {
                        for (int i=job.r0; i<job.r1; ++i) {
                            for (int j=job.c0; j<job.c1; ++j) {
                                long left=a[i*k+t], right=b[t*n+j];
                                long product=left*right, old=c[i*n+j];
                                c[i*n+j]=old+product; ++step;
                                block.updates.add(new Update(i,j,t,left,right,product,old,c[i*n+j],start+step));
                            }
                        }
                    }
                    for (int i=job.r0; i<job.r1; ++i)
                        for (int j=job.c0; j<job.c1; ++j) block.tileAfter.add(c[i*n+j]);
                    job.blocks.add(block);
                }
                job.finish=start+step; wave.jobs.add(job);
                result.multiplications+=step; result.additions+=step;
                result.span=Math.max(result.span,step); end=Math.max(end,job.finish);
            }
            wave.finish=end; result.waves.add(wave);
            // The modeled wave barrier waits for its longest job.
            result.ticks=end;
        }
        result.output=c; return result;
    }
    static String arrayText(List<Long> values) {
        StringBuilder out=new StringBuilder("[");
        for (int i=0; i<values.size(); ++i) { if (i!=0) out.append(','); out.append(values.get(i)); }
        return out.append(']').toString();
    }
    static String matrixText(long[] a, int rows, int cols) {
        StringBuilder out=new StringBuilder("[");
        for (int i=0; i<rows; ++i) {
            if (i!=0) out.append(','); out.append('[');
            for (int j=0; j<cols; ++j) { if (j!=0) out.append(','); out.append(a[i*cols+j]); }
            out.append(']');
        }
        return out.append(']').toString();
    }
    static String formatResult(Result r, int m, int n) {
        StringBuilder out=new StringBuilder();
        for (Wave wave : r.waves) {
            out.append("wave=").append(wave.number).append(" ticks=[").append(wave.start)
               .append(',').append(wave.finish).append(")\n");
            for (Job job : wave.jobs) {
                out.append("  job=").append(job.number).append(" worker=").append(job.worker)
                   .append(" rows=[").append(job.r0).append(',').append(job.r1)
                   .append(") cols=[").append(job.c0).append(',').append(job.c1)
                   .append(") ticks=[").append(job.start).append(',').append(job.finish)
                   .append(") macs=").append(job.finish-job.start).append('\n');
                for (Block block : job.blocks) {
                    out.append("    k=[").append(block.lo).append(',').append(block.hi).append(")\n");
                    for (Update u : block.updates) {
                        out.append("      C[").append(u.i).append(',').append(u.j).append("] t=").append(u.t)
                           .append(" tick=").append(u.tick).append(" a=").append(u.left).append(" b=").append(u.right)
                           .append(" product=").append(u.product).append(" old=").append(u.old)
                           .append(" new=").append(u.next).append('\n');
                    }
                    out.append("    tile_after=").append(arrayText(block.tileAfter)).append('\n');
                }
            }
        }
        return out.append("C=").append(matrixText(r.output,m,n)).append('\n')
                  .append("multiplications=").append(r.multiplications).append(" additions=").append(r.additions)
                  .append(" model_ticks=").append(r.ticks).append(" tile_span=").append(r.span)
                  .append(" jobs=").append(r.jobs).append(" waves=").append(r.waves.size()).toString();
    }
    public static void main(String[] args) {
        int m=3, k=3, n=3, tile=2, processors=2;
        long[] a={1,2,3,4,5,6,7,8,9}, b={9,8,7,6,5,4,3,2,1};
        System.out.println("m="+m+" k="+k+" n="+n+" tile="+tile+" processors="+processors);
        System.out.println("A="+matrixText(a,m,k)); System.out.println("B="+matrixText(b,k,n));
        System.out.println(formatResult(multiply(a,b,m,k,n,tile,processors),m,n));
    }
}

Identical expected output (each item is one printed line; events are grouped by owner rather than globally sorted by tick):

• m=3 k=3 n=3 tile=2 processors=2
• A=[[1,2,3],[4,5,6],[7,8,9]]
• B=[[9,8,7],[6,5,4],[3,2,1]]
• wave=0 ticks=[0,12)
• job=0 worker=0 rows=[0,2) cols=[0,2) ticks=[0,12) macs=12
• k=[0,2)
• C[0,0] t=0 tick=1 a=1 b=9 product=9 old=0 new=9
• C[0,1] t=0 tick=2 a=1 b=8 product=8 old=0 new=8
• C[1,0] t=0 tick=3 a=4 b=9 product=36 old=0 new=36
• C[1,1] t=0 tick=4 a=4 b=8 product=32 old=0 new=32
• C[0,0] t=1 tick=5 a=2 b=6 product=12 old=9 new=21
• C[0,1] t=1 tick=6 a=2 b=5 product=10 old=8 new=18
• C[1,0] t=1 tick=7 a=5 b=6 product=30 old=36 new=66
• C[1,1] t=1 tick=8 a=5 b=5 product=25 old=32 new=57
• tile_after=[21,18,66,57]
• k=[2,3)
• C[0,0] t=2 tick=9 a=3 b=3 product=9 old=21 new=30
• C[0,1] t=2 tick=10 a=3 b=2 product=6 old=18 new=24
• C[1,0] t=2 tick=11 a=6 b=3 product=18 old=66 new=84
• C[1,1] t=2 tick=12 a=6 b=2 product=12 old=57 new=69
• tile_after=[30,24,84,69]
• job=1 worker=1 rows=[0,2) cols=[2,3) ticks=[0,6) macs=6
• k=[0,2)
• C[0,2] t=0 tick=1 a=1 b=7 product=7 old=0 new=7
• C[1,2] t=0 tick=2 a=4 b=7 product=28 old=0 new=28
• C[0,2] t=1 tick=3 a=2 b=4 product=8 old=7 new=15
• C[1,2] t=1 tick=4 a=5 b=4 product=20 old=28 new=48
• tile_after=[15,48]
• k=[2,3)
• C[0,2] t=2 tick=5 a=3 b=1 product=3 old=15 new=18
• C[1,2] t=2 tick=6 a=6 b=1 product=6 old=48 new=54
• tile_after=[18,54]
• wave=1 ticks=[12,18)
• job=2 worker=0 rows=[2,3) cols=[0,2) ticks=[12,18) macs=6
• k=[0,2)
• C[2,0] t=0 tick=13 a=7 b=9 product=63 old=0 new=63
• C[2,1] t=0 tick=14 a=7 b=8 product=56 old=0 new=56
• C[2,0] t=1 tick=15 a=8 b=6 product=48 old=63 new=111
• C[2,1] t=1 tick=16 a=8 b=5 product=40 old=56 new=96
• tile_after=[111,96]
• k=[2,3)
• C[2,0] t=2 tick=17 a=9 b=3 product=27 old=111 new=138
• C[2,1] t=2 tick=18 a=9 b=2 product=18 old=96 new=114
• tile_after=[138,114]
• job=3 worker=1 rows=[2,3) cols=[2,3) ticks=[12,15) macs=3
• k=[0,2)
• C[2,2] t=0 tick=13 a=7 b=7 product=49 old=0 new=49
• C[2,2] t=1 tick=14 a=8 b=4 product=32 old=49 new=81
• tile_after=[81]
• k=[2,3)
• C[2,2] t=2 tick=15 a=9 b=1 product=9 old=81 new=90
• tile_after=[90]
• C=[[30,24,18],[84,69,54],[138,114,90]]
• multiplications=27 additions=27 model_ticks=18 tile_span=12 jobs=4 waves=2

Complexity, Correctness, and Invariants

Assume constant-size exact values and constant-cost arithmetic. The asymptotic discussion describes the same algorithm without the small teaching limits. Use M rows, K shared coordinates, N output columns, tile side T, and P workers; keep them distinct instead of silently assuming square matrices.

Arithmetic work. Each of the $MN$ output entries processes K terms, giving

$$W_{\mathrm{MAC}}=MNK,\qquad W_{\times}=MNK,\qquad W_{+}=MNK.$$

This implementation performs $2MNK$ counted arithmetic operations because it adds each product to an initially zero accumulator. The common alternative count $MNK$ multiplications plus $MN(K-1)$ additions applies to nonempty dot products initialized from their first product; it is not this program's counter convention. All supported values, including zeros, trigger the same loop structure. Best, average, and worst arithmetic work are identical for fixed shapes. For two square matrices with $M=K=N=s\ge1$, classical work is $\Theta(s^3)$, not $\Theta(s^2)$.

Tile jobs and independent span. Let

$$R=\left\lceil\frac{M}{T}\right\rceil,\qquad S=\left\lceil\frac{N}{T}\right\rceil,\qquad Q=RS.$$

Each tile job q owns $r_q$ rows and $c_q$ columns and has cost $d_q=Kr_qc_q$ MAC ticks. These jobs are independent but execute their own cells serially. Consequently, for nonempty output, the tile-job DAG has span

$$D_{\mathrm{tile}}=\max_qd_q=K\min(T,M)\min(T,N).$$

With no output entries define this span as zero. If $K=0$, it is also zero. This is the span of the independent arithmetic jobs before the scheduler's wave barriers; it is not the span of the already serialized host program, nor a claim about a GPU block whose threads cooperate inside a tile. With $T=1$, output cells are separate jobs and span is K. A single tile containing all output entries has serial job cost $MNK$, so this teaching job grouping can expose no arithmetic concurrency in that case.

Exact wave schedule. Consecutive groups of at most P tile jobs form waves. If wave w contains job set $G_w$, its start follows the preceding wave and its duration is $\max_{q\in G_w}d_q$. Thus

$$T_{\mathrm{wave}}=\sum_w\max_{q\in G_w}d_q,\qquad \text{wave count}=\left\lceil\frac{Q}{P}\right\rceil.$$

For positive work, $T_{\mathrm{wave}}\ge\max(W_{\mathrm{MAC}}/P,D_{\mathrm{tile}})$ and $T_{\mathrm{wave}}\le W_{\mathrm{MAC}}$. One worker takes exactly $W_{\mathrm{MAC}}$ ticks; P at least Q takes exactly $D_{\mathrm{tile}}$. If both output dimensions are multiples of T, all tiles have cost $KT^2$, so $T_{\mathrm{wave}}=KT^2\lceil Q/P\rceil$. Partial edges produce unequal costs, so multiplying the wave count by one assumed uniform tile cost is generally wrong. A suitable dynamic work-conserving list schedule on independent jobs can obtain $O(W_{\mathrm{MAC}}/P+D_{\mathrm{tile}})$ arithmetic time; do not assume our fixed barrier groups always attain that bound tightly.

On the worked input, $W_{\mathrm{MAC}}=27$, $D_{\mathrm{tile}}=12$, and $T_{\mathrm{wave}}=18$ for two workers. Arithmetic-only speedup against the same one-worker job schedule is $27/18=1.5$ and efficiency is $1.5/2=0.75$. These are modeled quantities. Actual initialization, input/output, memory bandwidth, cache reuse, parallel startup, and barriers can alter measured speedup substantially.

A different decomposition gives a different span. Assigning one worker per cell but accumulating its K products sequentially gives span $\Theta(K)$. Computing a cell's products in parallel and reducing them with a balanced tree can lower arithmetic span to $\Theta(1+\log K)$ for positive K, while retaining $\Theta(MNK)$ work for positive dimensions. A straightforward all-products temporary array uses $O(MNK)$ space. Private partial tiles plus reduction offer intermediate tradeoffs; a shared unsynchronized accumulator is not a valid replacement. Changing the reduction tree also changes floating-point grouping.

Space and real simulator cost. C occupies $\Theta(MN)$ values. Input storage is $\Theta(MK+KN)$ and is supplied by the caller. Without the saved history, output ownership needs Q task descriptors and small loop state, with no extra arrays of A/B tiles in this particular code. With full history, there are $MNK$ update records and $MN\lceil K/T\rceil$ saved owned-tile entries. Since T is at least one, the history plus output occupies $O(MN(K+1)+1)$ records; no full-matrix snapshot is saved per update. Job/wave descriptors are bounded by $O(MN)$ when output exists. Including validation and output initialization, valid-input sequential host work is

$$\Theta(MK+KN+MN+MNK+1)$$

in a word/record model, before accounting for text lengths. The host does not run in $T_{\mathrm{wave}}$ time. For $K=0$, creating a nonempty zero output still costs $\Theta(MN)$ real work even though the arithmetic model reports zero ticks. For zero output size, validating a nonempty input can still take time.

Correctness invariant. For each owned cell, immediately before processing shared index t, its accumulator equals

$$C_{ij}^{(t)}=\sum_{u=0}^{t-1}A_{iu}B_{uj}.$$

Initialization establishes the empty prefix sum at t=0. One update adds exactly $A_{it}B_{tj}$, establishing the formula for t+1. Bands only group consecutive t values; they do not alter this induction or reset its accumulator. After all K values, the cell equals the required dot product. The tile grid covers each output coordinate exactly once, different owners write disjoint cells, and all inputs stay immutable. Therefore independently executing the tile jobs preserves every cell invariant and yields the correct matrix. Zero shared dimension follows directly from the empty-sum convention.

Variants, Pitfalls, and Optimizations

Do not share an unprotected output reduction. Parallelizing t while multiple workers execute C[i,j] += product can lose updates. For $C_{00}$, the first products are 9 and 12. If two workers both read zero, one can write 9 and the other write 12 from its stale read; the final product 9 then produces 21 instead of 30. A normal read-modify-write is not atomic. Keep one output owner, or compute private partial sums/tiles and combine them through a defined reduction. Atomics may serialize or introduce numerical nondeterminism; they are a design choice, not a free substitute for ownership.

Initialize once and clamp every dimension. Resetting C at each K band discards previous terms. Forgetting the row, column, or shared-index edge guard can access past an array or multiply unintended data. For T=2 and dimension 3, final ranges have length 1. The loops use actual bounds and do no padded arithmetic. This differs from kernels that deliberately compute masked or zero-padded full tiles; do not reuse their operation counts for this code.

Use the correct strides and dimensions. A uses stride K, B and C use stride N. Accidentally using M as a stride may pass square tests and fail rectangular matrices. A transpose changes both shape and indexing; it is not achieved by replacing one stride. The API's explicit dimensions reject missing/excess entries, including shapes with zero shared dimension. If introducing leading dimensions or submatrix views, add bounds and stride contracts rather than assuming dense contiguous rows.

Blocking is about reuse as well as scheduling. Our K bands organize loop traversal, but the code does not explicitly pack input tiles, load shared memory, use vector instructions, or measure cache traffic. For a full T-by-T-by-T tile product, $T^3$ MACs can reuse two $T^2$ input tiles; this is the locality opportunity. Effective reuse still depends on cache/register capacity, loop order, layout, and hardware. One commonly used CPU cache-sizing heuristic considers space for A, B, and C tiles, about $3T^2$ elements, but real kernels often use different rectangular blocking levels and packing buffers. A larger tile can reduce task count and reuse input while also reducing exposed parallelism in our serial-tile model. These effects must be measured separately.

Cooperative GPU kernels need their own barriers. A GPU block commonly lets threads load an A/B tile, synchronize before consuming shared data, accumulate private register results, and synchronize before overwriting that shared tile for the next K band. Those synchronization points protect the shared buffers; they are different from our logical wave barrier and from the absence of inter-tile data dependence. Respect the scope and participation requirements of the platform's barriers, including boundary threads. This lesson's sequential loops are not a ready-to-launch GPU kernel.

Use a deliberate numeric policy. Exact bounded integers avoid overflow here. With arbitrary C++ signed integers, overflow is undefined; Java fixed-width integer operations wrap. Promote operands before multiplication and validate bounds or use checked/wider arithmetic as required. Floating-point addition is not associative: a tree or split-K reduction can differ from a serial dot product. This particular owner loop preserves increasing-t accumulation order for each cell across T and P, but real libraries may reorder operations or use fused/mixed-precision instructions. Define tolerances, reproducibility requirements, and accumulator precision when comparing floating-point results.

Choose the larger algorithm to fit the data. Tuned GEMM typically computes $C\leftarrow\alpha\,\operatorname{op}(A)\operatorname{op}(B)+\beta C$ and supports layouts/transposition beyond this lesson's $C=AB$. Sparse multiplication changes representation, work, and load balance; do not apply dense MNK counts to it. Batched small products need overhead amortization. Distributed 2D block algorithms must also route/broadcast panels and manage communication. Strassen-style methods reduce square multiplication's arithmetic exponent but introduce extra additions, temporary storage, cutoffs, and numerical tradeoffs; they do not follow automatically from parallelizing these loops. For application code, a suitable tuned linear-algebra implementation is usually the baseline worth measuring against.

Practice Problems and Real-World Use

Start with a direct sequential row-by-column reference, then inspect the recorded prefix sum after every shared index. On rectangular cases, verify the exact row-major offsets for all three matrices. Count output owners and confirm their rectangles neither overlap nor omit a cell. Changing tile side or worker count must preserve the exact integer output and the multiset of coordinate/product contributions, while scheduling and trace grouping can change.

Useful properties include $AI=A$, $IA=A$ for compatible identity sizes, multiplication by a zero matrix, and $(AB)^{\mathsf T}=B^{\mathsf T}A^{\mathsf T}$. Over exact arithmetic, distributivity also provides tests when shapes match. Keep generated test values within the API's bounds. Include $K=0$ with nonempty output, $M=0$, $N=0$, negative values, cancellation, singleton dimensions, skinny/tall products, and sizes just below/at/above a tile boundary. Test invalid lengths and dimensions before trusting randomized arithmetic tests.

For a concurrent extension, retain one output owner and replace the deterministic host loop with actual independent tasks. First compare outputs and ownership under repeated runs. Then add a task queue, tile reuse or packing, and an appropriate cutoff for tiny jobs. A split-K experiment should use private partial tiles plus a reduction, making its extra memory and different reduction order visible. Never interpret the simulator's tick field as a measured race test or a hardware benchmark.

For performance reports, record matrix shapes, data types, layouts, tile/block sizes, workers, arithmetic policy, reference implementation, and whether allocation/transfers are included. Compare arithmetic throughput alongside elapsed time, memory movement, and load imbalance. A larger tile's cache benefit cannot be inferred from this arithmetic-only scheduler, and a smaller tile's modeled improvement can be outweighed by task overhead on real systems.

Further reading: Netlib's GEMM specification defines the broader general matrix-multiplication operation and shape parameters. NVIDIA's matrix multiplication background discusses tiling and hardware performance considerations. CUTLASS's tiling explanation illustrates hierarchical reuse and accumulation. Their hardware/library discussions provide context; the exact counts, ownership, and wave timing above describe this lesson's serial-tile simulator.