% 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')