% Test Full Multigrid V-Cycle code for solving Poisson's equation on a 1D grid
% grid with zero Dirichlet boundary conditions
% Inputs: (run makemgdemo to initialize inputs for demo)
%   b = right hand side, an n=2^k+1 by 1 matrix with explicit zeros on boundary
%   x = initial solution guess
%   jac1, jac2 = number of weight Jacobi steps to do before and
%                after recursive call to mgv
%   iter = number of full multigrid cycles to perform
% Outputs:
%   xnew = improved solution
%   res = vector of residual after each iteration
%   costpi = number of flops per iteration per unknown
%            (should be constant independent of number of unknowns)
%   plot of approximate solution after each iteration
%   residual after each iteration
%   plot of residual(i+1)/residual(i), which should be constant independent of
%     right hand side and dimension
%   plot of residual versus iteration numbers
%
% Written by J. Demmel, UC Berkeley, April 22, 1995
saveerr=[];
res=[];
[n,m]=size(b);
tt=2*eye(n-2)-diag(ones(n-3,1),1)-diag(ones(n-3,1),-1);
f2=0;
figure(1), hold off, clg
figure(2), hold off, clg
figure(3), hold off, clg
figure(1)
xnew = x;
hold off, clg
scaler = 1;
bscalemin=floor(min(b)/scaler)*scaler;
bscalemax=ceil(max(b)/scaler)*scaler;
xtrue=[0;tt\b(2:n-1);0];
xscalemin=floor(min(xtrue)/scaler)*scaler;
xscalemax=ceil(max(xtrue)/scaler)*scaler;
subplot(2,2,1),plot(xtrue,'w'),title('True Solution')
axis([0,n,xscalemin,xscalemax]); grid
subplot(2,2,2),plot(b,'w'),title('Right Hand Side')
axis([0,n,bscalemin,bscalemax]); grid
% print -dps MG1Dinitial.ps
% print -dgif8 MG1Dinitial.gif
% disp('hit return to continue'), pause,
%  Loop over all iterations of Full multigrid
for i=1:iter
   figure(2)
   f1=flops;
%  Do a full multigrid cycle
   xnew=fmgv1D(xnew,b,jac1,jac2);
%  Accumulate the number of floating point operations
   f2=(flops-f1)+f2;
%  Compute and save residual
   tmp = b(2:n-1) - ( 2*xnew(2:n-1) - xnew(1:n-2) - xnew(3:n) );
   t=norm(tmp,1); res=[res;t];
   hold off, clg
   subplot(1,2,1), plot(xnew,'w'), 
   grid, title(['approximate solution ',int2str(i)])
   grid, axis([0,n,xscalemin,xscalemax]); grid
   err=xnew-[0;tt\b(2:n-1);0];
   subplot(1,2,2), semilogy(abs(err),'w'), grid
   title(['error, norm = ', num2str(norm(err))])
   grid; axis([0,n,1e-7,scaler]),grid;
%  plot(xnew(round(n/2),:)), title('approximate solution')
   disp('iteration, residual='),i,t, 
%  disp('hit return to continue'), pause
%  figure(3)
   saveerr=[saveerr,abs(err)];
   if (1==0)
   if i==1,
     hold off, clg
     axes('position',[.1 .1 .8 .4]),
     semilogy(abs(err),'w'), grid, title('error after each iteration m')
     grid; axis([0,n,1e-7,scaler]),grid;
     hold on
   elseif rem(i,3)==1,
     axes('position',[.1 .1 .8 .4]),
     semilogy(abs(err),'-w'), grid, 
   elseif rem(i,3)==2,
     axes('position',[.1 .1 .8 .4]),
     semilogy(abs(err),'--w'), grid, 
   else,
     axes('position',[.1 .1 .8 .4]),
     semilogy(abs(err),'-.w'), grid, 
   end
   end
   end
%  disp('hit return to continue'), pause
end
% Compute number of flops per iteration per unknown
costpi=f2/(iter*n^2);
disp('flops per iteration per unknown='),costpi
figure(2)
% print -dps MG1DErrorPerStep.ps
% print -dgif8 MG1DErrorPerStep.gif
figure(1)
% hold off, clg
% This plot should be nearly horizontal, less than 1, and dependent only on
% jac1 and jac2, not the matrix dimension
subplot(223), 
% axes('position',[.1,.6,.3,.3])
plot(res(2:iter)./res(1:iter-1),'w'), 
axis([1,iter,0,1]); 
title('norm(res(m+1))/norm(res(m))'), 
xlabel('iteration number m'), grid
% This plot should be a straight line with a negative slope
subplot(224), 
% axes('position',[.5,.6,.3,.3])
semilogy(res,'w'),
axis([1,iter,1e-7,1]);
title('norm(res(m))'), 
xlabel('iteration number m'), grid
% print -dps MG1DErrorSummary.ps
% print -dgif8 MG1DErrorSummary.gif
figure(3)
% axes('position',[.25 .1 .5 .35]) 
for ierr=1:iter,
   if (ierr==1),
     semilogy(saveerr(:,ierr),'-w')
     hold on
   elseif (rem(ierr,3)==1),
     semilogy(saveerr(:,ierr),'-w')
   elseif (rem(ierr,3)==2),
     semilogy(saveerr(:,ierr),'--w')
   else,
     semilogy(saveerr(:,ierr),'-.w')
   end
end
for i=0:7, semilogy([1,n],10^(-i)*[1 1],':w'), end
axis([1,n,1e-7,1])
title('Error of each iteration')
