088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 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 choices are forced by numbl's
19% compiler, which types a loop variable from the expressions assigned to it and
20% drops to a slow path when a loop is not uniform. Written the obvious way this
21% file runs about twenty times slower, which would make the timings meaningless.
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);
50ss0 = res_sumsq(U, N);
51res(1) = ss0;
52err(1) = maxsq(U - UE);
53tms(1) = 0;
54stop = tol * tol * ss0;
55tacc = 0;
56for it = 1:nit
57 t1 = tic;
58 U = jacobi_sweep(U, N);
59 tacc = tacc + toc(t1);
60 res(it+1) = res_sumsq(U, N);
61 err(it+1) = maxsq(U - UE);
62 tms(it+1) = tacc * 1000;
63 if res(it+1) < stop
64 break;
65 end
66end
67% The squares are handed back as they are; the page takes the square roots.
68% Doing it here would mean either a sqrt inside the loop, which widens the
69% inferred type of the history to complex and costs a factor of twenty, or a
70% call after it, which does the same thing to the loop that filled the arrays.
71H = [res / ss0; err; tms];
72end
74function H = loop_sor(N, U, UE, nit, tol, omega)
75% omega = 1 is Gauss-Seidel: new values are used as soon as they are available,
76% which costs nothing and halves the sweep count. Over-relaxing with the
77% optimal omega = 2/(1 + sin(pi h)) changes the rate from 1 - O(h^2) to
78% 1 - O(h). That optimal value is known only because the eigenvalues of this
79% particular matrix are known.
80res = NaN(1, nit+1);
81err = NaN(1, nit+1);
82tms = NaN(1, nit+1);
83ss0 = res_sumsq(U, N);
84res(1) = ss0;
85err(1) = maxsq(U - UE);
86tms(1) = 0;
87stop = tol * tol * ss0;
88tacc = 0;
89for it = 1:nit
90 t1 = tic;
91 U = sor_sweep(U, N, omega);
92 tacc = tacc + toc(t1);
93 res(it+1) = res_sumsq(U, N);
94 err(it+1) = maxsq(U - UE);
95 tms(it+1) = tacc * 1000;
96 if res(it+1) < stop
97 break;
98 end
99end
100% The squares are handed back as they are; the page takes the square roots.
101% Doing it here would mean either a sqrt inside the loop, which widens the
102% inferred type of the history to complex and costs a factor of twenty, or a
103% call after it, which does the same thing to the loop that filled the arrays.
104H = [res / ss0; err; tms];
105end
107function H = loop_cg(N, U, UE, nit, tol)
108% Conjugate gradients. A is symmetric positive definite, so CG applies; it
109% needs one operator application per iteration and four vectors of storage.
110% Its convergence is governed by sqrt(kappa) rather than kappa, which for this
111% matrix means O(N) iterations, like optimal SOR, but without having to know
112% anything about the spectrum in advance.
113res = NaN(1, nit+1);
114err = NaN(1, nit+1);
115tms = NaN(1, nit+1);
116R = neg_residual(U, N); % the right-hand side is zero at interior points
117P = R;
118rs = dot2(R, R);
119ss0 = rs;
120res(1) = ss0;
121err(1) = maxsq(U - UE);
122tms(1) = 0;
123stop = tol * tol * ss0;
124tacc = 0;
125for it = 1:nit
126 t1 = tic;
127 AP = applyA(P, N);
128 al = rs / dot2(P, AP);
129 U = U + al * P;
130 R = R - al * AP;
131 rsn = dot2(R, R);
132 P = R + (rsn / rs) * P;
133 rs = rsn;
134 tacc = tacc + toc(t1);
135 res(it+1) = rsn;
136 err(it+1) = maxsq(U - UE);
137 tms(it+1) = tacc * 1000;
138 if res(it+1) < stop
139 break;
140 end
141end
142% The squares are handed back as they are; the page takes the square roots.
143% Doing it here would mean either a sqrt inside the loop, which widens the
144% inferred type of the history to complex and costs a factor of twenty, or a
145% call after it, which does the same thing to the loop that filled the arrays.
146H = [res / ss0; err; tms];
147end
149function H = loop_mg(N, U, UE, nit, tol)
150% Geometric multigrid. Gauss-Seidel is a poor solver but an excellent smoother:
151% two sweeps remove the oscillatory part of the error and leave something
152% smooth, which is well represented on a grid twice as coarse and much cheaper
153% to correct there. Recursing on that idea gives a V-cycle whose cost is
154% proportional to the number of unknowns and whose error reduction per cycle is
155% a constant independent of N.
156nlev = 1;
157while mod(N, 2^nlev) == 0 && N / 2^nlev >= 2
158 nlev = nlev + 1;
159end
160ns = zeros(nlev, 1);
161offs = zeros(nlev, 1);
162tot = 0;
163nn = N;
164for l = 1:nlev
165 ns(l) = nn;
166 offs(l) = tot;
167 tot = tot + (nn+1)^2;
168 nn = nn / 2;
169end
171% one flat vector holds every level: the iterate in the first half, the
172% right-hand side in the second
173W = zeros(2*tot, 1);
174W = flatten_into(W, U, 0, N);
175WE = flatten_into(zeros((N+1)^2, 1), UE, 0, N);
177res = NaN(1, nit+1);
178err = NaN(1, nit+1);
179tms = NaN(1, nit+1);
180ss0 = res_sumsq_flat(W, 0, tot, N);
181res(1) = ss0;
182err(1) = maxsq_flat(W, WE, N);
183tms(1) = 0;
184stop = tol * tol * ss0;
185tacc = 0;
186for it = 1:nit
187 t1 = tic;
188 W = vcycle(W, offs, ns, nlev, tot, 2, 2);
189 tacc = tacc + toc(t1);
190 res(it+1) = res_sumsq_flat(W, 0, tot, N);
191 err(it+1) = maxsq_flat(W, WE, N);
192 tms(it+1) = tacc * 1000;
193 if res(it+1) < stop
194 break;
195 end
196end
197% The squares are handed back as they are; the page takes the square roots.
198% Doing it here would mean either a sqrt inside the loop, which widens the
199% inferred type of the history to complex and costs a factor of twenty, or a
200% call after it, which does the same thing to the loop that filled the arrays.
201H = [res / ss0; err; tms];
202end
204% ---------------------------------------------------------------- kernels ---
206function U = jacobi_sweep(U, N)
207V = U;
208for j = 2:N
209 for i = 2:N
210 U(i,j) = 0.25 * (V(i-1,j) + V(i+1,j) + V(i,j-1) + V(i,j+1));
211 end
212end
213end
215function U = sor_sweep(U, N, omega)
216for j = 2:N
217 for i = 2:N
218 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));
219 end
220end
221end
223function Y = applyA(U, N)
224Y = zeros(N+1, N+1);
225for j = 2:N
226 for i = 2:N
227 Y(i,j) = 4*U(i,j) - U(i-1,j) - U(i+1,j) - U(i,j-1) - U(i,j+1);
228 end
229end
230end
232function Y = neg_residual(U, N)
233Y = zeros(N+1, N+1);
234for j = 2:N
235 for i = 2:N
236 Y(i,j) = -(4*U(i,j) - U(i-1,j) - U(i+1,j) - U(i,j-1) - U(i,j+1));
237 end
238end
239end
241function s = dot2(A, B)
242s = 0;
243n = numel(A);
244for k = 1:n
245 s = s + A(k) * B(k);
246end
247end
249function ss = res_sumsq(U, N)
250% The sum of squares, not the norm: a square root taken here would widen the
251% inferred type of the accumulator and slow this loop by a factor of 20.
252ss = 0;
253for j = 2:N
254 for i = 2:N
255 r = 4*U(i,j) - U(i-1,j) - U(i+1,j) - U(i,j-1) - U(i,j+1);
256 ss = ss + r*r;
257 end
258end
259end
261function m = maxsq(E)
262m = max(max(E .* E));
263end
266% -------------------------------------------------------------- multigrid ---
268function W = flatten_into(W, U, off, N)
269m = N + 1;
270for j = 1:m
271 b = off + (j-1)*m;
272 for i = 1:m
273 W(b+i) = U(i,j);
274 end
275end
276end
278function ss = res_sumsq_flat(W, off, tot, N)
279m = N + 1;
280ss = 0;
281for j = 2:N
282 b = off + (j-1)*m;
283 for i = 2:N
284 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));
285 ss = ss + r*r;
286 end
287end
288end
290function m = maxsq_flat(W, WE, N)
291mm = N + 1;
292m = 0;
293for k = 1:mm*mm
294 d = W(k) - WE(k);
295 if d*d > m
296 m = d*d;
297 end
298end
299end
301function W = vcycle(W, offs, ns, nlev, tot, nu1, nu2)
302for l = 1:nlev-1
303 W = mg_down(W, offs(l), offs(l+1), tot, ns(l), nu1);
304end
305W = gs_flat(W, offs(nlev), tot, ns(nlev), 8); % the coarsest grid: one unknown
306for l = nlev-1:-1:1
307 W = mg_up(W, offs(l), offs(l+1), tot, ns(l), nu2);
308end
309end
311function W = mg_down(W, off, offc, tot, n, nu)
312% Smooth, then carry the residual to the next grid and start the coarse
313% iterate at zero.
314%
315% The Gauss-Seidel sweeps are written out here rather than delegated to
316% gs_flat, and the residual is formed inline rather than through a helper, for
317% the same compiler reason as everywhere else in this file: a loop that writes
318% W by index runs fast only if W did not arrive from a function call. Calling
319% gs_flat first and then entering the restriction loop costs a factor of 250 on
320% this function alone, which was measured, not guessed.
321m = n + 1;
322for s = 1:nu
323 for j = 2:n
324 b = off + (j-1)*m;
325 for i = 2:n
326 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));
327 end
328 end
329end
331nc = n / 2;
332mc = nc + 1;
333for jc = 2:nc
334 j = 2*jc - 1;
335 for ic = 2:nc
336 i = 2*ic - 1;
337 k = off + (j-1)*m + i;
338 r00 = W(tot+k) + W(k+1) + W(k-1) + W(k+m) + W(k-m) - 4*W(k);
339 rm0 = W(tot+k-1) + W(k) + W(k-2) + W(k-1+m) + W(k-1-m) - 4*W(k-1);
340 rp0 = W(tot+k+1) + W(k+2) + W(k) + W(k+1+m) + W(k+1-m) - 4*W(k+1);
341 r0m = W(tot+k-m) + W(k-m+1) + W(k-m-1) + W(k) + W(k-2*m) - 4*W(k-m);
342 r0p = W(tot+k+m) + W(k+m+1) + W(k+m-1) + W(k+2*m) + W(k) - 4*W(k+m);
343 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);
344 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);
345 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);
346 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);
347 s = 4*r00 + 2*(rm0 + rp0 + r0m + r0p) + (rmm + rpm + rmp + rpp);
348 % full weighting is s/16; the extra factor of 4 is because the coarse
349 % operator is written in the same scaled form and so stands for (2h)^2
350 % times the Laplacian. Leave it out and the correction comes back four
351 % times too small, which shows up as a V-cycle that reduces the
352 % residual by a factor of two instead of twenty.
353 W(tot + offc + (jc-1)*mc + ic) = s / 4;
354 end
355end
356for k = 1:mc*mc
357 W(offc + k) = 0;
358end
359end
361function W = mg_up(W, off, offc, tot, n, nu)
362% Interpolate the coarse correction, add it, smooth again. The correction is
363% built in a scratch array first: interpolating in place would mix it with the
364% fine iterate that is already there.
365m = n + 1;
366nc = n / 2;
367mc = nc + 1;
368Z = zeros(m, m);
369for jc = 1:mc
370 for ic = 1:mc
371 Z(2*ic-1, 2*jc-1) = W(offc + (jc-1)*mc + ic);
372 end
373end
374for j = 1:2:m
375 for i = 2:2:n
376 Z(i,j) = 0.5 * (Z(i-1,j) + Z(i+1,j));
377 end
378end
379for j = 2:2:n
380 for i = 1:m
381 Z(i,j) = 0.5 * (Z(i,j-1) + Z(i,j+1));
382 end
383end
384for j = 1:m
385 b = off + (j-1)*m;
386 for i = 1:m
387 W(b+i) = W(b+i) + Z(i,j);
388 end
389end
390for s = 1:nu
391 for j = 2:n
392 b = off + (j-1)*m;
393 for i = 2:n
394 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));
395 end
396 end
397end
398end
400function W = gs_flat(W, off, tot, n, nu)
401% Gauss-Seidel on one level of the flat hierarchy. Used only for the coarsest
402% grid, where there is a single unknown and eight sweeps solve it exactly; the
403% sweeps on the other levels are written out inside mg_down and mg_up.
404m = n + 1;
405for s = 1:nu
406 for j = 2:n
407 b = off + (j-1)*m;
408 for i = 2:n
409 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));
410 end
411 end
412end
413end