Matrix-free CG Hands-On Walkthrough (Fritz Icelake Node!)
MFCG Hands-on Solution (performed on Fritz Icelake node)
Time profiling
Perform runtime profile using gprof:
icx -Ofast -xhost -std=c99 -pg -o ./mfcg ./mfcg.c
Execute code:
likwid-pin -c S0:2 ./mfcg 8000 20000
Examine profile:
gprof ./mfcg
Flat profile:
Each sample counts as 0.01 seconds.
% cumulative self self total
time seconds seconds calls s/call s/call name
92.90 46.85 46.85 1 46.85 47.39 CG
4.76 49.25 2.40 main
1.07 49.79 0.54 1 0.54 0.54 computeResidual
0.68 50.14 0.35 __svml_sin4_l9
0.54 50.41 0.27 __intel_avx_rep_memcpy
0.07 50.44 0.04 __ocl_svml_h8__svml_dsin_cout_rare_internal_wrapper
Problem: Everything got inlined!
Solution: Add -fno-inline flag.
gprof ./mfcg
Flat profile:
Each sample counts as 0.01 seconds.
% cumulative self self total
time seconds seconds calls s/call s/call name
53.00 24.94 24.94 155 0.16 0.16 axpby
26.09 37.21 12.27 105 0.12 0.12 dotProduct
18.14 45.75 8.53 53 0.16 0.16 applyStencil
1.79 46.59 0.84 1 0.84 0.84 init
0.57 46.86 0.27 __intel_avx_rep_memcpy
0.38 47.04 0.18 __svml_sin4_l9
0.06 47.07 0.03 __ocl_svml_h8__svml_dsin_cout_rare_internal_wrapper
0.00 47.07 0.00 7 0.00 0.00 allocGrid
0.00 47.07 0.00 3 0.00 0.32 computeResidual
0.00 47.07 0.00 3 0.00 0.00 freeGrid
0.00 47.07 0.00 2 0.00 0.00 getTimeStamp
0.00 47.07 0.00 1 0.00 44.31 CG
0.00 47.07 0.00 1 0.00 0.00 copyGrid
Result: The hotspot is axpby.
Performance with one core:
Performance CG = 170.354034 [MLUP/s]
OpenMP parallelization
Compile with OpenMP enabled:
icx -Ofast -xhost -qopenmp -std=c99 -pg -o ./mfcg ./mfcg.c
Run on one NUMA domain (18 cores in this example):
Performance CG = 646.201593 [MLUP/s]
Hypothesis: The code is limited by main memory bandwidth.
Roofline model of hot loops
Measure sustained memory bandwidth in 1 NUMA domain using bwbench:
likwid-pin -c S0:0-17 ./bwbench
OpenMP enabled, running with 18 threads
----------------------------------------------------------------------------
Function Rate(MB/s) Rate(MFlop/s) Avg time Min time Max time
Init: 55101.90 - 0.0150 0.0145 0.0163
Copy: 75431.99 - 0.0213 0.0212 0.0215
Update: 74215.65 4638.48 0.0217 0.0216 0.0219
Triad: 75167.49 6263.96 0.0321 0.0319 0.0324
Daxpy: 74068.96 6172.41 0.0326 0.0324 0.0327
STriad: 76363.42 4772.71 0.0422 0.0419 0.0424
SDaxpy: 75074.62 4692.16 0.0428 0.0426 0.0430
----------------------------------------------------------------------------
Solution Validates
Apply stencil
The kernel is a standard 5-point 2D Laplacian stencil (mfcg.c:238-246):
R[j*size_x + i] = w_c * U[c]
- w_y * (U[N] + U[S])
- w_x * (U[E] + U[W]);
- Floating-point operations per grid point
| Operation | Count |
|---|---|
| Multiplies (w_c*, w_y*, w_x*) | 3 |
| Adds (U[N]+U[S], U[E]+U[W]) | 2 |
| Subtractions (… - … - …) | 2 |
| Total FLOPs/point | 7 |
(w_x, w_y, w_c are loop-invariant scalars held in registers — no traffic.)
- Memory traffic per grid point (double = 8 B)
The 5 reads of U overlap between neighboring points. With a good cache and the layer condition holding (3 consecutive rows of U fit in cache), each U element is fetched from memory exactly once:
| Stream | Bytes/point |
|---|---|
| Load U (1×, the rest served from cache) | 8 |
| Store R | 8 |
| Write-allocate (RFO) on R | 8 |
| Total | 24 B |
- Arithmetic (computational) intensity
$$I = \frac{7\ \text{FLOP}}{24\ \text{B}} \approx 0.292\ \text{FLOP/B}$$
- Roofline performance bound
With sustained bandwidth b = 75 GB/s, this kernel is memory-bound, so:
$$P = I \times b = 0.292 \times 75 \approx \mathbf{21.9\ GFLOP/s}$$
Sensitivity to assumptions
| Scenario | Bytes/pt | Intensity (F/B) | Bound (GFLOP/s) |
|---|---|---|---|
| Layer condition + write-allocate (typical) | 24 | 0.292 | 21.9 |
| Non-temporal (streaming) stores, no RFO | 16 | 0.438 | 32.8 |
| Layer condition broken (each U row reloaded ~2–3×) | up to 56 | 0.125 | 9.4 |
- The kernel is strongly memory-bound (intensity ≈ 0.3 F/B is far left of any modern CPU's machine-balance ridge), so achievable performance is I × 75 GB/s, not the peak FLOP rate.
- Realistic expectation: ~22 GFLOP/s if the working set of 3 grid rows stays in cache.
- Optimisation options: (1) non-temporal stores on R eliminate the write-allocate (≈ +50% → ~33 GFLOP/s), and (2) spatial blocking in j to guarantee the layer condition for large grids where 3 rows exceed the last-level cache.
axpby
The kernel is a fused vector triad (mfcg.c:186-189):
R[i] = (a * X[i]) + (b * Y[i]);
a and b are scalar arguments held in registers — no traffic. This is a pure streaming kernel with no data reuse.
- FLOPs per element
| Operation | Count |
|---|---|
| Multiplies (a*X, b*Y) | 2 |
| Add | 1 |
| Total FLOPs/element | 3 |
- Memory traffic per element (double = 8 B)
No reuse — every element of each stream is touched exactly once:
| Stream | Bytes/element |
|---|---|
| Load X | 8 |
| Load Y | 8 |
| Store R | 8 |
| Write-allocate (RFO) on R | 8 |
| Total | 32 B |
- Arithmetic intensity
$$I = \frac{3\ \text{FLOP}}{32\ \text{B}} \approx 0.094\ \text{FLOP/B}$$
- Roofline performance bound
$$P = I \times b = 0.094 \times 75 \approx \mathbf{7.0\ GFLOP/s}$$
Important: in-place aliasing in the CG loop
Most axpby calls in this code alias the output with an input, e.g. axpby(x, 1.0, x, lambda, &p) and axpby(&p, 1.0, &r, …, &p) (mfcg.c:295-300). When R and one input point to the same array, only two distinct arrays are touched, and the write-allocate line is already resident from the read — so traffic drops to 24 B/element.
| Scenario | Bytes/elem | Intensity (F/B) | Bound (GFLOP/s) |
|---|---|---|---|
| 3 distinct arrays + write-allocate | 32 | 0.094 | 7.0 |
| In-place (R aliases an input), or non-temporal stores | 24 | 0.125 | 9.4 |
- Even more strongly memory-bound than applyStencil — intensity ≈ 0.09–0.125 F/B vs. ~0.29 for the stencil, because there is zero data reuse.
- Expected performance: ~7 GFLOP/s (distinct arrays) up to ~9.4 GFLOP/s for the in-place calls that dominate the CG iteration.
- The only software lever is non-temporal stores to drop the write-allocate (32 → 24 B), worth ~+33%. Otherwise it runs at exactly the bandwidth limit regardless of the CPU's peak FLOP rate.
Data volume prediction
Problem parameters
- size_x = 8000, size_y = 20000 → N = 1.6 × 10⁸ grid points (each double = 8 B; one grid = 1.28 GB)
- itermax defaults to 50 (only 2 size args given, mfcg.c:322,335) — assume the CG loop runs the full 50 iterations.
Call counts (whole program)
applyStencil — called once per computeResidual (mfcg.c:259) plus once per CG iteration:
| Site | Calls |
|---|---|
| main line 355 + 370 (computeResidual) | 2 |
| CG line 284 (computeResidual, pre-loop) | 1 |
| CG line 290 (in loop) | 50 |
| Total | 53 |
axpby — mfcg.c:353,368 and three per CG iteration (mfcg.c:295,297,300):
| Site | Calls | Aliasing |
|---|---|---|
| main line 353 + 368 | 2 | 3 distinct arrays → 32 B/elem |
| CG loop lines 295/297/300 | 3 × 50 = 150 | output aliases an input → 24 B/elem |
| Total | 152 |
▎ Note: all three in-loop axpby calls write back into one of ▎ their inputs (x, r, p), so the write-allocate line is already ▎ resident → 24 B/element, not 32.
Data volume (write-allocate cache model)
applyStencil — 24 B/point (load U + store R + RFO): $$V = 53 \times 24,\text{B} \times 1.6\times10^8 \approx \mathbf{203.5\ GB}$$
axpby:
- CG loop (aliased, 24 B): $150 \times 24 \times 1.6\times10^8 = 576\ \text{GB}$
main (non-aliased, 32 B): $2 \times 32 \times 1.6\times10^8 = 10.2\ \text{GB}$
$$V \approx \mathbf{586.2\ GB}$$
Summary
| Kernel | Calls | B/elem | Total volume |
|---|---|---|---|
| applyStencil | 53 | 24 | ≈ 204 GB |
| axpby | 152 | 24 / 32 | ≈ 586 GB |
| Combined | ≈ 790 GB |
At the sustained 75 GB/s, the data-transfer lower bound on runtime for these two kernels is:
- applyStencil: 204 / 75 ≈ 2.7 s
- axpby: 586 / 75 ≈ 7.8 s
- combined ≈ 10.5 s
This makes clear that despite applyStencil being the more "compute-looking" kernel, axpby moves ~3× more data and dominates the bandwidth budget — because it is invoked 3× per CG iteration with no data reuse.
(With non-temporal stores eliminating the write-allocate, applyStencil drops to 16 B/pt → ~136 GB and axpby's non-aliased calls to 24 B; the aliased loop calls are unchanged. The combined volume would fall to ~700 GB.)
If no aliasing is assumed, every axpby call moves 3 distinct arrays + write-allocate = 32 B/element, for all 152 calls:
$$V = 152 \times 32,\text{B} \times 1.6\times10^8 = \mathbf{778.2\ GB}$$
| Scenario | B/elem | Calls | Volume |
|---|---|---|---|
| With aliasing (in-loop calls = 24 B) | 24 / 32 | 152 | 586 GB |
| No aliasing (all 32 B) | 32 | 152 | 778 GB |
At 75 GB/s the transfer lower bound rises to 778 / 75 ≈ 10.4 s for axpby alone (vs. ~7.8 s with aliasing).
Performance profiling
Use likwid-perfctr to measure memory bandwidth. Example code is already
instrumented, only applyStencil region has to be added.
Compile with likwid:
icx -qopenmp -Ofast -xhost -std=c99 -DLIKWID_PERFMON -o ./mfcg $LIKWID_INC ./mfcg-marker.c $LIKWID_LIB -llikwid
Region axpby, Group 1: MEM_DP
+-------------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| Region Info | HWThread 0 | HWThread 1 | HWThread 2 | HWThread 3 | HWThread 4 | HWThread 5 | HWThread 6 | HWThread 7 | HWThread 8 | HWThread 9 | HWThread 10 | HWThread 11 | HWThread 12 | HWThread 13 | HWThread 14 | HWThread 15 | HWThread 16 | HWThread 17 |
+-------------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| RDTSC Runtime [s] | 8.077774 | 8.107685 | 8.107348 | 8.107557 | 8.107642 | 8.107953 | 8.107596 | 8.107856 | 8.107682 | 8.107646 | 8.107734 | 8.107312 | 8.107751 | 8.107581 | 8.107589 | 8.107107 | 8.108067 | 8.107872 |
| call count | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 | 155 |
+-------------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
+----------------------------------------+------------+-----------+------------+-----------+
| Metric | Sum | Min | Max | Avg |
+----------------------------------------+------------+-----------+------------+-----------+
| Runtime (RDTSC) [s] STAT | 145.9079 | 8.0778 | 8.1081 | 8.1060 |
| Runtime unhalted [s] STAT | 189.0819 | 10.4393 | 10.6106 | 10.5045 |
| Clock [MHz] STAT | 56521.6280 | 3121.7440 | 3169.0521 | 3140.0904 |
| CPI STAT | 201.9758 | 10.5501 | 12.1314 | 11.2209 |
| Energy [J] STAT | 1394.5000 | 0 | 1394.5000 | 77.4722 |
| Power [W] STAT | 172.6342 | 0 | 172.6342 | 9.5908 |
| Energy DRAM [J] STAT | 115.5714 | 0 | 115.5714 | 6.4206 |
| Power DRAM [W] STAT | 14.3073 | 0 | 14.3073 | 0.7948 |
| DP [MFLOP/s] STAT | 9178.4173 | 509.7808 | 511.6926 | 509.9121 |
| AVX DP [MFLOP/s] STAT | 9178.4059 | 509.7802 | 511.6920 | 509.9114 |
| Packed [MUOPS/s] STAT | 2294.6017 | 127.4450 | 127.9230 | 127.4779 |
| Scalar [MUOPS/s] STAT | 0.0108 | 0.0006 | 0.0006 | 0.0006 |
| Memory read bandwidth [MBytes/s] STAT | 49289.6386 | 0 | 49289.6386 | 2738.3133 |
| Memory read data volume [GBytes] STAT | 398.1506 | 0 | 398.1506 | 22.1195 |
| Memory write bandwidth [MBytes/s] STAT | 24430.4265 | 0 | 24430.4265 | 1357.2459 |
| Memory write data volume [GBytes] STAT | 197.3435 | 0 | 197.3435 | 10.9635 |
| Memory bandwidth [MBytes/s] STAT | 73720.0651 | 0 | 73720.0651 | 4095.5592 |
| Memory data volume [GBytes] STAT | 595.4940 | 0 | 595.4940 | 33.0830 |
| Operational intensity [FLOP/Byte] STAT | 0.1242 | 0.0069 | 0.0069 | 0.0069 |
| Vectorization ratio [%] STAT | 1799.9910 | 99.9995 | 99.9995 | 99.9995 |
+----------------------------------------+------------+-----------+------------+-----------+
Region stencil, Group 1: MEM_DP
+-------------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| Region Info | HWThread 0 | HWThread 1 | HWThread 2 | HWThread 3 | HWThread 4 | HWThread 5 | HWThread 6 | HWThread 7 | HWThread 8 | HWThread 9 | HWThread 10 | HWThread 11 | HWThread 12 | HWThread 13 | HWThread 14 | HWThread 15 | HWThread 16 | HWThread 17 |
+-------------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| RDTSC Runtime [s] | 2.381165 | 2.390895 | 2.390732 | 2.390819 | 2.390760 | 2.391043 | 2.390815 | 2.390833 | 2.390835 | 2.390768 | 2.390789 | 2.390812 | 2.390906 | 2.390724 | 2.390896 | 2.390586 | 2.391266 | 2.391017 |
| call count | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 | 53 |
+-------------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
+----------------------------------------+------------+-----------+------------+-----------+
| Metric | Sum | Min | Max | Avg |
+----------------------------------------+------------+-----------+------------+-----------+
| Runtime (RDTSC) [s] STAT | 43.0256 | 2.3812 | 2.3913 | 2.3903 |
| Runtime unhalted [s] STAT | 54.6628 | 3.0105 | 3.0736 | 3.0368 |
| Clock [MHz] STAT | 56282.3196 | 3101.1055 | 3173.6957 | 3126.7955 |
| CPI STAT | 28.8189 | 1.5729 | 1.6380 | 1.6010 |
| Energy [J] STAT | 443.4335 | 0 | 443.4335 | 24.6352 |
| Power [W] STAT | 186.2254 | 0 | 186.2254 | 10.3459 |
| Energy DRAM [J] STAT | 34.1181 | 0 | 34.1181 | 1.8954 |
| Power DRAM [W] STAT | 14.3283 | 0 | 14.3283 | 0.7960 |
| DP [MFLOP/s] STAT | 24824.8929 | 1378.6107 | 1384.4589 | 1379.1607 |
| AVX DP [MFLOP/s] STAT | 0 | 0 | 0 | 0 |
| Packed [MUOPS/s] STAT | 0 | 0 | 0 | 0 |
| Scalar [MUOPS/s] STAT | 24824.8929 | 1378.6107 | 1384.4589 | 1379.1607 |
| Memory read bandwidth [MBytes/s] STAT | 46469.0200 | 0 | 46469.0200 | 2581.6122 |
| Memory read data volume [GBytes] STAT | 110.6504 | 0 | 110.6504 | 6.1472 |
| Memory write bandwidth [MBytes/s] STAT | 29213.3817 | 0 | 29213.3817 | 1622.9657 |
| Memory write data volume [GBytes] STAT | 69.5619 | 0 | 69.5619 | 3.8645 |
| Memory bandwidth [MBytes/s] STAT | 75682.4018 | 0 | 75682.4018 | 4204.5779 |
| Memory data volume [GBytes] STAT | 180.2123 | 0 | 180.2123 | 10.0118 |
| Operational intensity [FLOP/Byte] STAT | 0.3294 | 0.0183 | 0.0183 | 0.0183 |
| Vectorization ratio [%] STAT | 0 | 0 | 0 | 0 |
+----------------------------------------+------------+-----------+------------+-----------+
The stencil fits the prediction assuming active NUMA balancing. For axpby it is slightly too good if only taking the kernel into account.
Scalability across NUMA domains
Scalability across NUMA domains:
- 1 domain: 652 MLUPS
- 2 domain: 1190 MLUPS (x1.8)
- 3 domain: 1561 MLUPS (x2.3)
- 4 domain: 1810 MLUPS (x2.8)
Optimizations
Two main targets:
- Correct first touch placement (specifically copyGrid)
- Fuse loops to lower data volume
- axpyDot: Replaces axpby(r,1,r,-lambda,v) + dotProduct(r,r) -> saves a full re-read of r.
- stencilDot: Replaces applyStencil(v,p) + dotProduct(v,p) -> saves a re-read of v and p.
Result with correct first touch policy (and fused kernels)
- 1 domain: 909 MLUPS
- 2 domain: 1811 MLUPS (x1.99)
- 3 domain: 2697 MLUPS (x2.96)
- 4 domain: 3590 MLUPS (x3.95)
Total data volume:
- Before: 989 GB
- After: 757 GB (-24%)