1% Five ways to solve the five-point system without factorising anything.
2%
3% Every method works only through the action of
4%
5% (A u)(i,j) = 4 u(i,j) - u(i-1,j) - u(i+1,j) - u(i,j-1) - u(i,j+1)
6%
7% on the interior points, with the boundary values held fixed in the array. No
8% matrix is formed, so the storage is the grid and nothing else. What separates
9% the methods is how fast the error goes away: Jacobi and Gauss-Seidel need
10% O(N^2) sweeps, optimal SOR and conjugate gradients O(N), and multigrid a
11% number of cycles that does not grow with N at all.
12%
13% UE is passed in only so that the true error can be recorded next to the
14% residual. No method uses it.
15%
16% A note on the shape of this file. Each method has its own loop rather than
17% one loop with a switch inside it, and multigrid keeps every level in one flat
18% vector rather than a cell array of grids. Both keep the *timed* loop
19% type-stable for numbl's compiler: an iterate assigned from several different
20% functions, or a body that indexes a cell, is not something it can specialise,
21% and here that loop is the measurement.
23function names = method_names()
24names = {'Jacobi', 'Gauss-Seidel', 'SOR', 'conjugate gradients', 'multigrid V-cycles'};
25end
27function H = run_iter(mid, N, U, UE, nit, tol)
28if mid == 1
29 H = loop_jacobi(N, U, UE, nit, tol);
30elseif mid == 2
31 H = loop_sor(N, U, UE, nit, tol, 1);
32elseif mid == 3
33 H = loop_sor(N, U, UE, nit, tol, 2 / (1 + sin(pi/N)));
34elseif mid == 4
35 H = loop_cg(N, U, UE, nit, tol);
36else
37 H = loop_mg(N, U, UE, nit, tol);
38end
39end
41% ------------------------------------------------------------ the methods ---
43function H = loop_jacobi(N, U, UE, nit, tol)
44% Every point replaced by the average of its four old neighbours. Error
45% components with wavelength comparable to the domain are damped by a factor
46% 1 - O(h^2) per sweep, so the sweep count grows like N^2.
47res = NaN(1, nit+1);
48err = NaN(1, nit+1);
49tms = NaN(1, nit+1);
50nf = 1;
51ss0 = res_sumsq(U, N);
52res(1) = ss0;
53err(1) = maxsq(U - UE);
54tms(1) = 0;
55stop = tol * tol * ss0;
56tacc = 0;
57for it = 1:nit
58 t1 = tic;
59 U = jacobi_sweep(U, N);
60 tacc = tacc + toc(t1);
61 res(it+1) = res_sumsq(U, N);
62 err(it+1) = maxsq(U - UE);
63 tms(it+1) = tacc * 1000;
64 nf = it + 1;
65 if res(it+1) < stop
66 break;
67 end
68end
69H = hist_pack(res, err, tms, ss0, nf);
70end
72function H = loop_sor(N, U, UE, nit, tol, omega)
73% omega = 1 is Gauss-Seidel: new values are used as soon as they are available,
74% which costs nothing and halves the sweep count. Over-relaxing with the
75% optimal omega = 2/(1 + sin(pi h)) changes the rate from 1 - O(h^2) to
76% 1 - O(h). That optimal value is known only because the eigenvalues of this
77% particular matrix are known.
78res = NaN(1, nit+1);
79err = NaN(1, nit+1);
80tms = NaN(1, nit+1);
81nf = 1;
82ss0 = res_sumsq(U, N);
83res(1) = ss0;
84err(1) = maxsq(U - UE);
85tms(1) = 0;
86stop = tol * tol * ss0;
87tacc = 0;
88for it = 1:nit
89 t1 = tic;
90 U = sor_sweep(U, N, omega);
91 tacc = tacc + toc(t1);
92 res(it+1) = res_sumsq(U, N);
93 err(it+1) = maxsq(U - UE);
94 tms(it+1) = tacc * 1000;
95 nf = it + 1;
96 if res(it+1) < stop
97 break;
98 end
99end
100H = hist_pack(res, err, tms, ss0, nf);
101end
103function H = loop_cg(N, U, UE, nit, tol)
104% Conjugate gradients. A is symmetric positive definite, so CG applies; it
105% needs one operator application per iteration and four vectors of storage.
106% Its convergence is governed by sqrt(kappa) rather than kappa, which for this
107% matrix means O(N) iterations, like optimal SOR, but without having to know
108% anything about the spectrum in advance.
109res = NaN(1, nit+1);
110err = NaN(1, nit+1);
111tms = NaN(1, nit+1);
112nf = 1;
113R = neg_residual(U, N); % the right-hand side is zero at interior points
114P = R;
115rs = dot2(R, R);
116ss0 = rs;
117res(1) = ss0;
118err(1) = maxsq(U - UE);
119tms(1) = 0;
120stop = tol * tol * ss0;
121tacc = 0;
122for it = 1:nit
123 t1 = tic;
124 AP = applyA(P, N);
125 al = rs / dot2(P, AP);
126 U = U + al * P;
127 R = R - al * AP;
128 rsn = dot2(R, R);
129 P = R + (rsn / rs) * P;
130 rs = rsn;
131 tacc = tacc + toc(t1);
132 res(it+1) = rsn;
133 err(it+1) = maxsq(U - UE);
134 tms(it+1) = tacc * 1000;
135 nf = it + 1;
136 if res(it+1) < stop
137 break;
138 end
139end
140H = hist_pack(res, err, tms, ss0, nf);
141end
143function H = loop_mg(N, U, UE, nit, tol)
144% Geometric multigrid. Gauss-Seidel is a poor solver but an excellent smoother:
145% two sweeps remove the oscillatory part of the error and leave something
146% smooth, which is well represented on a grid twice as coarse and much cheaper
147% to correct there. Recursing on that idea gives a V-cycle whose cost is
148% proportional to the number of unknowns and whose error reduction per cycle is
149% a constant independent of N.
150nlev = 1;
151while mod(N, 2^nlev) == 0 && N / 2^nlev >= 2
152 nlev = nlev + 1;
153end
154ns = zeros(nlev, 1);
155offs = zeros(nlev, 1);
156tot = 0;
157nn = N;
158for l = 1:nlev
159 ns(l) = nn;
160 offs(l) = tot;
161 tot = tot + (nn+1)^2;
162 nn = nn / 2;
163end
165% one flat vector holds every level: the iterate in the first half, the
166% right-hand side in the second
167W = zeros(2*tot, 1);
168W = flatten_into(W, U, 0, N);
169WE = flatten_into(zeros((N+1)^2, 1), UE, 0, N);
171res = NaN(1, nit+1);
172err = NaN(1, nit+1);
173tms = NaN(1, nit+1);
174nf = 1;
175ss0 = res_sumsq_flat(W, 0, tot, N);
176res(1) = ss0;
177err(1) = maxsq_flat(W, WE, N);
178tms(1) = 0;
179stop = tol * tol * ss0;
180tacc = 0;
181for it = 1:nit
182 t1 = tic;
183 W = vcycle(W, offs, ns, nlev, tot, 2, 2);
184 tacc = tacc + toc(t1);
185 res(it+1) = res_sumsq_flat(W, 0, tot, N);
186 err(it+1) = maxsq_flat(W, WE, N);
187 tms(it+1) = tacc * 1000;
188 nf = it + 1;
189 if res(it+1) < stop
190 break;
191 end
192end
193H = hist_pack(res, err, tms, ss0, nf);
194end
196% ---------------------------------------------------------------- kernels ---
198function U = jacobi_sweep(U, N)
199V = U;
200for j = 2:N
201 for i = 2:N
202 U(i,j) = 0.25 * (V(i-1,j) + V(i+1,j) + V(i,j-1) + V(i,j+1));
203 end
204end
205end
207function U = sor_sweep(U, N, omega)
208for j = 2:N
209 for i = 2:N
210 U(i,j) = (1-omega)*U(i,j) + 0.25*omega*(U(i-1,j) + U(i+1,j) + U(i,j-1) + U(i,j+1));
211 end
212end
213end
215function Y = applyA(U, N)
216Y = zeros(N+1, N+1);
217for j = 2:N
218 for i = 2:N
219 Y(i,j) = 4*U(i,j) - U(i-1,j) - U(i+1,j) - U(i,j-1) - U(i,j+1);
220 end
221end
222end
224function Y = neg_residual(U, N)
225Y = zeros(N+1, N+1);
226for j = 2:N
227 for i = 2:N
228 Y(i,j) = -(4*U(i,j) - U(i-1,j) - U(i+1,j) - U(i,j-1) - U(i,j+1));
229 end
230end
231end
233function s = dot2(A, B)
234s = 0;
235n = numel(A);
236for k = 1:n
237 s = s + A(k) * B(k);
238end
239end
241function ss = res_sumsq(U, N)
242% The sum of squares, not the norm: a square root taken here would widen the
243% inferred type of the accumulator and slow this loop by a factor of 20.
244ss = 0;
245for j = 2:N
246 for i = 2:N
247 r = 4*U(i,j) - U(i-1,j) - U(i+1,j) - U(i,j-1) - U(i,j+1);
248 ss = ss + r*r;
249 end
250end
251end
253function H = hist_pack(res, err, tms, ss0, nf)
254% Both histories are accumulated as squares and rooted here, over the first nf
255% columns — the part the iteration actually filled; the rest stays NaN, which
256% is what the page draws as a gap. Rooting inside the loop instead costs about
257% 60 per cent on the sweep itself, because sqrt may return a complex value and
258% the history array is then no longer a plain real array to numbl.
259r = res;
260e = err;
261r(1:nf) = sqrt(res(1:nf) / ss0);
262e(1:nf) = sqrt(err(1:nf));
263H = [r; e; tms];
264end
266function m = maxsq(E)
267m = max(max(E .* E));
268end
271% -------------------------------------------------------------- multigrid ---
273function W = flatten_into(W, U, off, N)
274m = N + 1;
275for j = 1:m
276 b = off + (j-1)*m;
277 for i = 1:m
278 W(b+i) = U(i,j);
279 end
280end
281end
283function ss = res_sumsq_flat(W, off, tot, N)
284m = N + 1;
285ss = 0;
286for j = 2:N
287 b = off + (j-1)*m;
288 for i = 2:N
289 r = W(tot + b + i) - (4*W(b+i) - W(b+i-1) - W(b+i+1) - W(b+i-m) - W(b+i+m));
290 ss = ss + r*r;
291 end
292end
293end
295function m = maxsq_flat(W, WE, N)
296mm = N + 1;
297m = 0;
298for k = 1:mm*mm
299 d = W(k) - WE(k);
300 if d*d > m
301 m = d*d;
302 end
303end
304end
306function W = vcycle(W, offs, ns, nlev, tot, nu1, nu2)
307for l = 1:nlev-1
308 W = mg_down(W, offs(l), offs(l+1), tot, ns(l), nu1);
309end
310W = gs_flat(W, offs(nlev), tot, ns(nlev), 8); % the coarsest grid: one unknown
311for l = nlev-1:-1:1
312 W = mg_up(W, offs(l), offs(l+1), tot, ns(l), nu2);
313end
314end
316function W = mg_down(W, off, offc, tot, n, nu)
317% Smooth, then carry the residual to the next grid and start the coarse
318% iterate at zero. The residual is formed inline, one 9-point full-weighting
319% stencil per coarse point, rather than materialised on the fine grid first.
320W = gs_flat(W, off, tot, n, nu);
321m = n + 1;
322nc = n / 2;
323mc = nc + 1;
324for jc = 2:nc
325 j = 2*jc - 1;
326 for ic = 2:nc
327 i = 2*ic - 1;
328 k = off + (j-1)*m + i;
329 r00 = W(tot+k) + W(k+1) + W(k-1) + W(k+m) + W(k-m) - 4*W(k);
330 rm0 = W(tot+k-1) + W(k) + W(k-2) + W(k-1+m) + W(k-1-m) - 4*W(k-1);
331 rp0 = W(tot+k+1) + W(k+2) + W(k) + W(k+1+m) + W(k+1-m) - 4*W(k+1);
332 r0m = W(tot+k-m) + W(k-m+1) + W(k-m-1) + W(k) + W(k-2*m) - 4*W(k-m);
333 r0p = W(tot+k+m) + W(k+m+1) + W(k+m-1) + W(k+2*m) + W(k) - 4*W(k+m);
334 rmm = W(tot+k-m-1) + W(k-m) + W(k-m-2) + W(k-1) + W(k-2*m-1) - 4*W(k-m-1);
335 rpm = W(tot+k-m+1) + W(k-m+2) + W(k-m) + W(k+1) + W(k-2*m+1) - 4*W(k-m+1);
336 rmp = W(tot+k+m-1) + W(k+m) + W(k+m-2) + W(k+2*m-1) + W(k-1) - 4*W(k+m-1);
337 rpp = W(tot+k+m+1) + W(k+m+2) + W(k+m) + W(k+2*m+1) + W(k+1) - 4*W(k+m+1);
338 s = 4*r00 + 2*(rm0 + rp0 + r0m + r0p) + (rmm + rpm + rmp + rpp);
339 % full weighting is s/16; the extra factor of 4 is because the coarse
340 % operator is written in the same scaled form and so stands for (2h)^2
341 % times the Laplacian. Leave it out and the correction comes back four
342 % times too small, which shows up as a V-cycle that reduces the
343 % residual by a factor of two instead of twenty.
344 W(tot + offc + (jc-1)*mc + ic) = s / 4;
345 end
346end
347for k = 1:mc*mc
348 W(offc + k) = 0;
349end
350end
352function W = mg_up(W, off, offc, tot, n, nu)
353% Interpolate the coarse correction, add it, smooth again. The correction is
354% built in a scratch array first: interpolating in place would mix it with the
355% fine iterate that is already there.
356m = n + 1;
357nc = n / 2;
358mc = nc + 1;
359Z = zeros(m, m);
360for jc = 1:mc
361 for ic = 1:mc
362 Z(2*ic-1, 2*jc-1) = W(offc + (jc-1)*mc + ic);
363 end
364end
365for j = 1:2:m
366 for i = 2:2:n
367 Z(i,j) = 0.5 * (Z(i-1,j) + Z(i+1,j));
368 end
369end
370for j = 2:2:n
371 for i = 1:m
372 Z(i,j) = 0.5 * (Z(i,j-1) + Z(i,j+1));
373 end
374end
375for j = 1:m
376 b = off + (j-1)*m;
377 for i = 1:m
378 W(b+i) = W(b+i) + Z(i,j);
379 end
380end
381W = gs_flat(W, off, tot, n, nu);
382end
384function W = gs_flat(W, off, tot, n, nu)
385% Gauss-Seidel on one level of the flat hierarchy: nu sweeps over the interior
386% of the level at `off`, reading its right-hand side from the second half of W.
387m = n + 1;
388for s = 1:nu
389 for j = 2:n
390 b = off + (j-1)*m;
391 for i = 2:n
392 W(b+i) = 0.25 * (W(b+i-1) + W(b+i+1) + W(b+i-m) + W(b+i+m) + W(tot+b+i));
393 end
394 end
395end
396end