24. Closing the Gap: A Matrix Product That Walks Memory in a Better Order

What you will understand: how a change to the order of three loops, with no change to the arithmetic, can make a program several times faster, and how to add that change to a compiler without breaking anything that depends on the old behavior. Chapter 23 measured Mountain Goat's matrix product, found it more than twice as slow as a simple hand-written C++ loop at N = 512, and listed "isolating why" as unfinished. This chapter finishes the job: it adds a second loop order to the lowering, measures it, looks at the machine code, and shows that the two orders give identical answers to the last bit.

What you need to know first: Chapter 23 (the benchmark, the noise, -O and --passes), Chapter 21 (how mg.matmul is lowered), and Chapter 7 (loop tiling). The new idea (loop order and memory access) is explained from scratch.

Compile and run

cd docs/part24/code
./build.sh                                 # builds ./build/mg-opt (Chapter 24's: adds the matmul-order option)
./mgc run examples/matmul_odd.mg --matmul-order ikj       # one program with the new loop order
./run_bench.sh > bench_out_1.txt           # the benchmark, pass 1 (about a minute)
./run_bench.sh > bench_out_2.txt           # pass 2
python3 summarize.py bench_out_1.txt bench_out_2.txt > summary.txt
./inspect_asm.sh > inspect_asm_out.txt     # vector instructions in the machine code
./bits_check.sh > bits_check_out.txt       # are the two orders identical to the last bit?
./machine_info.sh > machine_info_out.txt
./mgc_mutation.py > mgc_mutation_out.txt   # breaks mgc on purpose
./lowering_mutation.sh > lowering_mutation_out.txt   # breaks the lowering on purpose (a rebuild each; several minutes)
cd ../../part15/code && ./run_lit.sh       # the whole test suite

Timings will differ on your machine.

Why this page shows FAIL lines, and why that is good

The two mutation outputs near the end come from negative controls: deliberately broken code run against the tests. There a "caught by" line or a FAIL is the expected, wanted result. The baseline lines must say "all tests pass", and they do.

Primer: why loop order matters

The matrix product of A (m by k) and B (k by n) is C[i][j] = sum over k of A[i][k] * B[k][j]. The textbook way to compute it is three nested loops. Which loop is outermost, which is in the middle and which is innermost can be chosen freely, because the arithmetic is the same. The two orders in this chapter are named by the loop variables from outermost to innermost:

• ijk (what Chapters 21 to 23 generate): for each result element C[i][j], run a loop over k that adds up the products. The innermost loop reads A[i][k] (walking along a row of A: consecutive addresses) and B[k][j] (walking down a column of B).
• ikj: for each row i of A and each k, take the single number A[i][k] and add A[i][k] * B[k][j] into C[i][j] for every j. The innermost loop now reads B[k][j] and writes C[i][j] while j increases, and both of those walk along a row: consecutive addresses.

Why that matters. Memory is stored one row after another ("row-major"), so walking along a row visits neighboring addresses, and walking down a column jumps by a whole row's length each step (for a 512-wide matrix, 4096 bytes). A processor reads memory in chunks called cache lines (typically 64 bytes, eight numbers); a walk along a row uses all eight numbers of each chunk it fetches, while a walk down a column uses one number per chunk and wastes the rest. There is a second benefit: a loop that touches consecutive addresses is the kind of loop a compiler can turn into vector instructions (SIMD: one instruction operating on several numbers), which a loop jumping through memory is not.

What does not change. For every result element C[i][j], both orders add the products in the same order (k = 0, 1, 2, …). Floating-point addition is not associative (adding the same numbers in a different order can change the last bits), but because the order per element is unchanged, the two loop orders give bit-for-bit identical results. That is an argument; this chapter also tests it (below).

The change

