%Using PCG to solve Tx=b with different circulant preconditioners
clear
maxit=10000;tol=1e-7;
fprintf('n\t None\t\t Strang \t T. Chan\t R. Chan\t Tridiag\n');
for n=1e3*(1:5)%matrix size
%First row of T
t(1)= pi^2/3; t(2:n) = 2*(-1).^(1:n-1)./((1:n-1).^2); %Sinc matrix
t1 = conj(t(n:-1:2)); % last column of T
gev = real(fft([t 0 t1].'));%the eigenvalues of the embeding circulant matrix
b=ones(n,1);
Proj=@(z) z(1:length(z)/2);%Take first half of the long vector
Afun=@(v) Proj(ifft(fft([v;zeros(n,1)]).*gev));%Matrix-free product T*v

% no preconditioner
[x,flag,relres,iter] = mypcgfun(Afun,b,tol,maxit);

% Strang's preconditioner.
evS = fft([t(1:n/2) 0 conj(t(n/2:-1:2))].');
Mfun=@(r) ifft(fft(r)./evS); %preconditioning step
[xS,flagS,relresS,iterS] = mypcgfun(Afun,b,tol,maxit,Mfun);

% T. Chan's preconditioner
coef = 1/n:1/n:1-1/n;
evT = fft([t(1) (1-coef).*t(2:n)+coef.*t1])';
Mfun=@(r) ifft(fft(r)./evT); %preconditioning step
[xT,flagT,relresT,iterT] = mypcgfun(Afun,b,tol,maxit,Mfun);

% R. Chan's preconditioner
evR = fft([t(1) t(2:n)+t1].');
Mfun=@(r) ifft(fft(r)./evR); %preconditioning step
[xR,flagR,relresR,iterR] = mypcgfun(Afun,b,tol,maxit,Mfun);

% tridiagonal preconditioner
Ptr=gallery('tridiag',n,-1,2,-1);
Mfun=@(r) Ptr\r; %tridiagonal, Thomas algorithm solving O(n) operations
[xTr,flagTr,relresTr,iterTr] = mypcgfun(Afun,b,tol,maxit,Mfun);

fprintf('%d\t %d\t\t%d\t\t%d\t\t%d\t\t %d\n',n,iter,iterS,iterT,iterR,iterTr);

%[norm(x-xT) norm(x-xS) norm(x-xR) norm(x-xTr)]
end
figure(1)
loglog(relres,'--');hold on
loglog(relresT,'-x');hold on
loglog(relresS,'-s');hold on
loglog(relresR,'-o');hold on
loglog(relresTr,'-*');hold on

legend('None','T. Chan','Strang','R. Chan','Tridiagonal')
xlabel('Iterations')
ylabel('Residual norm')