function u=VCycleOC(A,f,nu,u,l)
% VCYCLEOC performs a V-cycle of multigrid for the Laplace equation
%   u=VCycleOC(A,f,nu,u,l) computes for the given linear system Au=f
%   representing a an optimality system of a linear quadratic elliptic
%   optimal control problem (with an l2 regularization parameter nu) 
%   defined on a square domain a V-cycle approximation of the solution 
%   starting with the guess u using l levels.
%   The resolution size m=sqrt(length(f)/2) must be such that m+1 can
%   be divided by 2^l. The method uses a Gauss-Seidel collective smoothing
%   method, and the number of pre- and post-smoothing steps nu1 and nu2 
%   can be set in the code.

nu1=1; nu2=1;                        % pre- and post-smoothing steps
if l==1                              % direct solve on coarsest level
  u=A\f;
else   
  u=GS_CollectiveS(u,f,A,nu,eps,nu1);% pre-smoothing
  r=f-A*u;                           % compute residual
  m=sqrt(length(f)/3);
  R=RMatrix(m);                      % compute restriction matrix R
  Rt=blkdiag(R,R,R);
  M=(m+1)/2-1;
  AH=4*Rt*A*Rt';                     % compute (Galerkin) coarse matrix
  e=VCycleOC(AH,Rt*r,nu,zeros(3*M^2,1),l-1);% compute coarse correction
  u=u+4*Rt'*e;                       % use full weighting 
  u=GS_CollectiveS(u,f,A,nu,eps,nu2);% post-smoothing
end
