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 preconditioned BiCGSTAB (van der Vorst), with the round-sphere operator
6% M = 1 + dtD*lam as the preconditioner — the same M solvers/richardson.m
7% inverts, applied here as M\v inside the Krylov recurrence. Same operator,
8% different solver: where Richardson converges only while the geometric
9% correction stays small against M, BiCGSTAB builds a Krylov space and
10% converges on configurations well outside that radius, at the price of two
11% dlap evaluations per iteration instead of one, plus three dot products.
12%
13% The scalars (rho, alpha, omega) are GPU-resident 1-element values: `dot` is
14% a reduction dispatch and the recurrences on its results compile to
15% 1-element kernels, so the whole solve is still one command stream with no
16% CPU in the loop. Inner products are taken in the half-spectrum weighting
17% wlm (m > 0 counts twice), so they are the real L2 inner products on the
18% sphere; the weight is folded into the fixed shadow residual once, and into
19% t per iteration.
20%
21% niter is fixed at compile time: no residual test, and so no early exit —
22% which for BiCGSTAB matters at the *converged* end, where the residual is
23% fp32 noise and every textbook ratio is 0/0. Each ratio a/b is therefore
24% written in the guarded form a*b/(b*b + 1e-30): identical to a/b (to a
25% couple of ulp) whenever b is meaningfully sized, and 0 when it is not, so
26% a converged or broken-down iteration degrades to a stationary one
27% (coefficients 0, X carried unchanged) instead of poisoning the state with
28% NaNs. X0 is the round-sphere answer, so niter = 0 computes richardson's
29% starting divide.
31function X = bicgstab(B, dtD, lam, filt, wlm, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, niter)
32 M = 1 + dtD * lam;
33 X = B ./ M;
34 dL0 = dlap(X, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, lam);
35 R = B - (M .* X - dtD * dL0);
36 Rh = R .* wlm;
37 P = 0 * R;
38 V = 0 * R;
39 rho = 1;
40 alpha = 1;
41 omega = 1;
42 for k = 1:niter
43 rho1 = dot(Rh, R);
44 beta = (rho1 * rho) / (rho * rho + 1e-30) * ((alpha * omega) / (omega * omega + 1e-30));
45 P = R + beta * (P - omega * V);
46 Ph = P ./ M;
47 dLp = dlap(Ph, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, lam);
48 V = M .* Ph - dtD * dLp;
49 rv = dot(Rh, V);
50 alpha = (rho1 * rv) / (rv * rv + 1e-30);
51 S = R - alpha * V;
52 Sh = S ./ M;
53 dLs = dlap(Sh, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, lam);
54 T = M .* Sh - dtD * dLs;
55 Tw = T .* wlm;
56 ts = dot(Tw, S);
57 tt = dot(Tw, T);
58 omega = (ts * tt) / (tt * tt + 1e-30);
59 X = X + alpha * Ph + omega * Sh;
60 R = S - omega * T;
61 rho = rho1;
62 end
63end