function [y2] = csi(M,b,i)

% Generates the ith chebyshev semi iteration vector using Jacobi iteration for a 2D problem
%
% Input : A - Matrix to be solved
%       : b - Right hand side
%       : i - number of iterations of the Cheyshev semi iteration
%
% Output : y2 - approximation to y satisfying A*y = b;

l = length(M);
I = speye(l);                     % define identity
invM = spdiags(1./diag(M), 0, l, l);

w = 4/5;                                % define the relaxation parameter

g = w*invM*b;                           % define relaxed Jacobi g
S = (I-w*invM*M);                       % define the iteration Matrix
y1 = zeros(length(b),1);
y2 = S*y1 + g;                          % x_(m+1) = Sx_m + g
mu = 5/4;                               % largest eig for this relaxation parameter

ck0 = 2;
ck1 = mu;

for j = 2:i
   ck2 = mu*ck1-0.25*ck0;
   w = mu*ck1/ck2;
   y3 = w*(S*y2+g-y1)+y1;

   y1 = y2;
   y2 = y3;
   ck0 = ck1;
   ck1 = ck2;
end

Published with MATLAB® R2019a