function u=VCycle(A,f,u,l)
% VCYCLE performs a V-cycle of multigrid for the Laplace equation
%   u=VCycle(A,f,u,l) computes for the given linear system Au=f
%   representing a discretized Laplacian on a square domain using
%   the five point finite difference stencil a V-cycle approximation
%   of the solution starting with the guess u using l levels.
%   The resolution size m=sqrt(length(f)) must be such that m+1 can
%   be divided by 2^l. The method uses damped Jacobi, and the number
%   of pre- and post-smoothing steps nu1 and nu2 as well as the damping
%   parameter w can be set in the code.

nu1=1; nu2=1;                        % pre- and post-smoothing steps
w=2/3;                               % damping parameter for Jacobi
if l==1                              % direct solve on coarsest level
  u=A\f;
else   
  for i=1:nu1                        % pre-smoothing with relaxed Jacobi    
    u=u+w*(f-A*u)./diag(A);          
  end
  r=f-A*u;                           % compute residual
  m=sqrt(length(f));
  R=RMatrix(m);                      % compute restriction matrix R
  M=(m+1)/2-1;
  H=1/(M+1);
  AH=Laplacian(M,2)/H^2;             % compute coarse matrix
  % AH=4*R*A*R';                     % can use Galerkin instead
  e=VCycle(AH,R*r,zeros(M^2,1),l-1); % compute coarse correction
  u=u+4*R'*e;                        % use full weighting 
  for i=1:nu2                        % post-smoothing with relaxed Jacobi    
    u=u+w*(f-A*u)./diag(A);          
  end
end
