Make the cross-language benchmark measure the languages
The first version of the offline benchmarks had Julia about ten times faster
than MATLAB and Python, and almost none of that was a property of Julia. Three
things were wrong, and all three are the kind of thing that makes a
cross-language comparison measure the author instead.
The kernel matrix was assembled in a different amount of work per language.
Both boundary integral scripts formed their coordinate differences by
multiplying against ones(1,n), which does rank-one matrix multiplies and
allocates four n-by-n temporaries. At n = 1896 that cost real MATLAB 62 ms
against 19 ms for implicit expansion and 9 ms for a scalar double loop; under
numbl the same three shapes are 117, 99 and 50 ms, and in Julia the loop beats
broadcasting by more than tenfold. So the loop is now what ch10/bie.m and
ch11/pacman.m use, which speeds up the browser too: Chapter 10's Nystrom solve
at n = 512 drops from 13 ms to 6 ms, and Chapter 11's deepest panels from 53
to 25. The old idiom was never required by numbl, which supports implicit
expansion and is itself faster with it.
NumPy cannot run that loop at any acceptable speed, so the Python port keeps
whole-array expressions and now performs them in place, reusing two buffers
instead of allocating five; and it asks for the same linear algebra as the
others, a pivoted QR through SciPy's gelsy driver rather than the SVD behind
numpy.linalg.lstsq. Julia's assembly gains @inbounds and column-major order.
Holding the code shape fixed across languages would have been the unfair
choice, not the fair one, and the scripts now say so at the top.
Timing was measured differently in each environment: Julia took the best of
two runs while the others took one whenever a solve exceeded 5 ms. All four
environments, the page included, now keep the fastest of several runs, at
least three whenever a solve is cheap enough to allow it, which matters
because a single measurement of a twenty-millisecond solve on this hardware
varies by a factor of two with whatever else the machine is doing. Chapter 10
adopts the same rule, since its prose quotes its timings; its table is
updated, and every assertion still passes.
What is left is a comparison worth reading. On a 10-core i7-1355U the four
environments now land within a factor of two of each other, MATLAB and Julia
within about ten per cent up to a 1300-square system, with Python's remaining
1.3 to 1.7 attributable to the assembly stage and MATLAB's lead at the largest
sizes to its BLAS taking ten threads where Julia's default takes six. Each
label now records that thread count, and the chapter has a new section on the
three disciplines the comparison rests on.
13 changed files+348−99
.gitignoremodified+1−0View file
@@ -4,3 +4,4 @@ dist
44 *.local
55 tmp
66 bench-pacman-*.json
7+__pycache__
CLAUDE.mdmodified+39−4View file
@@ -111,6 +111,24 @@ form allocates several temporaries per sweep. The kernels in
111111 `ch08/iterative.m` are written as loops for that reason, which also happens to
112112 make them readable.
113113
114+The same holds for building a dense kernel matrix, and there the measurement is
115+worth keeping. Three shapes for the Nyström matrix, at n = 1896:
116+
117+| shape | numbl | MATLAB | Julia |
118+|---|---|---|---|
119+| scalar double loop | 50 ms | 9 ms | 7 ms |
120+| implicit expansion, `x - x'` | 99 ms | 13 ms | 82 ms |
121+| outer products, `x * ones(1,n)` | 117 ms | 62 ms | — |
122+
123+So the loop wins everywhere the tutorial runs, and `ch10/bie.m` and
124+`ch11/pacman.m` both use it. The third row is what those files used to
125+contain; the idiom is not required by numbl, which supports implicit expansion
126+and is itself faster with it, and it was costing real MATLAB a factor of
127+seven in that stage. If you find it anywhere else, it is worth an experiment.
128+NumPy is the exception that proves the rule: it cannot run the loop at all, so
129+the Python port in `bench/` assembles with whole-array expressions in place,
130+and the chapter's prose makes that asymmetry part of the point.
131+
114132 ### The browser is slower than the CLI
115133
116134 The CLI uses a native LAPACK addon; the browser falls back to a TypeScript
@@ -123,10 +141,27 @@ were what a reader will see.
123141
124142 Timings come from `tic`/`toc` **inside** the MATLAB layer, around the solve only, so
125143 they exclude session boot, JSON encoding and rendering. Each method gets a warm-up
126-call before its first timed solve, because the first call also pays for compilation,
127-and any solve under 5 ms is repeated a few times with the minimum kept. Below about
128-0.1 ms the measurement is dominated by fixed overhead and the work-precision curves
129-flatten out on the left; that is the clock, not the method.
144+call before its first timed solve, because the first call also pays for compilation.
145+Below about 0.1 ms the measurement is dominated by fixed overhead and the
146+work-precision curves flatten out on the left; that is the clock, not the method.
147+
148+Chapters 10 and 11, whose numbers are compared against other environments or quoted
149+in the prose, keep the **fastest of several runs**: at least three whenever a solve
150+is cheap enough to allow it, many more when it is very cheap, one when it is
151+expensive, with the loop bounded to about 0.6 s per point.
152+
153+```matlab
154+while rep == 0 || (rep < 50 && tot < 0.6 && (rep < 3 || tot < 0.1))
155+```
156+
157+That rule is not cosmetic. A single measurement of a twenty-millisecond solve on
158+this hardware varies by a factor of two with whatever else the machine is doing,
159+which is enough to reorder neighbouring points on a work-precision plot and, in
160+the offline benchmarks, to make one language look several times faster than
161+another. The three scripts under `bench/` implement the same rule, and any new
162+cross-environment comparison must too. The earlier chapters still use the older
163+rule (repeat only when under 5 ms), which is fine because nothing outside those
164+chapters is compared against their numbers.
130165
131166 ## Testing
132167
README.mdmodified+5−3View file
@@ -35,7 +35,7 @@ error of the method.
3535 ## The chapters
3636
3737 **Introduction** — Laplace's and Poisson's equations, the Dirichlet, Neumann and
38-Robin boundary conditions, what accuracy-against-cost means here, and the four
38+Robin boundary conditions, what accuracy-against-cost means here, and the five
3939 domains the chapters work on.
4040
4141 **One dimension**
@@ -107,8 +107,10 @@ npm run test:site # the deployed site, in headless Chrome
107107 ```
108108
109109 The offline benchmarks for Chapter 11 live in `bench/`, with their own README;
110-they need only whichever of MATLAB, Python + NumPy, or Julia is at hand, and each
111-writes a JSON file that the benchmark chapter's plot loads.
110+they need only whichever of MATLAB, Python (NumPy and SciPy), or Julia is at
111+hand, and each writes a JSON file that the benchmark chapter's plot loads. That
112+README also records the three rules that keep a cross-language comparison from
113+measuring the author instead of the languages.
112114
113115 `test:matlab` needs a numbl clone; set `NUMBL_DIR` if it is not at `~/src/numbl`,
114116 and build its CLI (`npm run build:cli` in that checkout) after updating it — the
bench/README.mdmodified+37−5View file
@@ -3,13 +3,12 @@
33 The benchmark chapter of the tutorial runs its four methods in the browser
44 through numbl. The scripts here run the two best of those methods, the
55 harmonic basis and the graded-panel double layer, outside the browser, on the
6-same Pac-Man problem, the same error lattice, and the same timing protocol
7-(solve only, warmed up, minimum over repeats when a solve is under 5 ms).
6+same Pac-Man problem and the same error lattice.
87
98 Run whichever of these your machine has, from the repository root:
109
1110 ```
12-python3 bench/python/bench_pacman.py # needs NumPy
11+python3 bench/python/bench_pacman.py # needs NumPy and SciPy
1312 julia bench/julia/bench_pacman.jl # standard library only
1413 matlab -batch "run('bench/matlab/bench_pacman.m')"
1514 ```
@@ -26,8 +25,41 @@ results then share a clock and can honestly share an axis.
2625 The MATLAB harness runs the chapter's own driver and editable script,
2726 unchanged, in real MATLAB. The Python and Julia scripts are ports; the check
2827 that they are faithful is that their error columns land on top of the
29-browser's curves, since the error depends only on the mathematics. All three
30-were verified to reproduce the browser's errors to plotting accuracy.
28+browser's curves, since the error depends only on the mathematics. On the
29+double layer the three environments agree bit for bit; on the basis method
30+they agree to within a factor of two at the $10^{-15}$ roundoff floor, where
31+the order of the QR pivoting is all that separates them.
32+
33+## Comparing languages, not authors
34+
35+A cross-language benchmark measures how carefully each port was written unless
36+that is deliberately controlled for, so three rules are held to here. The
37+chapter's prose discusses them; in the code they look like this.
38+
39+**Each port uses the fastest shape its own language offers, and they differ.**
40+The kernel matrix is assembled by a scalar double loop in the MATLAB script
41+(hence in numbl too) and in Julia, because that is about twice as fast as
42+implicit expansion in MATLAB and numbl and much faster than broadcasting in
43+Julia. NumPy cannot run that loop at any acceptable speed, so the Python port
44+assembles with whole-array expressions, in place: the obvious version, one
45+broadcast statement per line, needs five n-by-n temporaries and is about 1.6
46+times slower for that reason alone. Insisting on one shape everywhere would
47+make the benchmark less fair, not more.
48+
49+**The same linear algebra is requested everywhere.** `numpy.linalg.lstsq` is
50+an SVD, while MATLAB's and Julia's backslash are pivoted QR on a rectangular
51+system, so the Python port calls SciPy's QR driver (`gelsy`) instead. For the
52+square systems all three call the same LAPACK and land within ten per cent.
53+
54+**One timing protocol, everywhere, including in the page.** Keep the fastest
55+of several runs: at least three whenever a solve is cheap enough to allow it,
56+many more when it is very cheap, one when it is expensive. A single
57+measurement of a fifty-millisecond solve varies by a factor of two with
58+whatever else the machine is doing.
59+
60+Each label records the BLAS thread count, since that is a real difference
61+between environments and otherwise an invisible one: on the machine these
62+were developed on, MATLAB takes ten threads, NumPy twelve, and Julia six.
3163
3264 ## The file format
3365
bench/julia/bench_pacman.jlmodified+45−20View file
@@ -11,6 +11,14 @@
1111 # Load that file from the benchmark chapter's "load offline results" button.
1212 # Because it ran on the same machine as your browser, the timings share a
1313 # clock and can share an axis. Uses only the standard library.
14+#
15+# Each port is written the way its own language wants to be written, since the
16+# comparison is meant to be between languages and not between the author's
17+# habits in each. Here that means scalar loops with @inbounds for the matrix
18+# assembly, which is both the fastest and the most readable shape available;
19+# the same loop in NumPy would be hopeless, which is why the Python port
20+# assembles with whole-array expressions instead. Where the work is one
21+# LAPACK call, all three ports do the same thing and should agree.
1422
1523 using LinearAlgebra
1624 using Dates
@@ -166,26 +174,34 @@ end
166174 function bie_solve(d, px, py)
167175 qx, qy, qnx, qny, qw, qkap = pac_panels(d)
168176 n = length(qx)
169- A = zeros(n, n)
170- for j in 1:n, i in 1:n
171- if i == j
172- A[i, j] = -0.5 - qkap[j] / (4pi) * qw[j]
173- else
174- dx = qx[i] - qx[j]
175- dy = qy[i] - qy[j]
176- A[i, j] = (qnx[j] * dx + qny[j] * dy) / (2pi * (dx^2 + dy^2)) * qw[j]
177+ # column by column, so the writes run down memory, and with the bounds
178+ # checks elided: every index here is a loop variable over 1:n
179+ A = Matrix{Float64}(undef, n, n)
180+ c = 1 / (2pi)
181+ @inbounds for j in 1:n
182+ xj = qx[j]; yj = qy[j]
183+ nxj = qnx[j]; nyj = qny[j]; wj = qw[j]
184+ for i in 1:n
185+ dx = qx[i] - xj
186+ dy = qy[i] - yj
187+ A[i, j] = c * (nxj * dx + nyj * dy) / (dx * dx + dy * dy) * wj
177188 end
189+ # the smooth limit of the kernel on the diagonal, and the jump relation
190+ A[j, j] = -0.5 - qkap[j] * wj / (4pi)
178191 end
192+
179193 sigma = A \ exact(qx, qy)
194+
180195 v = zeros(length(px))
181- for k in eachindex(px)
196+ @inbounds for k in eachindex(px)
197+ xk = px[k]; yk = py[k]
182198 s = 0.0
183199 for j in 1:n
184- dx = px[k] - qx[j]
185- dy = py[k] - qy[j]
186- s += (qnx[j] * dx + qny[j] * dy) / (2pi * (dx^2 + dy^2)) * qw[j] * sigma[j]
200+ dx = xk - qx[j]
201+ dy = yk - qy[j]
202+ s += (qnx[j] * dx + qny[j] * dy) / (dx * dx + dy * dy) * qw[j] * sigma[j]
187203 end
188- v[k] = s
204+ v[k] = c * s
189205 end
190206 return v
191207 end
@@ -193,18 +209,27 @@ end
193209 # ------------------------------------------------------------- the sweep ----
194210
195211 function sweep(solve, ns, px, py, ve)
212+ # Warm up, then time each resolution the way the chapter driver does.
213+ # Keep the fastest of several runs. Three runs at least whenever a solve
214+ # is cheap enough to allow it, many more when it is very cheap, and one
215+ # when it is expensive: the budget below caps the wall clock per point.
216+ # The minimum over runs is what makes these numbers comparable at all,
217+ # here and in the offline scripts, which apply the same rule. A single
218+ # measurement of a fifty-millisecond solve reports whatever else the
219+ # machine happened to be doing at the time.
196220 E = Float64[]
197221 T = Float64[]
222+ local v
198223 solve(ns[1], px[1:1], py[1:1]) # warm-up, which in Julia is compilation
199224 for n in ns
200- el = @elapsed v = solve(n, px, py)
201- # a second timed run always: stray compilation and GC pauses
202- # otherwise land in single measurements
203- el = min(el, @elapsed v = solve(n, px, py))
225+ el = Inf
226+ tot = 0.0
204227 rep = 0
205- while el < 0.005 && rep < 4
228+ while rep == 0 || (rep < 50 && tot < 0.6 && (rep < 3 || tot < 0.1))
229+ dt = @elapsed v = solve(n, px, py)
230+ el = min(el, dt)
231+ tot += dt
206232 rep += 1
207- el = min(el, @elapsed v = solve(n, px, py))
208233 end
209234 push!(E, maximum(abs.(v .- ve)))
210235 push!(T, el * 1000)
@@ -234,7 +259,7 @@ function main()
234259 push!(entries,
235260 """{"name":"$(name)","slot":"$(slot)","ns":$(json_arr(ns)),"E":$(json_arr(E)),"T":$(json_arr(T))}""")
236261 end
237- label = "Julia $(VERSION.major).$(VERSION.minor)"
262+ label = "Julia $(VERSION.major).$(VERSION.minor), $(BLAS.get_num_threads()) threads"
238263 open("bench-pacman-julia.json", "w") do f
239264 write(f,
240265 """{"schema":"laplace-dirichlet-bench-1","label":"$(label)","kind":"$(KIND)","date":"$(Dates.today())","methods":[$(join(entries, ","))]}""")
bench/matlab/bench_pacman.mmodified+9−1View file
@@ -7,6 +7,14 @@
77 % double layer), with a deeper panel ladder than the page uses, since a native
88 % LAPACK has no trouble with the larger dense systems.
99 %
10+% Because this runs the chapter's own script, whatever that script does badly
11+% it does badly here too, which is a virtue for comparing MATLAB against numbl
12+% and a hazard for comparing MATLAB against the Python and Julia ports. The
13+% shape the script uses, a scalar double loop for the kernel matrix, happens to
14+% be the fastest of the three shapes available in both MATLAB and numbl, so
15+% nothing is being given away; if that ever stops being true, the honest fix is
16+% to change the chapter, not to keep a faster copy here.
17+%
1018 % matlab -batch "run('bench/matlab/bench_pacman.m')"
1119 %
1220 % from the repository root writes bench-pacman-matlab.json into the
@@ -76,7 +84,7 @@ end
7684
7785 out = struct();
7886 out.schema = 'laplace-dirichlet-bench-1';
79-out.label = ['MATLAB ' version('-release')];
87+out.label = sprintf('MATLAB %s, %d threads', version('-release'), maxNumCompThreads);
8088 out.kind = kind;
8189 out.date = char(datetime('now', 'Format', 'yyyy-MM-dd'));
8290 out.methods = methods;
bench/python/bench_pacman.pymodified+94−24View file
@@ -1,5 +1,5 @@
11 #!/usr/bin/env python3
2-"""Offline benchmark of the Pac-Man problem in Python (NumPy).
2+"""Offline benchmark of the Pac-Man problem in Python (NumPy and SciPy).
33
44 A port of the two best methods from the in-browser benchmark, the harmonic
55 basis and the graded-panel double layer, on the same problem, the same error
@@ -11,6 +11,24 @@ from anywhere; it writes bench-pacman-python.json into the current directory.
1111 Load that file from the benchmark chapter's "load offline results" button.
1212 Because it ran on the same machine as your browser, the timings share a clock
1313 and can share an axis.
14+
15+The point of the comparison is the language and its libraries, not the
16+author's fluency in each, so each port is written the way its own language
17+wants to be written. For NumPy that means whole-array expressions, with two
18+details that matter more than they look:
19+
20+ * the kernel matrix is assembled in place. Written as the obvious sequence
21+ of broadcast expressions it allocates five n-by-n temporaries and spends
22+ more time on memory traffic than the dense solve itself spends on
23+ arithmetic; reusing two buffers cuts that stage by nearly half. Writing
24+ it as a scalar double loop, which is the fastest shape in MATLAB and
25+ Julia, would be catastrophic here, and that asymmetry is the honest
26+ reason the ports do not look alike.
27+
28+ * the least-squares fit goes through a QR driver. numpy.linalg.lstsq is an
29+ SVD (LAPACK gelsd), while MATLAB's backslash and Julia's backslash are
30+ both pivoted QR, so using lstsq would have Python solving a harder
31+ problem than the others and calling it the same benchmark.
1432 """
1533
1634 import json
@@ -20,6 +38,7 @@ import time
2038 from datetime import date
2139
2240 import numpy as np
41+from scipy.linalg import lstsq as qr_lstsq
2342
2443 KIND = "both"
2544 Z0 = 1.6 + 1.3j
@@ -58,13 +77,22 @@ def lattice():
5877
5978
6079 def pac_basis(x, y, M, K):
80+ """Re z^m and Im z^m up to degree M, then K corner functions.
81+
82+ Written into a preallocated Fortran-ordered array: LAPACK wants columns
83+ contiguous, and building a list of columns for column_stack would copy the
84+ whole matrix once more than it has to.
85+ """
6186 z = x + 1j * y
62- cols = [np.ones_like(x)]
87+ A = np.empty((x.size, 2 * M + 1 + K), order="F")
88+ A[:, 0] = 1.0
89+ col = 1
6390 p = z.copy()
6491 for _ in range(M):
65- cols.append(p.real)
66- cols.append(p.imag)
67- p = p * z
92+ A[:, col] = p.real
93+ A[:, col + 1] = p.imag
94+ col += 2
95+ p *= z
6896 th = np.arctan2(y, x)
6997 th = np.where(th < -1e-13, th + 2 * np.pi, th)
7098 r = np.abs(z)
@@ -73,8 +101,9 @@ def pac_basis(x, y, M, K):
73101 k += 1
74102 if k % 3 != 0:
75103 q += 1
76- cols.append(r ** (2.0 * k / 3.0) * np.sin(2.0 * k * th / 3.0))
77- return np.column_stack(cols)
104+ A[:, col] = r ** (2.0 * k / 3.0) * np.sin(2.0 * k * th / 3.0)
105+ col += 1
106+ return A
78107
79108
80109 def pac_bpoints(ntot):
@@ -92,8 +121,10 @@ def mps_solve(n, px, py):
92121 nb = 3 * n + 1
93122 bx, by = pac_bpoints(3 * nb)
94123 A = pac_basis(bx, by, n, n)
95- # QR-based least squares, the closest match to MATLAB's backslash
96- c, *_ = np.linalg.lstsq(A, exact(bx, by), rcond=None)
124+ # gelsy is pivoted QR: what MATLAB's and Julia's backslash do on a
125+ # rectangular system. The default lstsq driver is an SVD and does more
126+ # work than either.
127+ c = qr_lstsq(A, exact(bx, by), lapack_driver="gelsy")[0]
97128 return pac_basis(px, py, n, n) @ c
98129
99130
@@ -141,14 +172,27 @@ def pac_panels(d):
141172 def bie_solve(d, px, py):
142173 qx, qy, qnx, qny, qw, qkap = pac_panels(d)
143174 n = qx.size
144- DX = qx[:, None] - qx[None, :]
145- DY = qy[:, None] - qy[None, :]
146- R2 = DX**2 + DY**2
147- np.fill_diagonal(R2, 1.0)
148- K = (qnx[None, :] * DX + qny[None, :] * DY) / (2 * np.pi * R2)
149- np.fill_diagonal(K, -qkap / (4 * np.pi))
150- A = -0.5 * np.eye(n) + K * qw[None, :]
175+
176+ # A is built in two buffers and finished in place: every operation below
177+ # writes into memory that already exists. The obvious version, one
178+ # broadcast expression per line, needs five n-by-n arrays and is about
179+ # 1.6 times slower for that reason alone.
180+ A = qx[:, None] - qx[None, :] # dx
181+ B = qy[:, None] - qy[None, :] # dy
182+ R2 = A * A
183+ R2 += B * B
184+ np.fill_diagonal(R2, 1.0) # keep 0/0 off the diagonal
185+ A *= qnx[None, :]
186+ B *= qny[None, :]
187+ A += B
188+ A /= R2
189+ A *= qw[None, :] / (2 * np.pi)
190+ # the diagonal is the smooth limit of the kernel, plus the jump relation
191+ diag = np.arange(n)
192+ A[diag, diag] = -0.5 - qkap * qw / (4 * np.pi)
193+
151194 sigma = np.linalg.solve(A, exact(qx, qy))
195+
152196 EX = px[:, None] - qx[None, :]
153197 EY = py[:, None] - qy[None, :]
154198 KE = (qnx[None, :] * EX + qny[None, :] * EY) / (2 * np.pi * (EX**2 + EY**2))
@@ -159,23 +203,45 @@ def bie_solve(d, px, py):
159203
160204
161205 def sweep(solve, ns, px, py, ve):
206+ """Warm up, then time each resolution the way the chapter driver does.
207+
208+ Keep the fastest of several runs. Three runs at least whenever a solve
209+ is cheap enough to allow it, many more when it is very cheap, and one
210+ when it is expensive: the budget below caps the wall clock per point.
211+ The minimum over runs is what makes these numbers comparable at all,
212+ here and in the offline scripts, which apply the same rule. A single
213+ measurement of a fifty-millisecond solve reports whatever else the
214+ machine happened to be doing at the time.
215+ """
162216 E, T = [], []
163217 solve(ns[0], px[:1], py[:1]) # warm-up
164218 for n in ns:
165- t0 = time.perf_counter()
166- v = solve(n, px, py)
167- el = time.perf_counter() - t0
168- rep = 0
169- while el < 0.005 and rep < 4:
170- rep += 1
219+ el, tot, rep = float("inf"), 0.0, 0
220+ while rep == 0 or (rep < 50 and tot < 0.6 and (rep < 3 or tot < 0.1)):
171221 t0 = time.perf_counter()
172222 v = solve(n, px, py)
173- el = min(el, time.perf_counter() - t0)
223+ dt = time.perf_counter() - t0
224+ el = min(el, dt)
225+ tot += dt
226+ rep += 1
174227 E.append(float(np.max(np.abs(v - ve))))
175228 T.append(el * 1000.0)
176229 return E, T
177230
178231
232+def threads():
233+ """How many threads the BLAS behind NumPy will use, if we can find out."""
234+ try:
235+ from threadpoolctl import threadpool_info
236+
237+ n = max((p["num_threads"] for p in threadpool_info()), default=None)
238+ if n:
239+ return n
240+ except ImportError:
241+ pass
242+ return None
243+
244+
179245 def main():
180246 px, py = lattice()
181247 ve = exact(px, py)
@@ -190,9 +256,13 @@ def main():
190256 for n, e, t in zip(ns, E, T):
191257 print(f" n={n:3d} err={e:.3e} {t:8.2f} ms")
192258
259+ nt = threads()
260+ label = f"Python {sys.version_info.major}.{sys.version_info.minor} + NumPy {np.__version__}"
261+ if nt:
262+ label += f", {nt} threads"
193263 out = {
194264 "schema": "laplace-dirichlet-bench-1",
195- "label": f"Python {sys.version_info.major}.{sys.version_info.minor} + NumPy {np.__version__}",
265+ "label": label,
196266 "kind": KIND,
197267 "date": date.today().isoformat(),
198268 "machine": platform.machine(),
src/chapters/ch10/bie.mmodified+35−13View file
@@ -58,22 +58,44 @@ nx = yp ./ sp; % outward normal, counterclockwise curve
5858 ny = -xp ./ sp;
5959 kappa = (xp .* ypp - yp .* xpp) ./ sp.^3;
6060
61-% the Nystrom matrix: A = -I/2 + h * K * diag(speed)
62-DX = x * ones(1, n) - ones(n, 1) * x';
63-DY = y * ones(1, n) - ones(n, 1) * y';
64-R2 = DX.^2 + DY.^2;
65-K = ((ones(n, 1) * nx') .* DX + (ones(n, 1) * ny') .* DY) ./ (2*pi*R2);
66-K(1:n+1:end) = -kappa / (4*pi); % the diagonal limit, not a singularity
67-A = -0.5*eye(n) + h * K .* (ones(n, 1) * sp');
61+% The Nystrom matrix, A = -I/2 + h * K * diag(speed), entry by entry. A
62+% double loop is about twice as fast here as the same expression written with
63+% implicit expansion, and several times faster than forming the coordinate
64+% differences with outer products against ones(1,n); it also writes each
65+% entry once instead of passing over four n-by-n temporaries.
66+A = zeros(n, n);
67+c = h / (2*pi);
68+for j = 1:n
69+ xj = x(j);
70+ yj = y(j);
71+ nxj = nx(j);
72+ nyj = ny(j);
73+ sj = sp(j);
74+ for i = 1:n
75+ dx = x(i) - xj;
76+ dy = y(i) - yj;
77+ A(i, j) = c * (nxj*dx + nyj*dy) / (dx*dx + dy*dy) * sj;
78+ end
79+ % the diagonal limit, not a singularity, plus the jump relation's -1/2
80+ A(j, j) = -0.5 - c * kappa(j) * sj / 2;
81+end
6882
6983 sigma = A \ lp_bexact(x, y);
7084
71-% evaluate the layer potential at the target points
72-np = numel(px);
73-EX = px(:) * ones(1, n) - ones(np, 1) * x';
74-EY = py(:) * ones(1, n) - ones(np, 1) * y';
75-KE = ((ones(np, 1) * nx') .* EX + (ones(np, 1) * ny') .* EY) ./ (2*pi*(EX.^2 + EY.^2));
76-v = (h * KE .* (ones(np, 1) * sp')) * sigma;
85+% evaluate the layer potential at the target points, by the same loop
86+nt = numel(px);
87+v = zeros(nt, 1);
88+for k = 1:nt
89+ xk = px(k);
90+ yk = py(k);
91+ s = 0;
92+ for j = 1:n
93+ dx = xk - x(j);
94+ dy = yk - y(j);
95+ s = s + (nx(j)*dx + ny(j)*dy) / (dx*dx + dy*dy) * sp(j) * sigma(j);
96+ end
97+ v(k) = c * s;
98+end
7799 end
78100
79101 function v = mfs(n, a, nlobe, par, px, py)
src/chapters/ch10/driver.mmodified+13−6View file
@@ -82,15 +82,22 @@ else
8282 vw = solve_bie(k, 16, a, nlobe, p.par, px(1), py(1)); % warm-up
8383 for j = 1:nn
8484 n = ns(j);
85- t0 = tic;
86- v = solve_bie(k, n, a, nlobe, p.par, px, py);
87- el = toc(t0);
85+ % Keep the fastest of several runs: three at least whenever a
86+ % solve is cheap enough to allow it, many more when it is very
87+ % cheap, one when it is expensive. A single measurement of a
88+ % twenty-millisecond solve varies by a factor of two with whatever
89+ % else the machine is doing, and the numbers here are quoted in the
90+ % prose. Chapter 11 times itself the same way.
91+ el = Inf;
92+ tot = 0;
8893 rep = 0;
89- while el < 0.005 && rep < 4
90- rep = rep + 1;
94+ while rep == 0 || (rep < 50 && tot < 0.6 && (rep < 3 || tot < 0.1))
9195 t0 = tic;
9296 v = solve_bie(k, n, a, nlobe, p.par, px, py);
93- el = min(el, toc(t0));
97+ dt = toc(t0);
98+ el = min(el, dt);
99+ tot = tot + dt;
100+ rep = rep + 1;
94101 end
95102 q = solve_bie(k, n, a, nlobe, p.par, qx, qy);
96103 E(k, j) = max(abs(v - ve));
src/chapters/ch10/text.mdmodified+2−2View file
@@ -38,8 +38,8 @@ On the default geometry, measured at interior points out to three quarters of th
3838
3939 | method | error | boundary unknowns | time |
4040 |---|---|---|---|
41-| double-layer Nyström | $5 \times 10^{-16}$ | 512 | 13 ms |
42-| fundamental solutions | $2 \times 10^{-14}$ | 128 | 12 ms |
41+| double-layer Nyström | $2 \times 10^{-15}$ | 512 | 6 ms |
42+| fundamental solutions | $2 \times 10^{-14}$ | 128 | 3 ms |
4343 | five-point on the square (Ch. 6) | $1 \times 10^{-8}$ | $10^{6}$ | 440 ms |
4444
4545 The first row is the point of the whole tutorial in one line: full double precision, on a curved domain, from a few hundred unknowns and a few milliseconds, because the formulation put the unknown where the information is and the quadrature matched the smoothness of the integrand.
src/chapters/ch11/driver.mmodified+14−6View file
@@ -116,15 +116,23 @@ else
116116 if ~isfinite(n) || n <= 0
117117 continue;
118118 end
119- t0 = tic;
120- v = solve_pac(k, n, g, px, py);
121- el = toc(t0);
119+ % Keep the fastest of several runs. Three runs at least whenever a solve
120+ % is cheap enough to allow it, many more when it is very cheap, and one
121+ % when it is expensive: the budget below caps the wall clock per point.
122+ % The minimum over runs is what makes these numbers comparable at all,
123+ % here and in the offline scripts, which apply the same rule. A single
124+ % measurement of a fifty-millisecond solve reports whatever else the
125+ % machine happened to be doing at the time.
126+ el = Inf;
127+ tot = 0;
122128 rep = 0;
123- while el < 0.005 && rep < 4
124- rep = rep + 1;
129+ while rep == 0 || (rep < 50 && tot < 0.6 && (rep < 3 || tot < 0.1))
125130 t0 = tic;
126131 v = solve_pac(k, n, g, px, py);
127- el = min(el, toc(t0));
132+ dt = toc(t0);
133+ el = min(el, dt);
134+ tot = tot + dt;
135+ rep = rep + 1;
128136 end
129137 T(k, j) = el * 1000;
130138 E(k, j) = max(abs(v - ve));
src/chapters/ch11/pacman.mmodified+38−13View file
@@ -262,22 +262,47 @@ function v = bie_pac(d, g, px, py)
262262 [qx, qy, qnx, qny, qw, qkap] = pac_panels(d);
263263 n = numel(qx);
264264
265-DX = qx * ones(1, n) - ones(n, 1) * qx';
266-DY = qy * ones(1, n) - ones(n, 1) * qy';
267-R2 = DX.^2 + DY.^2;
268-R2(1:n+1:end) = 1; % keep 0/0 off the diagonal
269-K = ((ones(n, 1) * qnx') .* DX + (ones(n, 1) * qny') .* DY) ./ (2*pi*R2);
270-K(1:n+1:end) = -qkap / (4*pi); % the smooth-limit diagonal
271-A = -0.5 * eye(n) + K .* (ones(n, 1) * qw');
265+% The matrix is built by a double loop rather than by broadcasting the
266+% coordinate differences into n-by-n arrays. The loop is the faster shape in
267+% every environment this benchmark runs in: about twice as fast as implicit
268+% expansion under numbl and in MATLAB, and several times faster than forming
269+% the differences with outer products against ones(1,n), which allocates four
270+% n-by-n temporaries and multiplies two of them into being. It is also the
271+% only shape that writes each entry exactly once.
272+A = zeros(n, n);
273+c = 1 / (2*pi);
274+for j = 1:n
275+ xj = qx(j);
276+ yj = qy(j);
277+ nxj = qnx(j);
278+ nyj = qny(j);
279+ wj = qw(j);
280+ for i = 1:n
281+ dx = qx(i) - xj;
282+ dy = qy(i) - yj;
283+ A(i, j) = c * (nxj*dx + nyj*dy) / (dx*dx + dy*dy) * wj;
284+ end
285+ % the diagonal is the smooth limit of the kernel, -kappa/(4 pi), not a
286+ % singularity, and the jump relation contributes the -1/2
287+ A(j, j) = -0.5 - qkap(j) * wj / (4*pi);
288+end
272289
273290 sigma = A \ g(qx, qy);
274291
275-np = numel(px);
276-EX = px(:) * ones(1, n) - ones(np, 1) * qx';
277-EY = py(:) * ones(1, n) - ones(np, 1) * qy';
278-KE = ((ones(np, 1) * qnx') .* EX + (ones(np, 1) * qny') .* EY) ...
279- ./ (2*pi*(EX.^2 + EY.^2));
280-v = (KE .* (ones(np, 1) * qw')) * sigma;
292+% the layer potential at the target points, by the same loop
293+nt = numel(px);
294+v = zeros(nt, 1);
295+for k = 1:nt
296+ xk = px(k);
297+ yk = py(k);
298+ s = 0;
299+ for j = 1:n
300+ dx = xk - qx(j);
301+ dy = yk - qy(j);
302+ s = s + (qnx(j)*dx + qny(j)*dy) / (dx*dx + dy*dy) * qw(j) * sigma(j);
303+ end
304+ v(k) = c * s;
305+end
281306 end
282307
283308 function [qx, qy, qnx, qny, qw, qkap] = pac_panels(d)
src/chapters/ch11/text.mdmodified+16−2View file
@@ -39,6 +39,20 @@ julia bench/julia/bench_pacman.jl
3939 matlab -batch "run('bench/matlab/bench_pacman.m')"
4040 ```
4141
42-Each writes a small JSON file of resolutions, errors and times; the "load offline results" button under the work-precision plot draws them into it as dashed curves. The MATLAB harness does not port anything: it runs this chapter's driver and script, unchanged, in real MATLAB, which is the most direct numbl comparison available. The Python and Julia scripts are line-by-line ports, and the check that they are honest ports is in the plot itself: their error curves must land on top of the browser's, since error depends only on the mathematics, while their time curves may not.
42+Each writes a small JSON file of resolutions, errors and times; the "load offline results" button under the work-precision plot draws them into it as dashed curves. The MATLAB harness does not port anything: it runs this chapter's driver and script, unchanged, in real MATLAB, which is the most direct numbl comparison available. The Python and Julia scripts are ports, and the check that they are honest ports is in the plot itself: their error curves must land on the browser's, since error depends only on the mathematics, while their time curves may not.
4343
44-The reason for the upload, rather than a table of times measured on the authors' machine, is comparability. A timing from someone else's hardware shares no clock with the numbers your browser just produced; a timing from the scripts above, run on the machine you are reading this on, does. Loading a file recorded elsewhere is still allowed, and the legend will say what it was, but the one comparison the plot can defend is browser and native on the same machine. In general, on this problem, expect the native runs to shift the dashed curves left by a factor that is large for the dense-solve-heavy double layer and modest for the small least-squares fits of the harmonic basis; how large is precisely what the plot is for.
44+The reason for the upload, rather than a table of times measured on the authors' machine, is comparability. A timing from someone else's hardware shares no clock with the numbers your browser just produced; a timing from the scripts above, run on the machine you are reading this on, does. Loading a file recorded elsewhere is still allowed, and the legend will say what it was, but the one comparison the plot can defend is browser and native on the same machine.
45+
46+### What it takes for such a comparison to mean anything
47+
48+A cross-language benchmark measures the author's fluency in each language unless care is taken, and the first version of these scripts measured exactly that. Three disciplines are what make the curves comparable, and they are worth stating because they are easy to get wrong.
49+
50+**Each port is written in the fastest shape its own language offers, and those shapes are not the same shape.** The Nyström matrix can be built by a scalar double loop over its entries, or by broadcasting the coordinate differences into $n \times n$ arrays. The loop is the faster shape in three of the four environments here, by about a factor of two against implicit expansion in both MATLAB and numbl and by more than ten against it in Julia, and it is what the chapter's script and both of those ports use. In NumPy the same loop would be hopeless, so the Python port assembles with whole-array expressions, and does so in place: written as the obvious sequence of broadcast statements it allocates five $n \times n$ temporaries and spends longer moving memory than the dense solve spends on arithmetic. Holding the *code* fixed across languages would have been the unfair choice, not the fair one.
51+
52+**The same linear algebra must be asked for.** NumPy's `lstsq` is an SVD, while MATLAB's and Julia's backslash are both pivoted QR on a rectangular system. Leaving that alone had Python solving a harder problem than the others under the same label; the port now calls a QR driver. The errors then agree exactly on the double layer, and to within a factor of two at the $10^{-15}$ roundoff floor of the basis method, where the pivoting order is all that separates them.
53+
54+**Timing needs a protocol, applied everywhere.** Every environment here, including the page, keeps the fastest of several runs, with at least three runs whenever a solve is cheap enough to allow it. A single measurement of a fifty-millisecond solve varies by a factor of two with whatever else the machine is doing, which is enough to reorder neighbouring points and to make a language look slow that is not.
55+
56+With those three in place, on the author's machine (a 10-core i7-1355U), the four environments land within a factor of two of each other, and the remaining spread has causes one can name. MATLAB and Julia are within about ten per cent up to a $1300 \times 1300$ system; beyond that MATLAB pulls ahead, and its BLAS is using ten threads where Julia's default is six, which is why the legend records the thread count. Python trails by 1.3 to 1.7 times, essentially all of it in the assembly, which is the one stage where NumPy cannot express the fast shape. numbl runs about twice the time of native MATLAB on the same script, and much more than that on the smallest problems, where what is being measured is the fixed cost of a call rather than the arithmetic in it.
57+
58+Before the three disciplines were in place, the same scripts made Julia look about ten times faster than the other two. Nothing about that number was a property of Julia.