function Smoothers4multigrid()
% this program compares performances of three different smoothers for
% multigrid:
% method 1: Red-Black Gauss-Seidel
% method 2: weighted Jacobi, omega = 4/5
% method 3: SOR, omega = 2/(1 + sin(pi*h))

global n r b nr nb nr1 nb1 f h2 h I Ip Im
figure(1);
clf;
hold on;

col(1) = 'k';
col(2) = 'r';
col(3) = 'b';

p = 4;
n = 2^p + 1;
n1 = n - 1;
t = linspace(0,1,n);
[x,y] = meshgrid(t,t);
f = 2*(x.^2 - x + y.^2 - y);
h = 1/n1;
h2 = h^2;
u1 = rand(n - 2);
uinit = zeros(n);
uinit(2 : n1,2:n1) = u1;
u_exact = x.*y.*(1 - x).*(1 - y);

% for Red-Black Gauss-Seidel  and SOR methods
r = 1 : 2 : n;
nr = length(r);
nr1 = nr - 1;
b = 2 : 2 : n1;
nb = length(b);
nb1 = nb - 1;

% for weighted Jacobi method
I = [2 : n1];
Ip = [3 : n];
Im = [1 : n - 2];



itermax = 1000;
for method = 1 : 3
    u = uinit;
    tic
    for k = 1 : itermax
        if method == 1
            u = RBGS(u);            
        else if method == 2
                u = WJac(u);
            else if method == 3
                    u = SOR(u);
                end
            end
        end

        % calculate error
        e(k) = max(max(abs(u - u_exact)));
        if k < 10 & mod(k,2) == 1
            err = u - u_exact;
            figure(2);
            subplot(5,3,3*(k-1)/2 + method)
            mesh(err)
            if method == 1
                str = 'Method: RBGS';
            else if method == 2
                    str = 'Method: WJac';
                else if method == 3
                        str = 'Method: SOR';
                    end
                end
            end
            title(strcat(str,sprintf(', iter # %d\n',k)));
        end
    end
    toc
    figure(1);
    plot(e,'Linewidth',2,'color',col(method))
    set(gca,'Yscale','log');
end
end

%%
function u = RBGS(u);
global n r b nr nb nr1 nb1 f h2
    % update all red points
    u(r(2 : nr1),r(2 : nr1)) = 0.25*(u(b(1 : nb1),r(2 : nr1)) + ...
        u(b(2 : nb),r(2 : nr1)) + u(r(2 : nr1),b(1 : nb1)) + ...
        u(r(2 : nr1),b(2 : nb)) - h2*f(r(2 : nr1),r(2 : nr1)));
    
    u(b(1 : nb),b(1 : nb)) = 0.25*(u(r(1 : nr1),b(1 : nb)) + ...
        u(r(2 : nr),b(1 : nb)) + u(b(1 : nb),r(1 : nr1)) + ...
        u(b(1 : nb),r(2 : nr)) - h2*f(b(1 : nb),b(1 : nb)));
    
    % update all black points   
    u(b(1 : nb),r(2 : nr1)) = 0.25*(u(r(1 : nr1),r(2 : nr1)) + ...
        u(r(2 : nr),r(2 : nr1)) + u(b(1 : nb),b(1 : nb1)) + ...
        u(b(1 : nb),b(2 : nb)) - h2*f(b(1 : nb),r(2 : nr1)));
    
    u(r(2 : nr1),b(1 : nb)) = 0.25*(u(b(1 : nb1),b(1 : nb)) + ...
        u(b(2 : nb),b(1 : nb)) + u(r(2 : nr1),r(1 : nr1)) + ...
        u(r(2 : nr1),r(2 : nr)) - h2*f(r(2 : nr1),b(1 : nb)));
end

%%
function u = WJac(u)
global n f h2 I Ip Im
om = 4/5;
unew = u;
unew(I,I) = (1 - om)*u(I,I) + om*0.25*(u(I,Ip) + u(I,Im) + u(Ip,I) + u(Im,I) - h2*f(I,I));
u = unew;
end

%%
function u = SOR(u);
global n r b nr nb nr1 nb1 f h2 h
 om = 2/( 1 + sin(pi*h));
    % update all red points
    u(r(2 : nr1),r(2 : nr1)) = (1 - om)*u(r(2 : nr1),r(2 : nr1)) + ...
        om*0.25*(u(b(1 : nb1),r(2 : nr1)) + ...
        u(b(2 : nb),r(2 : nr1)) + u(r(2 : nr1),b(1 : nb1)) + ...
        u(r(2 : nr1),b(2 : nb)) - h2*f(r(2 : nr1),r(2 : nr1)));
    
    u(b(1 : nb),b(1 : nb)) = (1 - om)*u(b(1 : nb),b(1 : nb)) + ...
        om*0.25*(u(r(1 : nr1),b(1 : nb)) + ...
        u(r(2 : nr),b(1 : nb)) + u(b(1 : nb),r(1 : nr1)) + ...
        u(b(1 : nb),r(2 : nr)) - h2*f(b(1 : nb),b(1 : nb)));
    
    % update all black points   
    u(b(1 : nb),r(2 : nr1)) = (1 - om)*u(b(1 : nb),r(2 : nr1)) + ...
        om*0.25*(u(r(1 : nr1),r(2 : nr1)) + ...
        u(r(2 : nr),r(2 : nr1)) + u(b(1 : nb),b(1 : nb1)) + ...
        u(b(1 : nb),b(2 : nb)) - h2*f(b(1 : nb),r(2 : nr1)));
    
    u(r(2 : nr1),b(1 : nb)) = (1 - om)*u(r(2 : nr1),b(1 : nb)) + ...
        om*0.25*(u(b(1 : nb1),b(1 : nb)) + ...
        u(b(2 : nb),b(1 : nb)) + u(r(2 : nr1),r(1 : nr1)) + ...
        u(r(2 : nr1),r(2 : nr)) - h2*f(r(2 : nr1),b(1 : nb)));
end


