function SecondFD1dMG()
%solving -u''(x)=f(x) with given u(xmin) u(xmax)
%including boundary into system for easier interpolation/restriction
ffun=@(x)-exp(x); ufun=@(x)exp(x);
tol=1e-10;maxit=20;
for fineL=5:10
tic
%Setup operators
mg=[];
for Level=fineL:-1:2
[mg(Level).A]=genMatrix(2^Level); %generate system at each level
end
A=mg(fineL).A;
%Setup right-hand-side
xmin=0; xmax=pi; N=2^fineL; h=(xmax-xmin)/N;
xgrid=(0:N)'*h; %including boundary nodes
b=ffun(xgrid); %rhs
b(1)=ufun(xmin); b(end)=ufun(xmax);%boundary condition
truex=ufun(xgrid); %true solution

x=rand(size(b)); %initial guess
it=1;r0=norm(b-A*x);rk=r0;
err=[];res=[];
%Begin Multigrid Algorithm: V-cycle
while((it<=maxit)&&(rk>=r0*tol))
x=mg_iter(mg,x,b,fineL,2,1,'JAC'); %smoother: JAC, GS, RB-GS
rk=norm(b-A*x);
err(it)=max(abs(x-truex)); %max error
res(it)=rk/r0; %relative residual
it=it+1;
end
tcpu=toc;
fprintf('m=%d, V-cycle iter=%d, Error = %10.3e, Residual=%10.3e, CPU=%1.3f \n',N,it,err(end),res(end),tcpu)

figure(1)
subplot(1,2,1)
semilogy([1:it-1],res,'-o');xlabel('iteration'); ylabel('residual norms');
hold on

subplot(1,2,2)
semilogy([1:it-1],err,'-o');xlabel('iteration'); ylabel('error norms');
hold on
end

end

function [A]=genMatrix(N)
xmin=0; xmax=pi; h=(xmax-xmin)/N; e = ones(N-1,1);
%include boundary nodes
A=[sparse(1,1,1,1,N+1);
spdiags([-e 2*e -e]/h^2, 0:2, N-1, N+1);
sparse(1,N+1,1,1,N+1)];
end

function [x]=mg_iter(mg,x0,b,level,pre,post,smoother)
%Multigrid V-cycle Recursive code
%fprintf('+++++++++++++++++++Level[%d] Begin+++++++++++++++++++++\n',level);
A=mg(level).A;
if(level==2) %Coarest level
x=A\b;%direct solver, a very small system
else
% presmooth
%fprintf('Level[%d],Before Presmooth res=%1.2e\n',level,norm(b-A*x0));
x = mg_smooth(A,x0,b,pre,smoother);
%fprintf('Level[%d],After Presmooth res=%1.2e\n',level,norm(b-A*x));
% Restrict residual
r = b - A*x;
rc = mg_restrict(r);
% coarse grid correction
cc = mg_iter(mg,zeros(size(rc)),rc,level-1,pre,post,smoother);
x = x + mg_interp(cc);
% postsmooth
%fprintf('Level[%d],Before Postsmooth res=%1.2e\n',level,norm(b-A*x));
x = mg_smooth(A,x,b,post,smoother);
%fprintf('Level[%d],After Postsmooth res=%1.2e\n',level,norm(b-A*x));
end
%fprintf('+++++++++++++++++++Level[%d] End+++++++++++++++++++++\n',level);
end

function [x]=mg_smooth(A,x0,b,nv,smoother)
%smoother iterations (also called relaxations)
x=x0;
%Jacobi Smoother
switch smoother
case 'JAC'
w=2/3;%Weighted Jacobi Smoother
dinv = 1./diag(A);
for k = 1:nv
x=x + w*dinv.*(b - A*x);
end
case 'GS' %Gauss-Seidel Smoother
L = tril(A);
for k = 1:nv
x=x + L\(b - A*x);
end
case 'RB-GS' %Red-Black G-S Smoother
dinv = 1./diag(A);
for k = 1:nv
x(1:2:end)=x(1:2:end) + dinv(1:2:end).*(b(1:2:end) - A(1:2:end,:)*x);
x(2:2:end)=x(2:2:end) + dinv(2:2:end).*(b(2:2:end) - A(2:2:end,:)*x);
end
end
end

function [rc]=mg_restrict(rf)
N=length(rf)-1;
%Injection
%rc=[rf(1);rf(3:2:N-1);rf(N+1)];
%Full weight, usually better
rc=[rf(1);(rf(2:2:N-2)+2*rf(3:2:N-1)+rf(4:2:N))/4;rf(N+1)];
end

function [rf]=mg_interp(rc)
N=length(rc)-1;
rf=zeros(2*N+1,1);
rf(1:2:2*N+1)=rc;rf(2:2:2*N)=0.5*(rf(1:2:2*N-1)+rf(3:2:2*N+1));
end