mg.matmul already lowers to a zero-fill loop nest plus a three-deep accumulate nest (Chapter 21). The new code builds the accumulate nest with its bounds and induction variables in either order, selected by a pass option, --convert-mg-to-affine=matmul-order=ikj (the default stays ijk, so every earlier chapter's output is unchanged). The lowering pattern now remembers which order it was asked for:

struct MatmulOpLowering : public OpConversionPattern<mg::MatmulOp> {
  /// Chapter 24: `ikj` runs the three loops in the order i, k, j instead of i, j, k. Same arithmetic, same
  /// order of additions into each result element (k increases for every i, j), so the results are bit-for-bit
  /// identical; only the memory access pattern differs.
  bool ikj;
  MatmulOpLowering(TypeConverter &converter, MLIRContext *ctx, bool ikj)
      : OpConversionPattern(converter, ctx), ikj(ikj) {}
  LogicalResult matchAndRewrite(mg::MatmulOp op, OpAdaptor adaptor,
                                ConversionPatternRewriter &rewriter) const override {
    auto loc = op.getLoc();

The accumulate nest, with the order choice in the loop bounds and in which induction variable plays which role:

    Value n = bound(outType.getDimSize(1), rhs, 1, Value());
    int64_t kStatic = !lt.isDynamicDim(1) ? lt.getDimSize(1) : rt.getDimSize(0);
    Value k = bound(kStatic, lhs, 1, Value());
    Value c0 = rewriter.create<arith::ConstantIndexOp>(loc, 0);
    // Loop bounds in nest order: (m, n, k) for ijk, (m, k, n) for ikj. ivs follow the same order.
    llvm::SmallVector<Value, 3> bounds = ikj ? llvm::SmallVector<Value, 3>{m, k, n} : llvm::SmallVector<Value, 3>{m, n, k};
    affine::buildAffineLoopNest(rewriter, loc, ValueRange{c0, c0, c0}, bounds, {1, 1, 1},
        [&](OpBuilder &b, Location l, ValueRange ivs) {
          Value i = ivs[0], j = ikj ? ivs[2] : ivs[1], kk = ikj ? ivs[1] : ivs[2];
          Value acc = b.create<affine::AffineLoadOp>(l, alloc, ValueRange{i, j});
          Value x = b.create<affine::AffineLoadOp>(l, lhs, ValueRange{i, kk});
          Value y = b.create<affine::AffineLoadOp>(l, rhs, ValueRange{kk, j});
          Value prod = b.create<arith::MulFOp>(l, x, y);
          b.create<affine::AffineStoreOp>(l, b.create<arith::AddFOp>(l, acc, prod), alloc, ValueRange{i, j});
        });

An option on a pass is declared as a member of the pass class. Because a pass with options must be copyable (MLIR clones passes), the class also gets an explicit copy constructor, and the pass rejects any value other than ijk or ikj with an error instead of silently falling back:

  MLIR_DEFINE_EXPLICIT_INTERNAL_INLINE_TYPE_ID(ConvertMgToAffinePass)

  ConvertMgToAffinePass() = default;
  ConvertMgToAffinePass(const ConvertMgToAffinePass &other) : PassWrapper(other) {}   // needed because of the Option below

  Option<std::string> matmulOrder{*this, "matmul-order",
                                  llvm::cl::desc("loop order for mg.matmul: ijk (the default, as in Chapters 21 to 23) or ikj"),
                                  llvm::cl::init("ijk")};

  llvm::StringRef getArgument() const final { return "convert-mg-to-affine"; }
  llvm::StringRef getDescription() const final {
    return "Lower Mountain Goat's own mg dialect to affine/memref/arith, "
           "hand-bufferizing each op's own result as it is lowered.";
  }

  void getDependentDialects(DialectRegistry &registry) const override {
    registry.insert<affine::AffineDialect, memref::MemRefDialect,
                    arith::ArithDialect, cf::ControlFlowDialect,
                    tensor::TensorDialect>();
  }

  void runOnOperation() override {
    if (matmulOrder != "ijk" && matmulOrder != "ikj") {
      getOperation().emitError("convert-mg-to-affine: matmul-order must be ijk or ikj, got '") << matmulOrder << "'";
      return signalPassFailure();
    }
    MgToAffineTypeConverter typeConverter;

mgc gets a flag that passes it through: --matmul-order ikj. The complete change (this chapter's lowering and driver, against Chapter 22's and 23's):

The full diff for this chapter (click to expand)

--- ../../part22/code/LowerToAffine.cpp 2026-10-02 11:17:33.277001200 +0000
+++ LowerToAffine.cpp   2026-10-02 18:03:38.091631141 +0000
@@ -281,7 +281,12 @@
 /// mg.matmul: result[i,j] = sum over k of lhs[i,k] * rhs[k,j]. A zero-fill nest, then a 3-deep
 /// accumulate nest. For `?` dimensions the inner extents must agree at run time: one cf.assert.
 struct MatmulOpLowering : public OpConversionPattern<mg::MatmulOp> {
-  using OpConversionPattern::OpConversionPattern;
+  /// Chapter 24: `ikj` runs the three loops in the order i, k, j instead of i, j, k. Same arithmetic, same
+  /// order of additions into each result element (k increases for every i, j), so the results are bit-for-bit
+  /// identical; only the memory access pattern differs.
+  bool ikj;
+  MatmulOpLowering(TypeConverter &converter, MLIRContext *ctx, bool ikj)
+      : OpConversionPattern(converter, ctx), ikj(ikj) {}
   LogicalResult matchAndRewrite(mg::MatmulOp op, OpAdaptor adaptor,
                                 ConversionPatternRewriter &rewriter) const override {
     auto loc = op.getLoc();
@@ -315,9 +320,11 @@
     int64_t kStatic = !lt.isDynamicDim(1) ? lt.getDimSize(1) : rt.getDimSize(0);
     Value k = bound(kStatic, lhs, 1, Value());
     Value c0 = rewriter.create<arith::ConstantIndexOp>(loc, 0);
-    affine::buildAffineLoopNest(rewriter, loc, ValueRange{c0, c0, c0}, ValueRange{m, n, k}, {1, 1, 1},
+    // Loop bounds in nest order: (m, n, k) for ijk, (m, k, n) for ikj. ivs follow the same order.
+    llvm::SmallVector<Value, 3> bounds = ikj ? llvm::SmallVector<Value, 3>{m, k, n} : llvm::SmallVector<Value, 3>{m, n, k};
+    affine::buildAffineLoopNest(rewriter, loc, ValueRange{c0, c0, c0}, bounds, {1, 1, 1},
         [&](OpBuilder &b, Location l, ValueRange ivs) {
-          Value i = ivs[0], j = ivs[1], kk = ivs[2];
+          Value i = ivs[0], j = ikj ? ivs[2] : ivs[1], kk = ikj ? ivs[1] : ivs[2];
           Value acc = b.create<affine::AffineLoadOp>(l, alloc, ValueRange{i, j});
           Value x = b.create<affine::AffineLoadOp>(l, lhs, ValueRange{i, kk});
           Value y = b.create<affine::AffineLoadOp>(l, rhs, ValueRange{kk, j});
@@ -440,6 +447,13 @@
     : public PassWrapper<ConvertMgToAffinePass, OperationPass<ModuleOp>> {
   MLIR_DEFINE_EXPLICIT_INTERNAL_INLINE_TYPE_ID(ConvertMgToAffinePass)

+  ConvertMgToAffinePass() = default;
+  ConvertMgToAffinePass(const ConvertMgToAffinePass &other) : PassWrapper(other) {}   // needed because of the Option below
+
+  Option<std::string> matmulOrder{*this, "matmul-order",
+                                  llvm::cl::desc("loop order for mg.matmul: ijk (the default, as in Chapters 21 to 23) or ikj"),
+                                  llvm::cl::init("ijk")};
+
   llvm::StringRef getArgument() const final { return "convert-mg-to-affine"; }
   llvm::StringRef getDescription() const final {
     return "Lower Mountain Goat's own mg dialect to affine/memref/arith, "
@@ -453,6 +467,10 @@
   }

   void runOnOperation() override {
+    if (matmulOrder != "ijk" && matmulOrder != "ikj") {
+      getOperation().emitError("convert-mg-to-affine: matmul-order must be ijk or ikj, got '") << matmulOrder << "'";
+      return signalPassFailure();
+    }
     MgToAffineTypeConverter typeConverter;
     ConversionTarget target(getContext());
     target.addLegalDialect<affine::AffineDialect, memref::MemRefDialect,
@@ -486,8 +504,9 @@
                  ElementwiseBinaryLowering<mg::SubOp, arith::SubFOp>,
                  ElementwiseBinaryLowering<mg::MulOp, arith::MulFOp>,
                  ElementwiseBinaryLowering<mg::DivOp, arith::DivFOp>,
-                 ScalarOpLowering, NegOpLowering, MatmulOpLowering,
+                 ScalarOpLowering, NegOpLowering,
                  ReluOpLowering, ReduceOpLowering, BroadcastOpLowering>(typeConverter, &getContext());
+    patterns.add<MatmulOpLowering>(typeConverter, &getContext(), matmulOrder == "ikj");
     populateFunctionOpInterfaceTypeConversionPattern<func::FuncOp>(patterns,
                                                                    typeConverter);
     populateCallOpTypeConversionPattern(patterns, typeConverter);
--- ../../part23/code/mgc   2026-10-02 16:31:18.251194699 +0000
+++ mgc 2026-10-02 18:08:07.897995992 +0000
@@ -7,27 +7,28 @@
 #   mgc ptx   prog.mg            compile to PTX (no launch: this machine has no GPU) and print the PTX
 # Options (Chapter 23): -O0|-O1|-O2|-O3        optimization level for the final clang step (default -O0, as in every earlier chapter)
 #   --passes "--affine-loop-tile=tile-size=32"  extra mg-opt passes run on the affine IR, before lowering (quote as one argument)
-# Env: MG_OPT (default Chapter 22 build), MGC_KEEP=dir (keep every intermediate file there), MGC_CFLAGS (extra flags for clang, e.g. -march=native),
+#   --matmul-order ikj|ijk                      (Chapter 24) the loop order mg.matmul is lowered to (default ijk)
+# Env: MG_OPT (default Chapter 24 build), MGC_KEEP=dir (keep every intermediate file there), MGC_CFLAGS (extra flags for clang, e.g. -march=native),
 #      MGC_CLANG (the clang to call; tests set it to `echo` to see exactly which flags mgc passes).
 set -e
 HERE=$(cd "$(dirname "$0")" && pwd)
-MG_OPT=${MG_OPT:-$HERE/../../part22/code/build/mg-opt}
+MG_OPT=${MG_OPT:-$HERE/../../part24/code/build/mg-opt}
 LIBDIR=${MLIR_LIBDIR:-/usr/lib/llvm-18/lib}
-usage() { echo "usage: mgc run|build|lib|mlir|ptx prog.mg [-o exe] [-O0..-O3] [--passes \"...\"]" >&2; exit 2; }
+usage() { echo "usage: mgc run|build|lib|mlir|ptx prog.mg [-o exe] [-O0..-O3] [--passes \"...\"] [--matmul-order ikj|ijk]" >&2; exit 2; }
 [ $# -ge 2 ] || usage
-cmd=$1; src=$2; shift 2; out=; OPT=; PASSES=
-while [ $# -gt 0 ]; do case $1 in -o) out=$2; shift 2;; -O0|-O1|-O2|-O3) OPT=$1; shift;; --passes) PASSES=$2; shift 2;; *) usage;; esac; done
+cmd=$1; src=$2; shift 2; out=; OPT=; PASSES=; CONVOPT=
+while [ $# -gt 0 ]; do case $1 in -o) out=$2; shift 2;; -O0|-O1|-O2|-O3) OPT=$1; shift;; --passes) PASSES=$2; shift 2;; --matmul-order) CONVOPT="=matmul-order=$2"; shift 2;; *) usage;; esac; done
 [ -f "$src" ] || { echo "mgc: no such file: $src" >&2; exit 2; }
 W=${MGC_KEEP:-$(mktemp -d)}; mkdir -p "$W"
 base=$(basename "$src" .mg)
 LIBFLAG=; [ "$cmd" = lib ] && LIBFLAG=--lib
 python3 "$HERE/../../part22/code/mgfront.py" $LIBFLAG "$src" > "$W/$base.mlir"          # 1. surface syntax -> mg dialect (exit 1 + message on a syntax/shape error)
 [ "$cmd" = mlir ] && { cat "$W/$base.mlir"; exit 0; }
-[ -x "$MG_OPT" ] || { echo "mgc: no mg-opt at $MG_OPT (run part22/code/build.sh)" >&2; exit 2; }
+[ -x "$MG_OPT" ] || { echo "mgc: no mg-opt at $MG_OPT (run part24/code/build.sh)" >&2; exit 2; }
 case $cmd in
   build|run|lib)
     LOWER='--lower-affine --convert-scf-to-cf --mg-lower-assert-to-stderr --convert-arith-to-llvm --finalize-memref-to-llvm --convert-func-to-llvm --reconcile-unrealized-casts'
-    $MG_OPT "$W/$base.mlir" --convert-mg-to-affine $PASSES -o "$W/$base.affine.mlir"   # 2. mg -> affine/memref/arith (+ runtime shape checks)
+    $MG_OPT "$W/$base.mlir" --convert-mg-to-affine$CONVOPT $PASSES -o "$W/$base.affine.mlir"   # 2. mg -> affine/memref/arith (+ runtime shape checks)
     $MG_OPT "$W/$base.affine.mlir" $LOWER -o "$W/$base.llvm.mlir"               # 3. -> LLVM dialect, assertions -> write(2,..)+abort
     mlir-translate-18 --mlir-to-llvmir "$W/$base.llvm.mlir" -o "$W/$base.ll"    # 4. -> LLVM IR
     if [ "$cmd" = lib ]; then   # object + header; the object is linked by the caller

To rebuild mg-opt:

#!/bin/sh
# Rebuild mg-opt for Chapter 24: Chapter 22's build with this directory's LowerToAffine.cpp (mg.matmul gains an `ikj` loop order,
# chosen with --convert-mg-to-affine=matmul-order=ikj; the default stays ijk). Produces ./build/mg-opt.
# Needs mlir-18-tools, libmlir-18-dev, llvm-18-dev, cmake, make.
set -e
HERE=$(cd "$(dirname "$0")" && pwd); D=$HERE/../..; T=$HERE/tree
mkdir -p $T/include/mg $T/lib $T/tools
cp $HERE/CMakeLists.txt $T/
cp $D/part3/code/MgDialect.td $D/part22/code/MgOps.td $D/part3/code/MgDialect.h $D/part3/code/MgOps.h $T/include/mg/
cp $D/part14/code/DynamicShapes.h $T/include/mg/
cp $D/part22/code/MgDialect.cpp $T/lib/
cp $D/part14/code/MgBufferizableOpInterfaceImpl.cpp $T/lib/
cp $D/part19/code/AssertToStderr.cpp $T/lib/
cp $HERE/LowerToAffine.cpp $T/lib/
cp $HERE/mg-opt.cpp $T/tools/
mkdir -p $HERE/build/include/mg && cd $HERE/build
cmake $T -DMLIR_DIR=/usr/lib/llvm-18/lib/cmake/mlir -DLLVM_DIR=/usr/lib/llvm-18/lib/cmake/llvm \
  -DCMAKE_BUILD_TYPE=Release -DLLVM_ENABLE_ASSERTIONS=OFF > cmake_out.txt 2>&1
make -j8 2>&1 | tail -2

You can see the order in the lowered code. This is a 2 by 3 times 3 by 4 product in each order; the accumulate nest has bounds 2, 4, 3 for ijk (the reduction k innermost) and 2, 3, 4 for ikj (the output column j innermost):

### matmul-order=ijk  (the zero-fill nest is omitted; this is the accumulate nest)
    affine.for %arg2 = 0 to 2 {
      affine.for %arg3 = 0 to 4 {
        affine.for %arg4 = 0 to 3 {
          %0 = affine.load %alloc[%arg2, %arg3] : memref<2x4xf64>
          %1 = affine.load %arg0[%arg2, %arg4] : memref<2x3xf64>
          %2 = affine.load %arg1[%arg4, %arg3] : memref<3x4xf64>
          %3 = arith.mulf %1, %2 : f64
          %4 = arith.addf %0, %3 : f64
          affine.store %4, %alloc[%arg2, %arg3] : memref<2x4xf64>
        }
      }
    }
### matmul-order=ikj  (the zero-fill nest is omitted; this is the accumulate nest)
    affine.for %arg2 = 0 to 2 {
      affine.for %arg3 = 0 to 3 {
        affine.for %arg4 = 0 to 4 {
          %0 = affine.load %alloc[%arg2, %arg4] : memref<2x4xf64>
          %1 = affine.load %arg0[%arg2, %arg3] : memref<2x3xf64>
          %2 = affine.load %arg1[%arg3, %arg4] : memref<3x4xf64>
          %3 = arith.mulf %1, %2 : f64
          %4 = arith.addf %0, %3 : f64
          affine.store %4, %alloc[%arg2, %arg4] : memref<2x4xf64>
        }
      }
    }

Results

The experiment is Chapter 23's, on the same kind of machine, with the new order added:

CPU:      Intel(R) Xeon(R) Processor @ 2.10GHz
CPU(s): 4
Thread(s) per core: 1
Core(s) per socket: 4
Socket(s): 1
L1d cache: 192 KiB (4 instances)
L2 cache: 8 MiB (4 instances)
L3 cache: 260 MiB (1 instance)
vector extensions available: avx avx2 avx512f fma sse4_2 
memory:   16094 MiB
kernel:   Linux 6.18.44-fc-v51
compiler: Ubuntu clang version 18.1.3 (1ubuntu1)
mlir:     Ubuntu LLVM version 18.1.3
This is a shared cloud virtual machine, not a quiet benchmarking box; see the chapter for what that means for the numbers.

As before: a shared cloud virtual machine, pinned to one core, two complete passes, every result compared with a reference by exact equality. The variants this time:

Name	What it is
mg_O2	Chapter 23's baseline: loop order ijk, clang -O2
mg_O2_ikj	the new order, -O2
mg_O2_ikj_tile32	the new order, tiled 32 by 32, -O2
mg_O3_native_tile32	Chapter 23's best dynamic-size variant: ijk, tiled, -O3 -march=native
mg_O3_native_ikj	the new order, -O3 -march=native
mg_O3_native_ikj_tile32	the new order, tiled, -O3 -march=native
C++ ijk loop, clang -O2, C++ ikj loop, -O3 native	the hand-written loops from Chapter 23, for comparison

The raw output of both passes (click to expand):

Pass 1: bench_out_1.txt

machine: 4 CPUs, Intel(R) Xeon(R) Processor @ 2.10GHz, Ubuntu clang version 18.1.3 (1ubuntu1)
pinned with: taskset -c 3   load average before: 0.86 0.43 0.17
mg_O2 (dynamic sizes)      N=64    min     0.094 ms  median     0.104 ms     5.06 GFLOP/s (median)  result ok
mg_O2 (dynamic sizes)      N=128   min     2.211 ms  median     2.257 ms     1.86 GFLOP/s (median)  result ok
mg_O2 (dynamic sizes)      N=256   min    17.219 ms  median    17.788 ms     1.89 GFLOP/s (median)  result ok
mg_O2 (static 256)         N=256   min    16.667 ms  median    16.867 ms     1.99 GFLOP/s (median)  result ok
mg_O2 (dynamic sizes)      N=512   min   443.277 ms  median   490.884 ms     0.55 GFLOP/s (median)  result ok
mg_O2_ikj (dynamic sizes)  N=64    min     0.063 ms  median     0.064 ms     8.13 GFLOP/s (median)  result ok
mg_O2_ikj (dynamic sizes)  N=128   min     0.587 ms  median     0.646 ms     6.49 GFLOP/s (median)  result ok
mg_O2_ikj (dynamic sizes)  N=256   min     6.264 ms  median     6.296 ms     5.33 GFLOP/s (median)  result ok
mg_O2_ikj (static 256)     N=256   min     3.636 ms  median     3.645 ms     9.21 GFLOP/s (median)  result ok
mg_O2_ikj (dynamic sizes)  N=512   min    43.697 ms  median    61.176 ms     4.39 GFLOP/s (median)  result ok
mg_O2_ikj_tile32 (dynamic sizes) N=64    min     0.060 ms  median     0.061 ms     8.65 GFLOP/s (median)  result ok
mg_O2_ikj_tile32 (dynamic sizes) N=128   min     0.470 ms  median     0.470 ms     8.92 GFLOP/s (median)  result ok
mg_O2_ikj_tile32 (dynamic sizes) N=256   min     6.272 ms  median     7.019 ms     4.78 GFLOP/s (median)  result ok
mg_O2_ikj_tile32 (static 256) N=256   min     4.967 ms  median     5.065 ms     6.63 GFLOP/s (median)  result ok
mg_O2_ikj_tile32 (dynamic sizes) N=512   min    38.142 ms  median    38.334 ms     7.00 GFLOP/s (median)  result ok
mg_O3_native_tile32 (dynamic sizes) N=64    min     0.100 ms  median     0.103 ms     5.07 GFLOP/s (median)  result ok
mg_O3_native_tile32 (dynamic sizes) N=128   min     1.239 ms  median     1.309 ms     3.20 GFLOP/s (median)  result ok
mg_O3_native_tile32 (dynamic sizes) N=256   min    11.280 ms  median    11.903 ms     2.82 GFLOP/s (median)  result ok
mg_O3_native_tile32 (static 256) N=256   min     3.388 ms  median     3.478 ms     9.65 GFLOP/s (median)  result ok
mg_O3_native_tile32 (dynamic sizes) N=512   min    78.325 ms  median    81.522 ms     3.29 GFLOP/s (median)  result ok
mg_O3_native_ikj (dynamic sizes) N=64    min     0.035 ms  median     0.042 ms    12.50 GFLOP/s (median)  result ok
mg_O3_native_ikj (dynamic sizes) N=128   min     0.491 ms  median     0.509 ms     8.24 GFLOP/s (median)  result ok
mg_O3_native_ikj (dynamic sizes) N=256   min     3.865 ms  median     3.939 ms     8.52 GFLOP/s (median)  result ok
mg_O3_native_ikj (static 256) N=256   min     3.370 ms  median     3.580 ms     9.37 GFLOP/s (median)  result ok
mg_O3_native_ikj (dynamic sizes) N=512   min    34.144 ms  median    36.563 ms     7.34 GFLOP/s (median)  result ok
mg_O3_native_ikj_tile32 (dynamic sizes) N=64    min     0.061 ms  median     0.068 ms     7.76 GFLOP/s (median)  result ok
mg_O3_native_ikj_tile32 (dynamic sizes) N=128   min     0.366 ms  median     0.418 ms    10.04 GFLOP/s (median)  result ok
mg_O3_native_ikj_tile32 (dynamic sizes) N=256   min     4.007 ms  median     4.863 ms     6.90 GFLOP/s (median)  result ok
mg_O3_native_ikj_tile32 (static 256) N=256   min     2.601 ms  median     2.835 ms    11.84 GFLOP/s (median)  result ok
mg_O3_native_ikj_tile32 (dynamic sizes) N=512   min    38.855 ms  median    42.166 ms     6.37 GFLOP/s (median)  result ok
C++ ijk loop, clang -O2    N=64    min     0.130 ms  median     0.136 ms     3.86 GFLOP/s (median)  result ok
C++ ijk loop, clang -O2    N=128   min     2.337 ms  median     2.424 ms     1.73 GFLOP/s (median)  result ok
C++ ijk loop, clang -O2    N=256   min    17.151 ms  median    17.407 ms     1.93 GFLOP/s (median)  result ok
C++ ijk loop, clang -O2    N=512   min   435.574 ms  median   489.424 ms     0.55 GFLOP/s (median)  result ok
C++ ikj loop, -O3 native   N=64    min     0.037 ms  median     0.038 ms    13.95 GFLOP/s (median)  result ok
C++ ikj loop, -O3 native   N=128   min     0.428 ms  median     0.443 ms     9.47 GFLOP/s (median)  result ok
C++ ikj loop, -O3 native   N=256   min     4.076 ms  median     4.534 ms     7.40 GFLOP/s (median)  result ok
C++ ikj loop, -O3 native   N=512   min    33.838 ms  median    35.939 ms     7.47 GFLOP/s (median)  result ok
load average after: 0.91 0.48 0.19

Pass 2: bench_out_2.txt

machine: 4 CPUs, Intel(R) Xeon(R) Processor @ 2.10GHz, Ubuntu clang version 18.1.3 (1ubuntu1)
pinned with: taskset -c 3   load average before: 0.92 0.49 0.20
mg_O2 (dynamic sizes)      N=64    min     0.094 ms  median     0.094 ms     5.58 GFLOP/s (median)  result ok
mg_O2 (dynamic sizes)      N=128   min     2.156 ms  median     2.208 ms     1.90 GFLOP/s (median)  result ok
mg_O2 (dynamic sizes)      N=256   min    16.735 ms  median    17.191 ms     1.95 GFLOP/s (median)  result ok
mg_O2 (static 256)         N=256   min    16.598 ms  median    16.915 ms     1.98 GFLOP/s (median)  result ok
mg_O2 (dynamic sizes)      N=512   min   435.352 ms  median   463.103 ms     0.58 GFLOP/s (median)  result ok
mg_O2_ikj (dynamic sizes)  N=64    min     0.060 ms  median     0.061 ms     8.56 GFLOP/s (median)  result ok
mg_O2_ikj (dynamic sizes)  N=128   min     0.563 ms  median     0.580 ms     7.23 GFLOP/s (median)  result ok
mg_O2_ikj (dynamic sizes)  N=256   min     3.942 ms  median     5.028 ms     6.67 GFLOP/s (median)  result ok
mg_O2_ikj (static 256)     N=256   min     4.211 ms  median     5.394 ms     6.22 GFLOP/s (median)  result ok
mg_O2_ikj (dynamic sizes)  N=512   min    62.850 ms  median    66.111 ms     4.06 GFLOP/s (median)  result ok
mg_O2_ikj_tile32 (dynamic sizes) N=64    min     0.085 ms  median     0.092 ms     5.71 GFLOP/s (median)  result ok
mg_O2_ikj_tile32 (dynamic sizes) N=128   min     0.649 ms  median     0.674 ms     6.22 GFLOP/s (median)  result ok
mg_O2_ikj_tile32 (dynamic sizes) N=256   min     4.697 ms  median     4.941 ms     6.79 GFLOP/s (median)  result ok
mg_O2_ikj_tile32 (static 256) N=256   min     3.164 ms  median     3.334 ms    10.06 GFLOP/s (median)  result ok
mg_O2_ikj_tile32 (dynamic sizes) N=512   min    39.203 ms  median    42.462 ms     6.32 GFLOP/s (median)  result ok
mg_O3_native_tile32 (dynamic sizes) N=64    min     0.103 ms  median     0.179 ms     2.93 GFLOP/s (median)  result ok
mg_O3_native_tile32 (dynamic sizes) N=128   min     0.773 ms  median     0.875 ms     4.79 GFLOP/s (median)  result ok
mg_O3_native_tile32 (dynamic sizes) N=256   min     7.778 ms  median     8.036 ms     4.18 GFLOP/s (median)  result ok
mg_O3_native_tile32 (static 256) N=256   min     3.292 ms  median     3.706 ms     9.05 GFLOP/s (median)  result ok
mg_O3_native_tile32 (dynamic sizes) N=512   min    85.210 ms  median   110.995 ms     2.42 GFLOP/s (median)  result ok
mg_O3_native_ikj (dynamic sizes) N=64    min     0.035 ms  median     0.044 ms    11.80 GFLOP/s (median)  result ok
mg_O3_native_ikj (dynamic sizes) N=128   min     0.408 ms  median     0.450 ms     9.32 GFLOP/s (median)  result ok
mg_O3_native_ikj (dynamic sizes) N=256   min     3.301 ms  median     3.450 ms     9.73 GFLOP/s (median)  result ok
mg_O3_native_ikj (static 256) N=256   min     3.190 ms  median     3.332 ms    10.07 GFLOP/s (median)  result ok
mg_O3_native_ikj (dynamic sizes) N=512   min    37.478 ms  median    38.669 ms     6.94 GFLOP/s (median)  result ok
mg_O3_native_ikj_tile32 (dynamic sizes) N=64    min     0.049 ms  median     0.059 ms     8.93 GFLOP/s (median)  result ok
mg_O3_native_ikj_tile32 (dynamic sizes) N=128   min     0.352 ms  median     0.361 ms    11.62 GFLOP/s (median)  result ok
mg_O3_native_ikj_tile32 (dynamic sizes) N=256   min     3.759 ms  median     3.970 ms     8.45 GFLOP/s (median)  result ok
mg_O3_native_ikj_tile32 (static 256) N=256   min     2.417 ms  median     2.714 ms    12.36 GFLOP/s (median)  result ok
mg_O3_native_ikj_tile32 (dynamic sizes) N=512   min    32.267 ms  median    35.527 ms     7.56 GFLOP/s (median)  result ok
C++ ijk loop, clang -O2    N=64    min     0.092 ms  median     0.104 ms     5.05 GFLOP/s (median)  result ok
C++ ijk loop, clang -O2    N=128   min     2.219 ms  median     2.227 ms     1.88 GFLOP/s (median)  result ok
C++ ijk loop, clang -O2    N=256   min    16.938 ms  median    17.020 ms     1.97 GFLOP/s (median)  result ok
C++ ijk loop, clang -O2    N=512   min   404.418 ms  median   500.284 ms     0.54 GFLOP/s (median)  result ok
C++ ikj loop, -O3 native   N=64    min     0.035 ms  median     0.038 ms    13.74 GFLOP/s (median)  result ok
C++ ikj loop, -O3 native   N=128   min     0.419 ms  median     0.428 ms     9.80 GFLOP/s (median)  result ok
C++ ikj loop, -O3 native   N=256   min     3.256 ms  median     3.348 ms    10.02 GFLOP/s (median)  result ok
C++ ikj loop, -O3 native   N=512   min    34.176 ms  median    36.967 ms     7.26 GFLOP/s (median)  result ok
load average after: 0.95 0.54 0.23

And the combined table (median time per call in each pass, how much the passes differ, average GFLOP/s, and speedup versus mg_O2 at the same size; every result matched the reference):

implementation                        N  pass 1 ms  pass 2 ms  differ  GFLOP/s vs mg -O2  results
mg_O2 (dynamic sizes)                64      0.104      0.094     10%     5.32     1.00x  ok
mg_O2_ikj (dynamic sizes)            64      0.064      0.061      5%     8.35     1.58x  ok
mg_O2_ikj_tile32 (dynamic sizes)     64      0.061      0.092     41%     7.18     1.29x  ok
mg_O3_native_tile32 (dynamic sizes)   64      0.103      0.179     54%     4.00     0.70x  ok
mg_O3_native_ikj (dynamic sizes)     64      0.042      0.044      5%    12.15     2.30x  ok
mg_O3_native_ikj_tile32 (dynamic sizes)   64      0.068      0.059     14%     8.34     1.56x  ok
C++ ijk loop, clang -O2              64      0.136      0.104     27%     4.46     0.83x  ok
C++ ikj loop, -O3 native             64      0.038      0.038      0%    13.84     2.61x  ok

mg_O2 (dynamic sizes)               128      2.257      2.208      2%     1.88     1.00x  ok
mg_O2_ikj (dynamic sizes)           128      0.646      0.580     11%     6.86     3.64x  ok
mg_O2_ikj_tile32 (dynamic sizes)    128      0.470      0.674     36%     7.57     3.90x  ok
mg_O3_native_tile32 (dynamic sizes)  128      1.309      0.875     40%     4.00     2.04x  ok
mg_O3_native_ikj (dynamic sizes)    128      0.509      0.450     12%     8.78     4.66x  ok
mg_O3_native_ikj_tile32 (dynamic sizes)  128      0.418      0.361     15%    10.83     5.73x  ok
C++ ijk loop, clang -O2             128      2.424      2.227      8%     1.80     0.96x  ok
C++ ikj loop, -O3 native            128      0.443      0.428      3%     9.64     5.13x  ok

mg_O2 (dynamic sizes)               256     17.788     17.191      3%     1.92     1.00x  ok
mg_O2 (static 256)                  256     16.867     16.915      0%     1.98     1.04x  ok
mg_O2_ikj (dynamic sizes)           256      6.296      5.028     22%     6.00     3.09x  ok
mg_O2_ikj (static 256)              256      3.645      5.394     39%     7.71     3.87x  ok
mg_O2_ikj_tile32 (dynamic sizes)    256      7.019      4.941     35%     5.79     2.92x  ok
mg_O2_ikj_tile32 (static 256)       256      5.065      3.334     41%     8.35     4.16x  ok
mg_O3_native_tile32 (dynamic sizes)  256     11.903      8.036     39%     3.50     1.75x  ok
mg_O3_native_tile32 (static 256)    256      3.478      3.706      6%     9.35     4.87x  ok
mg_O3_native_ikj (dynamic sizes)    256      3.939      3.450     13%     9.12     4.73x  ok
mg_O3_native_ikj (static 256)       256      3.580      3.332      7%     9.72     5.06x  ok
mg_O3_native_ikj_tile32 (dynamic sizes)  256      4.863      3.970     20%     7.67     3.96x  ok
mg_O3_native_ikj_tile32 (static 256)  256      2.835      2.714      4%    12.10     6.30x  ok
C++ ijk loop, clang -O2             256     17.407     17.020      2%     1.95     1.02x  ok
C++ ikj loop, -O3 native            256      4.534      3.348     30%     8.71     4.44x  ok

mg_O2 (dynamic sizes)               512    490.884    463.103      6%     0.56     1.00x  ok
mg_O2_ikj (dynamic sizes)           512     61.176     66.111      8%     4.22     7.49x  ok
mg_O2_ikj_tile32 (dynamic sizes)    512     38.334     42.462     10%     6.66    11.81x  ok
mg_O3_native_tile32 (dynamic sizes)  512     81.522    110.995     31%     2.85     4.96x  ok
mg_O3_native_ikj (dynamic sizes)    512     36.563     38.669      6%     7.14    12.68x  ok
mg_O3_native_ikj_tile32 (dynamic sizes)  512     42.166     35.527     17%     6.96    12.28x  ok
C++ ijk loop, clang -O2             512    489.424    500.284      2%     0.55     0.96x  ok
C++ ikj loop, -O3 native            512     35.939     36.967      3%     7.37    13.09x  ok

Reading the numbers

1. The new loop order alone is a large win. At -O2, ikj is 3.6× faster than ijk at N = 128, 3.1× at N = 256 and 7.5× at N = 512 (4.22 against 0.56 GFLOP/s). Nothing about the arithmetic changed.

2. With the CPU's full instruction set, the gap from Chapter 23 closes. Chapter 23's result 6 said the best Mountain Goat variant at N = 512 was about 2.1× slower than a hand-written ikj loop. Now mg_O3_native_ikj reaches 7.14 GFLOP/s against the hand-written loop's 7.37, within the noise (the two differ by 6% and 3% between passes). At N = 256 with dynamic sizes it reaches 9.12 against the hand-written loop's 8.71, but that hand-written figure moved by 30% between passes, so the right reading is "comparable", not "faster".

3. Tiling no longer clearly helps. On top of ikj with all instructions, tiling gives 6.96 against 7.14 GFLOP/s at N = 512 and 7.67 against 9.12 at N = 256 with dynamic sizes: no gain, and the passes differ by up to 20%. At N = 128 the tiled version is faster (10.83 against 8.78), but the passes differ by 12% and 15%, so it is not a firm result. The loop order had already fixed most of what tiling was fixing in Chapter 23. There are two exceptions, and both are specific: at plain -O2, tiling still helps ikj (6.66 against 4.22 GFLOP/s at N = 512, 1.6×), and with a fixed shape at N = 256 the tiled ikj native build reaches 12.10 GFLOP/s against 9.72 untiled, the best number in the table and the only one clearly above the hand-written loop (and that comparison is at one size, against a hand-written figure that was itself noisy).

4. It is the loop order that unlocks vectorization. Chapter 23 left this open (it had only seen vector instructions in one build and could not say which ingredient mattered). inspect_asm.sh counts them again:

ijk, untiled, mm_dyn: 0 vector/FMA instructions
ijk, untiled, mm_256: 0 vector/FMA instructions
ijk, tiled (32), mm_dyn: 0 vector/FMA instructions
ijk, tiled (32), mm_256: 64 vector/FMA instructions
ikj, untiled, mm_dyn: 8 vector/FMA instructions
ikj, untiled, mm_256: 16 vector/FMA instructions
ikj, tiled (32), mm_dyn: 8 vector/FMA instructions
ikj, tiled (32), mm_256: 32 vector/FMA instructions

With ijk, the only function with vector instructions is the tiled fixed-shape one (64). With ikj, even the untiled dynamic-size function has them (8), and the fixed-shape ones have more (16 untiled, 32 tiled). So for dynamic sizes the loop order alone is enough; the fixed shape and the tiling are not what was needed. The consistent explanation is that the ikj inner loop reads and writes consecutive addresses, which is what a vectorizer wants. I did not test other loop orders (for example jki) or look at the vectorizer's own reports, so "consistent" is as far as this goes.

5. The numbers are noisy, and some spreads are large. Between the two passes the same experiment differs by up to 54% (mg_O3_native_tile32 at N = 64), 41% and 39% at N = 256. As in Chapter 23, differences of 2× and more are claims; differences of 1.0× to 1.3× are not.

Same answers, to the last bit

The argument above says the two loop orders add each element's products in the same order. bits_check.sh tests it with numbers where it would show: a 5 by 7 times 7 by 4 product of values that are not exactly representable, printed as hexadecimal floating point (every bit of every result), built once with each loop order:

// Prints every element of a matrix product in hexadecimal floating point (%a), which shows every bit. The inputs are not
// "nice" numbers (they are not exactly representable and the products round), so a different order of additions would change
// bits. Built once against an ijk library and once against an ikj library; the two outputs must be identical.
#include <cstdio>
#include "ops.h"
int main() {
  const int m = 5, k = 7, n = 4;
  std::vector<double> a(m * k), b(k * n);
  for (int i = 0; i < m * k; i++) a[i] = 0.1 * (i % 13) + 1.0 / (i + 3);
  for (int i = 0; i < k * n; i++) b[i] = 0.7 / (i % 11 + 1) - 0.01 * i;
  mg::Matrix r = mg::matmul(mg::Matrix(m, k, a), mg::Matrix(k, n, b));
  for (double v : r.data) std::printf("%a\n", v);
  return 0;
}

ijk and ikj outputs: identical (20 elements, every bit)
compiled ikj output equals an independent Python sum in increasing k: True
the same sum in DECREASING k order differs from it in 19 of 20 elements (so this check would notice a changed order)

The two outputs are identical in all 20 elements. The script also checks that the check is able to fail: summing the same products in decreasing k order, instead of increasing, changes 19 of the 20 results, so a loop order that changed the order of additions would have been caught. (The benchmark itself uses small whole numbers, where every sum is exact in any order, so on its own it could not have shown this.)

Tests

Six new test files in test/loop-order/ (the suite now has 99 tests). As in Chapter 23, none asserts speed.

Test	What it checks
structure	the default and matmul-order=ijk give loop bounds 2, 4, 3 for a 2×3 @ 3×4 product, and ikj gives 2, 3, 4 (with the zero-fill nest before it)
bad-option	matmul-order=jki is an error with a clear message, not a silent fallback
results-identical	the same 3×5 @ 5×2 product (no dimension a multiple of the tile size) prints the same at ijk, ikj, and ikj tiled by 2 and by 3, static and dynamic, equal to the hand-computed 12 15 / 32 35 / 52 55
cpp-ikj	Chapter 21's whole C++ program (every elementwise result, matmul against a plain C++ loop, the dirty-heap check, the abort message) against a library built with ikj
flag-reaches-pass	mgc --matmul-order really reaches the pass: the kept affine IR differs in loop order
bitwise-identical	the hexadecimal-float comparison above, as a test

Two bugs were injected into mgc and four into the lowering, one at a time:

#!/usr/bin/env python3
# READ THIS FIRST: this script injects deliberate bugs into the mgc driver and re-runs the Chapter 24 tests (test/loop-order).
#   "caught by <tests>" lines are the EXPECTED, wanted result: a test noticed the bug. "NOT CAUGHT" would be a test gap.
#   mgc is restored afterwards. Output: mgc_mutation_out.txt
import os, subprocess, sys
HERE = os.path.dirname(os.path.abspath(__file__)); MGC = os.path.join(HERE, "mgc")
LIT = os.path.join(HERE, "..", "..", "part15", "code", "run_lit.sh")
print("NOTE: this script deliberately breaks mgc. 'caught by' lines are EXPECTED: they show the tests can detect the bug.")
def run(label, baseline=False):
    r = subprocess.run(["sh", LIT, "--filter", ":: loop-order/"], capture_output=True, text=True,
                       env={**os.environ, "MG_OPT": os.path.join(HERE, "build", "mg-opt"), "MG_TEST_TMP": os.path.join(HERE, "work", "perfmut")})
    fails = [l.split("::")[1].strip().split(" (")[0] for l in r.stdout.splitlines() if l.startswith("FAIL")]
    if baseline: print(f"{label}: " + ("all tests pass (expected)" if not fails else "UNEXPECTED FAILURES: " + ", ".join(fails)))
    else: print(f"{label}: " + ("caught by " + ", ".join(fails) if fails else "NOT CAUGHT (a test gap)"))
orig = open(MGC).read()
base = "run baseline (no bug)"
run(base, baseline=True)
muts = [("--matmul-order ignored by mgc", '--convert-mg-to-affine$CONVOPT $PASSES -o', '--convert-mg-to-affine $PASSES -o'),
        ("--matmul-order value passed to the wrong option", 'CONVOPT="=matmul-order=$2"', 'CONVOPT="=matmul-order=ijk"')]
try:
    for label, a, b in muts:
        if a not in orig: print(label, ": MUTATION DID NOT APPLY"); continue
        open(MGC, "w").write(orig.replace(a, b, 1)); run(label)
finally:
    open(MGC, "w").write(orig)
run("restored (baseline again)", baseline=True)

NOTE: this script deliberately breaks mgc. 'caught by' lines are EXPECTED: they show the tests can detect the bug.
run baseline (no bug): all tests pass (expected)
--matmul-order ignored by mgc: caught by loop-order/flag-reaches-pass.mlir
--matmul-order value passed to the wrong option: caught by loop-order/flag-reaches-pass.mlir
restored (baseline again): all tests pass (expected)

#!/bin/sh
# Rebuild mg-opt with one bug injected into LowerToAffine.cpp and run the ops tests against it.
# Each build is a full rebuild in work/mutN (about two minutes). Output: lowering_mutation_out.txt
# READ THIS FIRST: this script rebuilds mg-opt with one deliberate bug injected into the lowering at a time (a full rebuild each) and runs the tests.
#   FAIL lines (and non-zero exit statuses) in this script's output are the EXPECTED result: each bug SHOULD make at least one test fail (a FAIL line means a test caught it). A section with NO FAIL line would be a test gap.
#   The unmodified, current build is the baseline and must show no failures. A mutation or old build that makes NO test fail would be the problem.
echo "NOTE: this script deliberately runs broken or older code. FAIL lines below are EXPECTED: they show the tests can detect the problem. The baseline (unmodified current build) must show none."
HERE=$(cd "$(dirname "$0")" && pwd); LIT=$HERE/../../part15/code; M=$HERE/work/mut; rm -rf $M; mkdir -p $M
n=0
try() {  # label sed-expr
  n=$((n+1)); d=$M/$n; mkdir -p $d/build; cp -r $HERE/tree $d/tree; rm -rf $d/tree/build
  cp $d/tree/lib/LowerToAffine.cpp $d/before.cpp; sed -i "$2" $d/tree/lib/LowerToAffine.cpp
  if cmp -s $d/before.cpp $d/tree/lib/LowerToAffine.cpp; then echo "### $1: MUTATION DID NOT APPLY"; return; fi
  ( cd $d/build && cmake $d/tree -DMLIR_DIR=/usr/lib/llvm-18/lib/cmake/mlir -DLLVM_DIR=/usr/lib/llvm-18/lib/cmake/llvm -DCMAKE_BUILD_TYPE=Release -DLLVM_ENABLE_ASSERTIONS=OFF >/dev/null 2>&1 && mkdir -p include/mg && make -j8 >/dev/null 2>&1 )
  [ -x $d/build/mg-opt ] || { echo "### $1: BUILD FAILED"; return; }
  printf '%s\n' "### lowering mutation: $1"
  MG_OPT=$d/build/mg-opt MG_TEST_TMP=$d/out sh $LIT/run_lit.sh --filter ':: loop-order/' 2>&1 | grep -E "^(FAIL)|  Passed|  Failed" | sed 's#IR Forge / Mountain Goat :: ##'
}
try "the ikj order is ignored (always lowers as ijk)"   's/: OpConversionPattern(converter, ctx), ikj(ikj) {}/: OpConversionPattern(converter, ctx), ikj(false) {}/'
try "ikj bounds not reordered (m, n, k used for i, k, j)"  's/ikj ? llvm::SmallVector<Value, 3>{m, k, n} : llvm::SmallVector<Value, 3>{m, n, k}/llvm::SmallVector<Value, 3>{m, n, k}/'
try "the right-hand matrix is read transposed"           's/Value y = b.create<affine::AffineLoadOp>(l, rhs, ValueRange{kk, j});/Value y = b.create<affine::AffineLoadOp>(l, rhs, ValueRange{j, kk});/'
try "the matmul-order option is not validated"           's/if (matmulOrder != "ijk" \&\& matmulOrder != "ikj") {/if (false) {/'

NOTE: this script deliberately runs broken or older code. FAIL lines below are EXPECTED: they show the tests can detect the problem. The baseline (unmodified current build) must show none.
### lowering mutation: the ikj order is ignored (always lowers as ijk)
FAIL: loop-order/flag-reaches-pass.mlir (2 of 6)
FAIL: loop-order/structure.mlir (3 of 6)
  Passed  :  4 (4.04%)
  Failed  :  2 (2.02%)
### lowering mutation: ikj bounds not reordered (m, n, k used for i, k, j)
FAIL: loop-order/flag-reaches-pass.mlir (2 of 6)
FAIL: loop-order/results-identical.mlir (3 of 6)
FAIL: loop-order/structure.mlir (4 of 6)
FAIL: loop-order/cpp-ikj.mlir (5 of 6)
FAIL: loop-order/bitwise-identical.mlir (6 of 6)
  Passed  :  1 (1.01%)
  Failed  :  5 (5.05%)
### lowering mutation: the right-hand matrix is read transposed
FAIL: loop-order/cpp-ikj.mlir (4 of 6)
FAIL: loop-order/results-identical.mlir (6 of 6)
  Passed  :  4 (4.04%)
  Failed  :  2 (2.02%)
### lowering mutation: the matmul-order option is not validated
FAIL: loop-order/bad-option.mlir (1 of 6)
  Passed  :  5 (5.05%)
  Failed  :  1 (1.01%)

Reading the mutation results:

• "The ikj order is ignored" produces correct answers (the two orders are identical by design), so only tests that look at the structure of the generated loops can catch it: structure and flag-reaches-pass do, and nothing else fails.
• "Bounds not reordered" (loops run in i, k, j order but with the i, j, k sizes) is caught by five of the six tests: it produces wrong answers or reads outside the matrices.
• "The right-hand matrix is read transposed" is caught by only two tests, results-identical (it checks a hand-computed answer) and cpp-ikj. It is not caught by bitwise-identical, and that is worth understanding: that test compares the ijk result with the ikj result, and this bug is in the code both orders share, so both are equally wrong and still agree with each other. A differential test (compare two implementations) cannot catch a bug that both implementations have; a test against a known right answer can.
• "The option is not validated" is caught by exactly one test, bad-option, the only test that passes a bad value.

The full suite with the Chapter 24 build:

Testing Time: 6.37s

Total Discovered Tests: 99
  Passed: 99 (100.00%)

Limits and what is not established

• Same machine, same caveats as Chapter 23: a shared cloud VM, one core, one session, noisy (up to 54% between passes here).
• Still one operation. Square f64 products from 64 to 512; nothing about other shapes (a tall-thin product has different memory behavior), and not about the other operations.
• Two loop orders only. ijk and ikj. The other four orders were not tried, so "the loop order matters" is shown for this pair, not as a ranking of all six.
• One tile size (32), still untuned, and tiling was measured on top of ikj only at that size.
• No hardware counters and no vectorizer reports. The explanations (contiguous inner loop, cache lines) are consistent with the data and with the instruction counts, not demonstrated by measuring cache misses.
• The default is still ijk. The faster order is opt-in (--matmul-order ikj), because changing the default would change every earlier chapter's generated code and recorded output. A real compiler would make the better order the default.
• No multithreading, no GPU measurement, and nothing like a tuned library's blocking, packing or register micro-kernel; the best result here is still modest next to what a tuned library reaches on this CPU.
• The benchmark's exact-equality check uses small whole numbers; the bit-for-bit claim rests on the separate hexadecimal-float check (one 5×7×4 product), on the argument about the order of additions, and on that check being shown able to fail.

Reproducing

cd docs/part24/code
./build.sh
./run_bench.sh > bench_out_1.txt && ./run_bench.sh > bench_out_2.txt
python3 summarize.py bench_out_1.txt bench_out_2.txt > summary.txt
./inspect_asm.sh > inspect_asm_out.txt && ./bits_check.sh > bits_check_out.txt && ./machine_info.sh > machine_info_out.txt
./mgc_mutation.py > mgc_mutation_out.txt
./lowering_mutation.sh > lowering_mutation_out.txt     # four rebuilds
cd ../../part15/code && ./run_lit.sh                    # 99 tests

Chapter summary

• A second loop order for mg.matmul, ikj, is selected with --convert-mg-to-affine=matmul-order=ikj or mgc --matmul-order ikj; the default (ijk) is unchanged.
• It is 3× to 7.5× faster than ijk at -O2 on this machine, and with all CPU instructions it matches a hand-written ikj loop at N = 512 (7.14 against 7.37 GFLOP/s), closing the gap Chapter 23 left.
• Tiling on top of ikj helps little or not at all, except at -O2 and for one fixed shape; the loop order did most of the work.
• Loop order alone is enough for the compiler to produce vector instructions for dynamic sizes; Chapter 23's open question is answered for this pair of orders.
• The two orders give bit-identical results, shown by a hexadecimal-float test that can fail; the benchmark's own check could not have shown it.
• Six tests (suite: 99) and six injected bugs, all caught.

Self-check questions

Each answer is collapsed; try the question first.

1. In ijk and ikj, which loop is innermost, and which memory does that innermost loop walk?

   Answer

   In ijk the innermost loop is k: it reads A[i][k] along a row of A (consecutive addresses) and B[k][j] down a column of B (jumping a whole row each step). In ikj the innermost loop is j: it reads B[k][j] and updates C[i][j], and both walk along a row, so consecutive addresses.

2. Why does walking down a column of a row-major matrix cost more than walking along a row?

   Answer

   Memory is fetched in cache lines of about 64 bytes (eight doubles). Along a row, every number of a fetched line is used. Down a column, each step lands in a different line, so one number of each fetched line is used and the rest is wasted, and for large matrices those lines are evicted before they are needed again. A walk along a row also lets the compiler use vector instructions; a walk with a large stride does not.

3. Why are the ijk and ikj results identical to the last bit, rather than merely close?

   Answer

   Floating-point addition is not associative, so the order of additions matters. But in both loop orders, each result element C[i][j] accumulates its products in the same order (k increasing from 0), so each element goes through exactly the same sequence of additions and the same roundings. What differs between the orders is only the order in which different elements are worked on, which does not affect any element's value.

4. Why does the benchmark's exact-equality check not prove that, and what does?

   Answer

   The benchmark uses small whole numbers, for which every product and sum is exact in any order, so any order would pass. The proof needs numbers that round: bits_check.sh multiplies values that are not exactly representable, prints every bit, and compares the two orders. It also shows the comparison could fail: adding in decreasing k order changes 19 of the 20 results.

5. Chapter 23 could not say whether tiling or fixed sizes unlocked vectorization. What did this chapter find, and what did it not test?

   Answer

   With the ikj order, even the untiled dynamic-size function contains vector instructions (8, against 0 with ijk), so for dynamic sizes the loop order alone is sufficient; neither tiling nor fixed shape was needed. It did not test other loop orders, did not read the compiler's vectorization reports, and did not measure why the fixed-shape builds have more vector instructions. The explanation (consecutive addresses in the inner loop) is consistent with the evidence, not proven.

6. Tiling gave 2.6× to 5.2× in Chapter 23 but little here. Why, and where does it still help?

   Answer

   Part of what tiling does is reduce cache misses caused by the column walk. The ikj order removed the column walk, so there is less left for tiling to fix. It still helps at plain -O2 (1.6× at N = 512) and with a fixed 256 shape and all instructions (12.10 against 9.72 GFLOP/s), where other effects are in play that this chapter did not isolate.

7. Why is the default still ijk if ikj is better?

   Answer

   Changing the default would change the generated code for every earlier chapter's matmul, and with it recorded outputs, structure tests and page text. Making the better order opt-in kept everything earlier valid. A production compiler would make the better order the default and keep the old one selectable, or choose per shape.

8. The test bad-option checks that matmul-order=jki is an error. What would the mutation "option not validated" do without it?

   Answer

   The pass would accept any string and treat everything other than ikj as ijk, so jki (or a typo like ijk or IKJ) would silently produce the default order. The user would think they had changed the order and would see no change. The test makes a wrong value fail loudly.

9. The mutation "the ikj order is ignored (always lowers as ijk)" produces correct results. How can any test catch it?

   Answer

   Only by looking at the structure of the lowered code, because the results are identical by design. structure and flag-reaches-pass check the loop bounds in the generated IR (2, 3, 4 for ikj). This is the same situation as Chapter 22's equivalent mutant: a bug with no effect on results can only be caught by checking implementation, not behavior.

10. What would you have to do to claim that ikj is "the" right loop order for Mountain Goat?

    Answer

    Try all six orders, many shapes (including tall, thin and small ones), more machines and compilers, and look at hardware counters and vectorizer reports to confirm the mechanism. Then compare against a tuned library. This chapter measured two orders on square products on one shared virtual machine.