% Matlab code mgv1D.m
% For "Applied Numerical Linear Algebra",  Question 6.16
% Written by James Demmel, Apr 22, 1995
%                Modified, Jun  2, 1997
%
% Multigrid V-cycle for Poisson's equation on a 1D grid
% with zero Dirichlet boundary conditions.
% Algorithm from Brigg's ``Multigrid Tutorial''
% Include zero boundary values in arrays for ease of programming.
% Assume dimension n = 2^k + 1
% Inputs:
%   x = initial guess (n by 1 matrix with zeros on boundary)
%   b = right hand side (n by 1 matrix)
%   jac1, jac2 = number of weight Jacobi steps to do before and
%                after recursive call to mgv
% Outputs:
%   z = improved solution (n by n matrix with zeros on boundary)
%
% function z=mgv1D(x,b,jac1,jac2)
%
function z=mgv1D(x,b,jac1,jac2)
[n,m]=size(b);
if n == 3,
%  Solve 1 by 1 problem explicitly
   z=zeros(3,1);
%  z(2)=b(2)/4;
   z(2)=b(2)/2;
else
%  Perform jac1 steps of weighted Jacobi with weight w
   w=2/3;
   for i=1:jac1;
      x(2:n-1)=(1-w)*x(2:n-1) + (w/2)*( x(1:n-2) + x(3:n) + b(2:n-1));
   end,
%  Compute residual on current grid
   r=zeros(n,1);
   tmp = b(2:n-1) - ( 2*x(2:n-1) - x(1:n-2) - x(3:n) );
   r(2:n-1) = tmp;
%  Restrict residual to coarse grid
   m=(n+1)/2;
   rhat=zeros(m,1);
   rhat(2:m-1)=.5*r(3:2:n-2) + .25*(r(2:2:n-3)+r(4:2:n-1) );
%  Solve recursively with zero initial guess
%  Multiply residual by 4 for coarser grid
   xhat = mgv1D(zeros(size(rhat)),4*rhat,jac1,jac2);
%  xhat = mgv1D(zeros(size(rhat)),2*rhat,jac1,jac2);
%  Interpolate coarse solution to fine grid, and add to x
   xcor=zeros(size(x));
   xcor(3:2:n-2)=xhat(2:m-1);
   xcor(2:2:n-1)=.5*(xhat(1:m-1)+xhat(2:m));
   z=x+xcor;
%  Perform jac2 steps of weighted Jacobi with weight w
   for i=1:jac2;
      z(2:n-1)=(1-w)*z(2:n-1) + (w/2)*( z(1:n-2) + z(3:n) + b(2:n-1));
   end,
end
return
% end