-
Notifications
You must be signed in to change notification settings - Fork 1
/
fullMultigrid.m
86 lines (67 loc) · 2.13 KB
/
fullMultigrid.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
function [ x, res, mse ] = fullMultigrid( x0, A, b, up, down, iterations, xgt )
%FULLMULTIGRID Summary of this function goes here
% Detailed explanation goes here
hasgt = exist('xgt','var');
N = size(A,1);
numLevels = length(up)+1;
Alist = cell(numLevels,1);
Alist{numLevels} = A;
for i=numLevels-1:-1:1
Alist{i} = down{i}*Alist{i+1}*up{i};
end
x = x0;
res = zeros(iterations,1);
mse = zeros(iterations,1);
resh = figure;
for i=1:iterations
if( hasgt )
mse(i) = norm(x-xgt(:),2)^2/N;
end
res(i) = norm(b-A*x);
x = fullMultigrid_internal( x, Alist, b, up, down, numLevels );
fprintf(1,'%d %.1e\n',iterations-i,res(i));
set(0,'CurrentFigure',resh);
loglog(res);
drawnow;
end
end
function [x] = vcycle_internal( x0, A, b, up, down, level )%+++++++++++++++
% Pre-smoothing:
x0 = gaussSeidel(x0,A{level},b,3);
% Solve on smaller grid:
r = b - A{level}*x0;
r_small = down{level-1}*r;
if level==2
% if( size(down{1},1) < 1e4 )
% % If the system size is small, can solve exactly in memory.
% e = full(A{1})\r_small;
% else
% Otherwise, just do relaxation for a while.
e = gaussSeidel( zeros(size(r_small)), A{1}, r_small, 10 );
%fprintf(1,'%.2e\n',norm(A{1}*e-r_small));
% end
else
e = vcycle_internal( zeros(size(r_small)), A, r_small, up, down, level-1 );
end
% Fine grid correction:
x = x0 + up{level-1}*e;
% Post-smoothing:
x = gaussSeidel(x,A{level},b,1);
end%-----------------------------------------------------------------------
function [ x ] = fullMultigrid_internal( x0, A, b, up, down, level )%++++++
% I think this should be `level > 1`, but vcycle_internal below requires
% level >= 2 when called.
if( level > 2 )
r = b - A{level}*x0;
r_small = down{level-1}*r;
err_small = fullMultigrid_internal(zeros(size(r_small)), A, r_small, up, down, level-1);
err = up{level-1}*err_small;
% Right?
x = x0 + err;
% Book says this, but I'm pretty sure it's wrong.
%x = err;
else
x = x0;
end
x = vcycle_internal(x, A, b, up, down, level);
end%-----------------------------------------------------------------------