function z=mgv1D(x,b,jac1,jac2)
% 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)
%
% Written by J. Demmel, UC Berkeley, April 22, 1995
[n,m]=size(b);
if n == 3,
% Solve 1 by 1 problem explicitly
z=zeros(3,1);
z(2)=b(2)/4;
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);
% 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