function qcgp = pcgssor ( N, omega, rcgp );
% PCGSSOR
%
% This m-file is part of the material of the lectures on
%   Numerical Mathematics
% at the Math Department of Halle University, Germany.
%
% Author :      Prof. Dr. M. Arnold, arnold@mathematik.uni-halle.de
% Version of :  Jul 13, 2005
%
% SSOR preconditioner for CG method applied to model problem
%
% Parameters:
%   N      (input)  : number of grid points in FD approximation
%   omega  (input)  : (optimal) relaxation parameter
%   rcgp   (input)  : residual vector
%   qcgp   (output) : solution of system of linear equations
%
% Example:
%   see cgiter.m
%
%   Reference: Algorithm 8.7.10 of
%   Stoer, J.; Bulirsch, R.: Numerische Mathematik 2, Dritte Auflage.
%   Springer-Verlag Berlin Heidelberg New York, 1990.

% -> initialization
h = 1/N;

% -> rewrite right hand side in matrix form
Uaux = zeros ( N+1, N+1 );
% for i=1:N-1,
%     Uaux(1+(1:N-1),1+i) = rcgp((i-1)*(N-1)+(1:N-1));
% end;

% -> SSOR step
for i=1:N-1,  for j=1:N-1,
    Uaux(1+i,1+j) = ( 1 - omega ) * Uaux(1+i,1+j)              ...
                    + omega * ( rcgp((i-1)*(N-1)+j)            ...
                    + Uaux(1+(i-1),1+j) + Uaux(1+(i+1),1+j)    ...
                    + Uaux(1+i,1+(j-1)) + Uaux(1+i,1+(j+1)) ) / 4;
end;  end;
for i=N-1:-1:1,  for j=N-1:-1:1,
    Uaux(1+i,1+j) = ( 1 - omega ) * Uaux(1+i,1+j)              ...
                    + omega * ( rcgp((i-1)*(N-1)+j)            ...
                    + Uaux(1+(i-1),1+j) + Uaux(1+(i+1),1+j)    ...
                    + Uaux(1+i,1+(j-1)) + Uaux(1+i,1+(j+1)) ) / 4;
end;  end;

% -> rewrite solution in vector form
qcgp = zeros ( (N-1)^2, 1 );
for i=1:N-1,
    qcgp((i-1)*(N-1)+(1:N-1)) = Uaux(1+(1:N-1),1+i);
end;
