% 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