%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstration program for solution of 
% 
%  u_xx = f(x)
%
% using O(h^2) finite difference tecnhiques and successive overrelaxation
% (or Gauss-Seidel relaxation for omega = 1.0.
%
% This version uses red-black sub-sweeps to facilitate vectorization.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
level = 5;
nx = 2^level + 1;
x = linspace(0, 1, nx);
h = x(2) - x(1);
f = -pi^2 * sin(pi * x);
u = zeros(1, nx);
% omega = 1      -> Gauss-Seidel 
% 1 < omega <= 2 -> SOR (Successive Overrelaxation)
%
% The following value of omega is close to optimal for nx = 33
% 
% omega = 1.86875
omega = 1.0;
clf;
hold on;

% Define r(e)d and b(lac)k arrays of length nx, using logical operations
rd = rem((1:nx),2) == 1
bk = rem((1:nx),2) == 0
input('Type ENTER to continue: ');

epsilon = 1.0e-6;
du = zeros(1,nx);
dunrm = 1.0;
sweep = 0;
% Perform relaxation sweeps until convergence.
while dunrm > epsilon
   sweep = sweep + 1;
   % Note the use of vector operations in the following.

   % Save old value so that change in solution after each sweep can be computed.
   uold = u;
   % Compute CHANGES to solution on red points.
   du(2:nx-1) = (0.5 * (u(1:nx-2) + u(3:nx) - h^2*f(2:nx-1)) - u(2:nx-1)) .* rd(2:nx-1);
   % Add changes to solution---note mutiplication of change du by overrelaxation parameter.
   u = u + omega * du;
   % Compute CHANGES to solution on black points.
   du(2:nx-1) = (0.5 * (u(1:nx-2) + u(3:nx) - h^2*f(2:nx-1)) - u(2:nx-1)) .* bk(2:nx-1);
   % Add changes to solution.
   u = u + omega * du;

   % Compute change in solution during this sweep ...
   du = uold - u;
   % ... and its RMS value.
   dunrm = sqrt(sum(dot(du,du)/nx));
   if rem(sweep,50) == 0
      fprintf('Norm of solution change after %d sweeps is %g\n', ...
         sweep, dunrm);
      plot(x, u); plot(x,u,'or');
   end
end
fprintf('Used %d sweeps to converge to tolerance of %g using omega=%g\n', sweep, epsilon, omega);
