1% Solve the implicit diffusion system
2%
3% (I - dtD*lap_g) X = B, i.e. A*X = B with A*x = M.*x - dtD*dlap(x)
4%
5% by right-preconditioned GMRES(niter) — one Arnoldi sweep of niter
6% iterations, no restart — with the round-sphere operator M = 1 + dtD*lam as
7% the preconditioner. Same operator and preconditioner as the other two
8% solvers; what GMRES adds over bicgstab is optimality (the residual is
9% minimized over the whole Krylov space, monotonically nonincreasing) at the
10% price of storing the basis and the O(niter^2) orthogonalization sweep.
11%
12% The bookkeeping no other solver needs — the basis, the Hessenberg matrix,
13% the Givens rotations — lives in banks and small matrices accessed through
14% getslab/setslab and getat/setat (src/mgpu/externals.ts): functional
15% updates the planner turns into static-offset buffer copies once the
16% unrolled loop's variable makes every index a literal. That is also why the
17% triangular loops below (`for i = 1:j`) compile: each unrolled j plans its
18% own inner trip count.
19%
20% Scalars are GPU-resident 1-element values throughout, as in bicgstab, and
21% every ratio is guarded the same way (a*b/(b*b + 1e-30), and 1/sqrt(x) as
22% 1/sqrt(x + 1e-30)), so a converged or broken-down iteration contributes
23% zero coefficients instead of NaNs. sqrt takes abs() of its argument
24% because the compiler cannot see that <w, w> is nonnegative.
25%
26% niter is fixed at compile time. X0 is the round-sphere answer, and the
27% correction added at the end is X = X0 + M \ (V*y) accumulated slab by
28% slab.
30function X = gmres(B, dtD, lam, filt, wlm, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, nlm, niter)
31 M = 1 + dtD * lam;
32 X = B ./ M;
33 dL0 = dlap(X, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, lam);
34 R = B - (M .* X - dtD * dL0);
36 % The Krylov basis (niter+1 spectral fields), the Hessenberg matrix, the
37 % rotated right-hand side, the Givens coefficients, and the solve's y.
38 VB = zeros(2, nlm * (niter + 1));
39 H = zeros(niter + 1, niter);
40 g = zeros(niter + 1, 1);
41 c = zeros(niter, 1);
42 s = zeros(niter, 1);
43 y = zeros(niter, 1);
45 Rw = R .* wlm;
46 r2 = dot(Rw, R);
47 nr = sqrt(abs(r2));
48 invnr = nr / (r2 + 1e-30);
49 V1 = R * invnr;
50 VB = setslab(VB, V1, 1);
51 g = setat(g, nr, 1);
53 for j = 1:niter
54 % w = A * M \ v_j
55 Vj = getslab(VB, j);
56 Zj = Vj ./ M;
57 dLz = dlap(Zj, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, lam);
58 W = M .* Zj - dtD * dLz;
60 % Modified Gram-Schmidt against every basis vector so far.
61 for i = 1:j
62 Vi = getslab(VB, i);
63 Viw = Vi .* wlm;
64 hij = dot(Viw, W);
65 H = setat(H, hij, i, j);
66 W = W - hij * Vi;
67 end
68 Ww = W .* wlm;
69 w2 = dot(Ww, W);
70 hn = sqrt(abs(w2));
71 H = setat(H, hn, j + 1, j);
72 invh = hn / (w2 + 1e-30);
73 Vn = W * invh;
74 VB = setslab(VB, Vn, j + 1);
76 % Apply the previous Givens rotations to column j, then form the new
77 % one that zeroes H(j+1, j), and rotate g with it.
78 for i = 1:j-1
79 a = getat(H, i, j);
80 b = getat(H, i + 1, j);
81 ci = getat(c, i);
82 si = getat(s, i);
83 t1 = ci * a + si * b;
84 t2 = ci * b - si * a;
85 H = setat(H, t1, i, j);
86 H = setat(H, t2, i + 1, j);
87 end
88 a = getat(H, j, j);
89 b = getat(H, j + 1, j);
90 rr = a * a + b * b;
91 invr = 1 / sqrt(abs(rr) + 1e-30);
92 cj = a * invr;
93 sj = b * invr;
94 c = setat(c, cj, j);
95 s = setat(s, sj, j);
96 t1 = cj * a + sj * b;
97 H = setat(H, t1, j, j);
98 gj = getat(g, j);
99 t3 = cj * gj;
100 t4 = -(sj * gj);
101 g = setat(g, t3, j);
102 g = setat(g, t4, j + 1);
103 end
105 % Back-substitute the rotated (upper triangular) system H*y = g.
106 for j = niter:-1:1
107 acc = getat(g, j);
108 for i = j+1:niter
109 rji = getat(H, j, i);
110 yi = getat(y, i);
111 acc = acc - rji * yi;
112 end
113 rjj = getat(H, j, j);
114 yj = (acc * rjj) / (rjj * rjj + 1e-30);
115 y = setat(y, yj, j);
116 end
118 % X = X0 + M \ (V * y), slab by slab.
119 for j = 1:niter
120 Vj = getslab(VB, j);
121 yj = getat(y, j);
122 X = X + yj * (Vj ./ M);
123 end
124end