%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% multigrid demo on simple -D^2 u = f  u(0)=u(1)=0
%
% finest grid as 2^m+1 points (including boundaries)
%
% requires: buildA.m, coarsen.m, interp.m, jacobi.m 
%
% Fabian Waleffe, 11/15/2000 for Math 801
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m=5; j=1:2^m+1;

A=buildA(m); f=randn(length(A),1); ue=A\f; %A, f and exact solution

% First relaxation:

u=jacobi(A,f); rh=f-A*u; 

e2hs=[]; r2hs=[];

for n=m-1:-1:2,

  r2h=coarsen(rh); 
  A=buildA(n); e2h=jacobi(A,r2h); rh=r2h-A*e2h; 
  r2hs=[r2h;r2hs]; e2hs=[e2h;e2hs];       % need to save results

end

 r2h=coarsen(rh); A=buildA(1); e2h=A\r2h; % EXACT SOLUTION ON COARSEST GRID

i2=0;
for n=2:m-1,

  eh=interp(e2h);  
  i1=i2+1;i2=i2+(2^n+1); eh=e2hs(i1:i2)+eh;   % correction step

  A=buildA(n); rh=r2hs(i1:i2); e2h=jacobi(A,rh,eh);

end

eh=interp(e2h); u=u+eh; A=buildA(m); u=jacobi(A,f,u); %u=jacobi(A,f,u,3,10); 

