% Matlab code testfmgv1D.m
% For "Applied Numerical Linear Algebra",  Question 6.16
% Written by James Demmel, Apr 22, 1995
%                Modified, Jun  2, 1997
%
% Test Full Multigrid V-Cycle code for solving Poisson's equation 
% on a 1D grid grid with zero Dirichlet boundary conditions
% Inputs: (run makemgdemo1D 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)
%%%%  REMOVED, 12/9/2004
%   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
%
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, clf
figure(2), hold off, clf
figure(3), hold off, clf
figure(1)
xnew = x;
hold off, clf
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,'k'),title('True Solution')
axis([0,n,xscalemin,xscalemax]); grid
subplot(2,2,2),plot(b,'k'),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
%%%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];
   figure(2)
%  subplot(1,2,1), 
   hold off, clf
   plot(xnew,'k'), 
   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), 
   figure(3)
   semilogy(abs(err),'k'), grid, hold on
   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, clf
     axes('position',[.1 .1 .8 .4]),
     semilogy(abs(err),'k'), 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),'-k'), grid, 
   elseif rem(i,3)==2,
     axes('position',[.1 .1 .8 .4]),
     semilogy(abs(err),'--k'), grid, 
   else,
     axes('position',[.1 .1 .8 .4]),
     semilogy(abs(err),'-.k'), 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, clf
% 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),'k'), 
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,'k'),
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), hold off, clf
% axes('position',[.25 .1 .5 .35]) 
for ierr=1:iter,
   if (ierr==1),
     semilogy(saveerr(:,ierr),'-k')
     hold on
   elseif (rem(ierr,3)==1),
     semilogy(saveerr(:,ierr),'-k')
   elseif (rem(ierr,3)==2),
     semilogy(saveerr(:,ierr),'--k')
   else,
     semilogy(saveerr(:,ierr),'-.k')
   end
end
for i=0:7, semilogy([1,n],10^(-i)*[1 1],':k'), end
axis([1,n,1e-7,1])
title('Error of each iteration')