591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 1# Reducing spherical-harmonic transforms in $\Delta_\Gamma$
3**Summary.** Algorithm 4 costs 12 transforms per matvec (8 syntheses, 4 analyses). A flux-form
4reformulation, with weights chosen so that every analyzed field is smooth on $S^2$, evaluates the
5same operator in **6 transforms** (4 syntheses, 2 analyses). Notation follows `algos.pdf`.
7---
9## 1. Where the current cost comes from
11| Algorithm 4 line | Work | Transforms |
12|---|---|---|
13| 1 | $\partial_\theta u,\ \partial_\varphi u$ | 2 $\mathcal{S}$ |
14| 5 | analysis of 3 Cartesian components of $\nabla_\Gamma u$ | 3 $\mathcal{A}$ |
15| 6 | $\partial_\theta$ and $\partial_\varphi$ of each of those 3 components | 6 $\mathcal{S}$ |
16| 8 | final analysis | 1 $\mathcal{A}$ |
17| | | **12** |
19Two sources of waste:
211. The gradient is carried in **ambient $\mathbb{R}^3$ components** — 3 fields for an intrinsically
22 2-dimensional object.
232. Line 6 takes **both** derivatives of **each** component, where the divergence needs only one
24 derivative of each of two fluxes.
26There is also a possible free win independent of everything below: Algorithm 1 as specified returns
27all five derivatives. The Laplacian path needs only $\partial_\theta u$ and $\partial_\varphi u$;
28the second-derivative and mixed-derivative machinery is used exclusively by Algorithm 3 (curvature).
29If `surface_screened_laplacian` calls Algorithm 1 wholesale, it is doing 5 syntheses where 2 suffice
30at lines 1 and 6.
32---
34## 2. The constraint that shapes the solution
36The Cartesian design in Algorithm 4 exists to avoid pole singularities, and it is correct to do so.
37The relevant property is smoothness **as a scalar function on $S^2$**, since that is what controls
38SH coefficient decay and hence whether $\mathcal{A}$ is meaningful.
40| Quantity | Smooth on $S^2$? |
41|---|---|
42| $\partial_\varphi u$ | yes — exactly band-limited, eq. (2.2) |
43| $\sin\theta\,\partial_\theta u$ | yes — exactly band-limited, eq. (2.4) |
44| $\partial_\theta u$ | **no** — bounded, but $\varphi$-dependent limit at the poles |
45| $g_{\theta\theta},\ g^{\theta\theta},\ V_\theta,\ V_\varphi$ | **no** in general |
46| $(\nabla_\Gamma u)_x,\ (\nabla_\Gamma u)_y,\ (\nabla_\Gamma u)_z$ | yes |
48Concretely, for the ellipsoid $X = (a\sin\theta\cos\varphi,\ b\sin\theta\sin\varphi,\ c\cos\theta)$,
49$|X_\theta|^2 \to a^2\cos^2\varphi + b^2\sin^2\varphi$ as $\theta\to 0$: no limit exists.
51Algorithm 4 never analyzes anything in the "no" rows — the non-smooth quantities appear only as
52pointwise grid factors. **Any replacement must preserve this property.** The naive flux form
53$P = \sqrt{g}\,(g^{\theta\theta}u_\theta + g^{\theta\varphi}u_\varphi)$,
54$Q = \sqrt{g}\,(g^{\varphi\theta}u_\theta + g^{\varphi\varphi}u_\varphi)$ does not:
55on the round sphere with $u = x$, $Q = -\sin\varphi$, which is not a function on $S^2$.
57### The correct weighting
59$P$ and $Q$ are $\sqrt{g}$ times the contravariant components of $G := \nabla_\Gamma u$. Using
60$\det[X_\theta, X_\varphi, n] = -\sqrt{g}$:
62$$P = -\,G\cdot(X_\varphi \times n), \qquad \sin\theta\,Q = -\,G\cdot\big(n \times \sin\theta\,X_\theta\big).$$
64Every factor on the right is smooth on $S^2$: $G$ smooth, $n$ smooth, $X_\varphi$ smooth, and
65$\sin\theta\,X_\theta$ smooth because it is exactly band-limited by the recurrence already
66implemented. Hence $P$ and $\tilde{Q} := \sin\theta\,Q$ are analyzable, on the same footing and for
67the same structural reason as the Cartesian gradient components.
69The cross products are the smoothness certificate only — they are not needed in the code.
71---
73## 3. Precompute (once per surface, grid space)
75Replaces `_precompute_metric_quantities()`. From the embedding coefficients $\hat{X}^m_\ell$, obtain
76$X_\varphi$ and $\sin\theta\,X_\theta$ componentwise via Algorithm 1, where the latter is the
77**undivided** output of Algorithm 1 line 4, i.e. $\mathcal{S}(v^m_\ell)$ with no $/\sin\theta$.
78Then, pointwise:
80$$\tilde{g}_{\theta\theta} := |\sin\theta\,X_\theta|^2, \qquad
81 \tilde{g}_{\theta\varphi} := (\sin\theta\,X_\theta)\cdot X_\varphi, \qquad
82 g_{\varphi\varphi} := |X_\varphi|^2$$
84$$J := \frac{\sqrt{\tilde{g}_{\theta\theta}\,g_{\varphi\varphi} - \tilde{g}_{\theta\varphi}^{\,2}}}{\sin^2\theta}
85 \qquad\text{so that } \sqrt{\det g} = J\sin\theta$$
87(The radicand is $\sin^2\theta\det g = J^2\sin^4\theta$, so the square root is $J\sin^2\theta$ — hence
88$\sin^2\theta$, not $\sin\theta$, in the denominator.)
90Store four scalar grid arrays:
92$$p_1 = \frac{g_{\varphi\varphi}}{J\sin^2\theta}, \qquad
93 p_2 = -\frac{\tilde{g}_{\theta\varphi}}{J\sin^2\theta}, \qquad
94 q_2 = \frac{\tilde{g}_{\theta\theta}}{J\sin^2\theta}, \qquad
95 r = \frac{1}{J\sin^2\theta}$$
97All four are bounded: the $\sin^2\theta$ denominators cancel against vanishing numerators
98($\tilde{g}_{\theta\varphi} = O(\sin^2\theta)$, $g_{\varphi\varphi} = O(\sin^2\theta)$), the same
99finite limits the current $V_\theta, V_\varphi$ have. Note $p_1, p_2, q_2$ are bounded where
100$V_\varphi = O(1/\sin\theta)$ is not.
102**Three scalar arrays replace the six components of $V_\theta, V_\varphi$.** The surface
103representation is unchanged: $X_\theta, X_\varphi$ still come componentwise from $\hat{X}^m_\ell$
104via Algorithm 1.
106---
108## 4. The per-matvec algorithm
110Input $\{u^m_\ell\}$; output $\{(\Delta_\Gamma u)^m_\ell\}$.
112| # | Step | Transforms |
113|---|---|---|
114| 1 | $v^m_\ell \leftarrow \alpha^+(\ell-1,m)u^m_{\ell-1} + \alpha^-(\ell+1,m)u^m_{\ell+1}$; $A \leftarrow \mathcal{S}(v^m_\ell)$ | $\mathcal{S}$ |
115| 2 | $B \leftarrow \mathcal{S}(im\,u^m_\ell)$ | $\mathcal{S}$ |
116| 3 | $P \leftarrow p_1 A + p_2 B$, $\tilde{Q} \leftarrow p_2 A + q_2 B$ — pointwise | — |
117| 4 | $\hat{P} \leftarrow \mathcal{A}(P)$, $\hat{\tilde{Q}} \leftarrow \mathcal{A}(\tilde{Q})$ | 2 $\mathcal{A}$ |
118| 5 | $s^m_\ell \leftarrow \alpha^+(\ell-1,m)\hat{P}^m_{\ell-1} + \alpha^-(\ell+1,m)\hat{P}^m_{\ell+1} + im\,\hat{\tilde{Q}}^m_\ell$ | — |
119| 6 | $\Delta_\Gamma u \leftarrow r \cdot \mathcal{S}(s^m_\ell)$ | $\mathcal{S}$ |
120| 7 | $\{(\Delta_\Gamma u)^m_\ell\} \leftarrow \mathcal{A}(\Delta_\Gamma u)$; zero $\ell \ge L-2$ | $\mathcal{A}$ |
122**4 syntheses + 2 analyses = 6**, versus 12.
124- $A$ and $B$ are exactly $\sin\theta\,\partial_\theta u$ and $\partial_\varphi u$.
125- Steps 1 and 5 use the **same** precomputed $\alpha^\pm$ table; step 5 is the adjoint-style reuse
126 of the shift already implemented for step 1. Adding the two flux contributions in coefficient
127 space before synthesizing is what saves the final pair of transforms.
128- The **only** division by $\sin\theta$ anywhere is folded into $p_1, p_2, q_2, r$ at precompute
129 time. The per-matvec path contains none.
131---
133## 5. Numerical trade-off
135Both schemes contain two powers of $\sin\theta$ division in total. What differs is **placement**.
137- **Algorithm 4** spends them in separate stages, one before line 5 and one after. The intervening
138 analysis suppresses the polar spike: Gauss–Legendre weights give $w_1 = O(L^{-2})$ at the polar
139 ring, so a grid error of $\varepsilon L$ there contributes
140 $\sim L^{-2}\cdot L^{1/2}\cdot \varepsilon L = \varepsilon L^{-1/2}$ to any coefficient. Stage 2
141 then starts from clean coefficients and incurs a *fresh* $\varepsilon L$. The two amplifications
142 never multiply. Net grid-space relative error: $\varepsilon L$.
143- **The new scheme** has no division at all through step 5, then pays for both powers at once in
144 $r = O(L^2)$ at step 6 — one event, with no intervening analysis to break it in half. Net
145 grid-space relative error: $\varepsilon L^2$.
147The mechanism: with $N \approx L+1$ Gauss–Legendre nodes, $1 - x_1 = O(N^{-2})$ so
148$\sin\theta_1 = O(N^{-1})$. Since $s = J\sin^2\theta\,\Delta_\Gamma u$ is $O(L^{-2})$ at the polar
149ring but $O(1)$ over the bulk, and synthesis commits roundoff scaled by the field's *global* size at
150every node alike, multiplying by $r \sim L^2$ recovers the signal and inflates the noise.
152| | divisions | placement | grid-space relative error |
153|---|---|---|---|
154| Algorithm 4 | $\sin\theta$, $\sin\theta$ | separated by $\mathcal{A}$ | $\varepsilon L$ |
155| Six-transform | $\sin^2\theta$ | all at the end | $\varepsilon L^2$ |
157**This likely does not reach the returned coefficients.** Step 7's analysis suppresses the spike
158exactly as line 5 does today: $L^{-2}\cdot L^{1/2}\cdot\varepsilon L^2 = \varepsilon L^{1/2}$,
159comparable to the ordinary $\varepsilon\sqrt{L}$ accumulation of a transform pair — and the new
160scheme runs half as many transforms, lowering that baseline. Inside the implicit solve, GMRES sees
161only coefficients, so the extra power should be invisible.
163It matters only if grid values of $\Delta_\Gamma u$ are consumed directly: a nonlinear reaction
164term, max-norm diagnostics, or an adaptive error estimator.
166Algorithm 1 line 7 already divides by $\sin^2\theta$, so the code is exposed to $\varepsilon L^2$
167today — just on the second-derivative path, which the Laplacian never touches.
169---
171## 5a. float32 / WebGPU
173Target is WebGPU, which is float32-only: $\varepsilon = 2^{-24} \approx 6\times10^{-8}$ (spacing
174$2^{-23} \approx 1.2\times10^{-7}$). No float64 fallback exists on device. All estimates in §5 are
175linear in $\varepsilon$, so they scale directly:
177| | $\varepsilon\sqrt{L}$ (coeffs) | $\varepsilon L$ (Alg. 4 grid) | $\varepsilon L^2$ (new, grid) |
178|---|---|---|---|
179| $L=64$ | $5\times10^{-7}$ | $4\times10^{-6}$ | $2\times10^{-4}$ |
180| $L=128$ | $7\times10^{-7}$ | $8\times10^{-6}$ | $1\times10^{-3}$ |
181| $L=256$ | $1\times10^{-6}$ | $1.5\times10^{-5}$ | $4\times10^{-3}$ |
183**Coefficient space is fine** (~$10^{-6}$), which is the floor a float32 iterative solve sits at
184anyway. **Grid space is not**: 0.1–0.4% relative on the polar rings at $L\ge128$. For a
185reaction–diffusion solver this matters only if grid-space $\Delta_\Gamma u$ is consumed outside the
186matvec. If the IMEX splitting evaluates $f(u)$ from $u$ on the grid (typical), it never is.
188**Two float32-specific arguments in favour of the new scheme:**
190- Baseline SHT roundoff accumulates per transform ($\sim\varepsilon\sqrt{L}$ to $\varepsilon L$
191 each). Running 6 transforms instead of 12 halves that accumulation. On the coefficient-space error
192 GMRES actually sees, this plausibly outweighs the polar term — the new scheme may be *more*
193 accurate end-to-end in float32. Not asserted without measurement.
194- 3 weight arrays instead of 6 halves per-matvec texture/buffer traffic. On GPU that is often the
195 real bottleneck, independent of arithmetic.
197**Mitigations available without float64:**
199- **Pairwise or blocked summation in the Legendre sum over $\ell$.** The single highest-value
200 float32 change, and it benefits the existing code too. See §5b.
201- **Double-float (`f32x2`) arithmetic** for the pointwise steps 3 and 6 if needed — cheap, no
202 transforms involved. Does not help with transform roundoff, which is the dominant term, so try
203 summation order first.
204- **CPU precompute in float64.** JS `Number` is float64, so §3 can run on the CPU regardless of
205 WebGPU's limits, with float32 weights uploaded. Cost is CPU-side SHTs plus upload, paid once per
206 surface update; viable if the surface evolves slowly or is prescribed analytically, likely too
207 slow if the metric is rebuilt every timestep. Per the correction below, this is probably
208 unnecessary.
209- **Cap $L$.** All the error terms grow with $L$; float32 sets a practical ceiling that float64
210 would not.
212---
214## 5b. Summation order in the Legendre transform
216This is orthogonal to the 12→6 change, applies equally to the current code, and in float32 is
217probably worth more than the transform-count reduction. Do it first and independently, so its effect
218can be measured on its own.
220**Why.** Every $\varepsilon L$ and $\varepsilon L^2$ in §5 rides on the per-transform roundoff
221floor, and in float32 that floor is set by *how the sums are accumulated*, not by the mathematics.
222For each $(m, \theta_i)$ the synthesis evaluates
224$$u^m(\theta_i) = \sum_{\ell=|m|}^{L} u^m_\ell\,\bar P^m_\ell(\cos\theta_i),$$
226an $O(L)$-term sum. Error growth by accumulation strategy, for an $N$-term sum:
228| Strategy | Worst case | Typical (random signs) |
229|---|---|---|
230| Sequential | $\varepsilon N$ | $\varepsilon\sqrt{N}$ |
231| Pairwise / tree | $\varepsilon\log_2 N$ | $\varepsilon\sqrt{\log_2 N}$ |
232| Kahan compensated | $\varepsilon$ (+ $O(\varepsilon^2 N)$) | $\varepsilon$ |
234At $L=256$ in float32 that is the difference between $\sim1.5\times10^{-5}$ and $\sim5\times10^{-7}$
235per transform — more than an order of magnitude, for no change in operation count.
237**On GPU this may already be partly free.** A workgroup tree reduction over $\ell$ *is* pairwise
238summation. The failure mode is a serial `for` loop over $\ell$ inside a single thread, which is the
239natural way to write the shader if each thread owns one $(m,\theta_i)$ pair and is exactly the
240$\varepsilon N$ row above. Check which shape the kernel has before assuming anything.
242**Where it applies.**
244- Synthesis $\mathcal{S}$: the sum over $\ell$, as above. The $\varphi$-direction FFT is already
245 tree-structured and needs no attention.
246- Analysis $\mathcal{A}$: the quadrature sum over latitude nodes $\theta_i$ carries the identical
247 problem and the identical fix. It also matters more here, because this is the step relied on in
248 §5 to suppress the polar spike — a noisy quadrature sum weakens exactly the mechanism the
249 six-transform scheme depends on.
251**Practical notes.**
253- Blocked summation (accumulate in blocks of 8–32, then combine) captures most of the pairwise
254 benefit with a simpler kernel and better register behaviour than a full tree.
255- Kahan costs ~4 flops per term and is usually bandwidth-hidden on GPU; worth benchmarking rather
256 than assuming it is too expensive.
257- For $m>0$ near the poles, $\bar P^m_\ell(\cos\theta)$ spans many orders of magnitude across $\ell$.
258 Summing smallest-magnitude-first helps, and is nearly free here because the terms are already
259 roughly ordered by $\ell$.
260- Standard stable recurrences for $\bar P^m_\ell$ (and guarding their under/overflow in float32's
261 narrower exponent range) are a separate prerequisite — no summation strategy rescues inaccurate
262 Legendre values.
264**Measurement.** Transform a band-limited field forward then back and compare to the input, in
265float32, sweeping $L\in\{64,128,256\}$. Sequential accumulation shows error growing roughly linearly
266in $L$; pairwise shows near-flat growth. This isolates the transform floor from everything else in
267§7 and should be run before the validation gate there, since it sets the baseline that gate is
268measured against.
270### Measured (2026-08-04, Dawn/Metal, `scripts/sht-accuracy.ts`)
272The sweep was run and the summation-order changes tried. Outcome: **withdrawn — the floor here is
273not summation-limited.**
275| $L$ | grid | rel-$L_2$ roundtrip | worst degree |
276|---|---|---|---|
277| 63 | 64×128 | $3.4\times10^{-6}$ | $\ell=62$: $4.5\times10^{-6}$ |
278| 127 | 128×256 | $4.7\times10^{-6}$ | $\ell=110$: $6.4\times10^{-6}$ |
279| 255 | 256×512 | $1.1\times10^{-5}$ | $\ell=246$: $1.4\times10^{-5}$ |
281- **The analysis side already sums pairwise.** The quadrature over latitudes is a workgroup
282 tree/subgroup reduction (`leg_analys`); only the synthesis has the serial per-thread $\ell$-loop.
283- **Kahan is unavailable on WebGPU in practice.** Dawn/Metal compiles WGSL with fast-math: a probe
284 kernel evaluates $((10^8 + 1) - 10^8) - 1$ to $0$, so the compensation folds away and Kahan
285 compiles to plain summation (bit-identical results, verified).
286- **Blocked summation (B=16) in the synthesis $\ell$-loop moved nothing**: $1.128\times10^{-5}
287 \to 1.128\times10^{-5}$ at $L=255$ (low digits shift, confirming the reordering was live), while
288 costing ~5% per round trip at $L=255$. Reverted.
289- **Diagnosis:** the worst error concentrates at the top degrees — the signature of the Legendre
290 *recurrence* error (chains of length $\sim\ell$), not of $\ell$-uniform accumulation noise. This
291 is the "standard stable recurrences are a separate prerequisite" caveat above: the floor is set
292 by the accuracy of the $\bar P^m_\ell$ values themselves, and no summation strategy touches it.
293- The measured floor ($\sim\varepsilon L^{0.85}$, $1.1\times10^{-5}$ at $L=255$) is what the §7
294 validation gate should be read against.
296---
298## 6. Code changes
300| Location | Change |
301|---|---|
302| `src/surface_gradient/partial_derivatives` | Expose $\mathcal{S}(v^m_\ell)$ **pre-division** (flag or separate entry point). Needed by both the precompute and step 1. |
303| `SurfaceDiffOperator._precompute_metric_quantities()` | Return `p1, p2, q2, r` instead of `V_theta, V_phi`, per §3. |
304| `src/surface_screened_laplacian::surface_screened_laplacian()` | Replace body with §4. Both `for i in {x,y,z}` loops disappear. |
305| `SurfaceDiffOperator._precompute_curvature()`, Algorithm 3 | **Unchanged.** Still needs $X_{\theta\theta}, X_{\theta\varphi}, X_{\varphi\varphi}$ and the full Algorithm 1. |
306| `src/timestepping::make_implicit_op()`, Algorithm 5 | **Unchanged.** Only what line 8 calls changes. |
307| `src/real_embedding.py` | **Unchanged.** |
309The deprecated `SurfaceDiffOperator` methods for $\Delta_\Gamma$ and $(I + c\Delta_\Gamma)$ are the
310natural place to keep the old path as a reference implementation for the validation below.
312---
314### Correction: precompute conditioning
316An earlier draft claimed the polar relative error in $\tilde g_{\theta\theta}$ is $\varepsilon L^2$,
317making float64 precompute essential. That was wrong by a factor of $L$, in the safe direction.
318$\tilde g_{\theta\theta}$ is not synthesized directly; it is the square of $\sin\theta\,X_\theta$,
319which *is* synthesized, is $O(\sin\theta)$ at the poles, and carries absolute error $\varepsilon$ —
320so relative error $\varepsilon L$, preserved (up to a factor 2) by squaring. Same for
321$g_{\varphi\varphi}$ and $\tilde g_{\theta\varphi}$. The determinant combination is $O(\sin^4\theta)$
322and so are both of its terms, so there is no extra cancellation generically; $J$ inherits
323$\sim\varepsilon L$, i.e. $\sim10^{-5}$ at $L=128$ in float32. Acceptable.
325Caveat: this assumes the difference is not small compared to its terms, which fails if $X_\theta$
326becomes nearly parallel to $X_\varphi$ (near-degenerate parametrization). Worth a runtime check on
327$\det g$ if the surface can deform that far.
329Precompute error is also a *fixed* perturbation, identical every matvec, so it perturbs which
330operator is being solved but injects no noise into the Krylov space — GMRES converges normally.
332---
334## 7. Validation, in order
336At $\varepsilon=6\times10^{-8}$ there is no margin for the $\varepsilon\sqrt{L}$ suppression estimate
337in §5 to be off by an order of magnitude. Step 2 is a **gate**, not a confirmation.
3391. **Smoothness check (do this first).** On a deformed, non-axisymmetric surface, form $P$ and
340 $\tilde{Q}$ on the grid and compare their SH coefficient decay against $(\nabla_\Gamma u)_x$ from
341 the current code. Matching tails confirm both are genuinely smooth on $S^2$. If this fails,
342 nothing else is worth doing. Run this in float64 on CPU — it is a mathematical check, not a
343 precision one.
3442. **Coefficient-space diff in float32 at production $L$**, against a float64 CPU reference
345 implementation of Algorithm 4. Landing near $10^{-6}$ means the suppression argument holds.
346 Landing near $10^{-4}$ means the polar spike is surviving the analysis and the polar rings need
347 separate handling.
3483. **Sweep $L \in \{64,128,256\}$** and fit the growth exponent of (2). Flat-ish confirms
349 suppression; growth like $L^2$ means it is not working.
3504. **Grid-space max-norm diff near the poles.** Expect the extra power of $L$ here. If only this
351 grows and (2) stays flat, the scheme is fine for use inside the implicit solve.
3525. **GMRES iteration count and final achieved residual**, float32, versus the current code. The
353 operator is the same, so iterations should be unchanged; a stall above tolerance that does not
354 occur in float64 indicates the matvec noise floor is binding.
356---
358## 8. Suggestions considered and withdrawn
360- **Splitting $r$ across steps 3 and 6** to keep $\sin^1$ scaling. Does not work: $P/\sin\theta$ is
361 not smooth on $S^2$ (round sphere, $u = x$: $P = \sin\theta\cos\theta\cos\varphi$, so
362 $P/\sin\theta \to \cos\varphi$ at the pole). That $\sin\theta$ must stay in $r$. Algorithm 4 *can*
363 split because its intermediate — the Cartesian gradient — is smooth; that smoothness is precisely
364 what the six extra transforms buy.
365- **Weighting by $J$ to get an SPD operator and use PCG.** Avoiding the $1/\sqrt{g}$ makes
366 $\mathcal{L}$ self-adjoint, but the mass term becomes multiplication by $J$, costing its own
367 synthesis/analysis pair. A wash on transform count; worth it only if the CG properties themselves
368 are wanted. Discrete symmetry would also hold only to quadrature accuracy unless products are
369 dealiased (3/2 rule).
370- **Mixed precision (float64 for step 6's synthesis and the $r$ multiply only).** Unavailable:
371 WebGPU is float32-only. Superseded by the mitigations in §5a.
372- If the $\varepsilon L^2$ ever does bind, the remaining remedies are a shifted or uniform-in-$\theta$
373 latitude grid (removing the $O(L^{-2})$ node clustering) or a separate local formula for the polar
374 rings — both more work than the transform savings justify without a specific reason. In float32
375 the ceiling on $L$ may bind first and be the cheaper accommodation.
377## 9. Related: external vector transforms
379Steps 1–2 and 4–5 together are a vector/spin-weighted spherical harmonic transform. If SHTns
380(`spat_to_SHsphtor`, `SHsphtor_to_spat`) or SPHEREPACK (`gradgs`, `divgs`) can be linked, the
381gradient and divergence each become a single library call, the pole divisions are handled internally,
382and the hand-rolled $\alpha^\pm$ recurrences on this path are no longer needed.
384## 10. Beyond transform count
386The other lever is iteration count rather than cost per iteration. Since $M^{-1}A$ approaches
387multiplication by $1/J$ at high $\ell$, folding a mean or smoothed $J$ into the preconditioner could
388reduce GMRES iterations by more than any of the above reduces transforms.