% Matlab code testfmgv.m
% For "Applied Numerical Linear Algebra",  Question 6.16
% Written by James Demmel, Jul 10, 1993, modified Jun 2, 1997.
%                Simplified Jan, 2009 by Ivar Gustafsson
%
% 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
%   itmax = maximum number of full multigrid cycles to perform
%   tolres = for stopping criterion, norm of residual smaller than tolres
% Outputs:
%   xnew = improved solution
%   res = vector of residual after each iteration
%   costpi = time 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
%
res=[];
[n,m]=size(b);
f2=0;
figure(1)
hold off, clf
subplot(1,2,1),mesh(xtrue),title('True Solution')
subplot(1,2,2),mesh(b),title('Right Hand Side')
disp('hit return to continue'),pause
xnew = x;
hold off, clf
%  Loop over all iterations of Full multigrid
fprintf(' Full multigrid \n jac1=%3.0f, jac2=%3.0f, tolres=%10.2e \n',jac1,jac2, tolres)
  fprintf(' Hit return to start on a new cycle! \n \n')
  fprintf('  it  Residual  Error norm \n')
  costpi=0;
for i=1:itmax,
   t1=cputime;
%  Do a full multigrid cycle
   xnew=fmgv(xnew,b,jac1,jac2);
%  Accumulate the cpu-time
   costpi=costpi+cputime-t1;
%  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)])
%  Choose one of these depending on whether true solution is known
%   subplot(1,2,2), mesh(tmp), title('residual')
   subplot(1,2,2), mesh(xnew-xtrue), title('error')
%  plot(xnew(round(n/2),:)), title('approximate solution')
%   disp('iteration, residual, error='),i,t,norm(xnew-xtrue)
   fprintf('%4.0f %10.2e %10.2e \n',i,t,norm(xnew-xtrue)); pause
%   disp('hit return to continue'),pause
% Stop if residual small enough
   if t<tolres, break; end;
end;
iter=length(res);
% Compute number of flops per iteration per unknown
costpi=costpi/(iter*n^2);
disp('Time per iteration per unknown='),disp(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, clf
subplot(121), 
plot(res(2:iter)./res(1:iter-1),'r'), 
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,'g'),
title('residual(k)'), 
xlabel('iteration number k'), grid