% Test Full Multigrid V-Cycle code for solving Poisson's equation on a square 
% grid with zero Dirichlet boundary conditions
% Inputs: (run makemgdemo to initialize inputs for demo)
%   b = right hand side, an n=2^k+1 by n matrix with explicit zeros on boundary
%   x = initial solution guess
%   xtrue = true solution
%   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, July 10, 1993
res=[];
[n,m]=size(b);
f2=0;
figure(1)
hold off, clg
subplot(1,2,1),mesh(xtrue),title('True Solution')
subplot(1,2,2),mesh(b),title('Right Hand Side')
% print -dps MG2Dinitial.ps
% print -dgif8 MG2Dinitial.gif
disp('hit return to continue'),pause
xnew = x;
hold off, clg
%  Loop over all iterations of Full multigrid
for i=1:iter
   f1=flops;
%  Do a full multigrid cycle
   xnew=fmgv(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:n-1) - ( 4*xnew(2:n-1,2:n-1) ...
          - xnew(1:n-2,2:n-1) - xnew(3:n,2:n-1) ...
          - xnew(2:n-1,1:n-2) - xnew(2:n-1,3:n) );
   t=norm(tmp,1); res=[res;t];
   subplot(1,2,1), mesh(xnew), title(['approximate solution ',int2str(i)])
   subplot(1,2,2), mesh(xtrue-xnew), title('error')
%  plot(xnew(round(n/2),:)), title('approximate solution')
   disp('iteration, residual, error='),i,t,norm(xnew-xtrue)
   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
% This plot should be nearly horizontal, less than 1, and dependent only on
% jac1 and jac2, not the matrix dimension
figure(2)
hold off, clg
subplot(121), 
plot(res(2:iter)./res(1:iter-1),'w'), 
axis([1,iter,0,1]); 
title('residual(k+1)/residual(k)'), 
xlabel('iteration number k'), grid
% This plot should be a straight line with a negative slope
subplot(122), axis; semilogy(res,'w'),
title('residual(k)'), 
xlabel('iteration number k'), grid
% print -dps MG2DErrorSummary.ps
% print -dgif8 MG2DErrorSummary.gif
