function [w,res]=GS_CollectiveS(w0,rhs,Ahat,nu,tol,kmax)
% GS_CollectiveS: Gauss-Seidel Collective Smoothing Iteration
%   [w,res]=GS_CollectiveS(v,rhs,Ahat,nu,tol,kmax) solves the optimality
%   system Ahat w=fhat  using the Gauss-Seidel Collective Smoothing 
%   method starting at the initial guess w0 up to a tolerance tol using 
%   at most kmax iterations. 
%   GS_CollectiveS returns in w the solution computed and in res the 
%   history of the norm of the residuals.

n=length(w0)/3; m=sqrt(n); h=1/(m+1);
w=w0;
y=zeros(m+2,m+2); p=y; u=y; f=y; yd=y;
f(2:end-1,2:end-1)=reshape(rhs(1:n),m,m);
yd(2:end-1,2:end-1)=reshape(rhs(n+1:2*n),m,m);
p(2:end-1,2:end-1)=reshape(w(1:n),m,m);
y(2:end-1,2:end-1)=reshape(w(n+1:2*n),m,m);
u(2:end-1,2:end-1)=reshape(w(2*n+1:3*n),m,m);
res(1)=norm(rhs-Ahat*w);
k=1;
vv=(4+h^4/(4*nu));
while k<=kmax && res(end)/res(1)>tol
  k=k+1;
  for i=2:m+1
    for j=2:m+1
      cy=(y(i-1,j)+y(i+1,j)+y(i,j-1)+y(i,j+1))+h^2*f(i,j);
      cp=(p(i-1,j)+p(i+1,j)+p(i,j-1)+p(i,j+1))+h^2*yd(i,j);
      p(i,j)=(cp-h^2*cy/4)/vv; 
      y(i,j)=(h^2*p(i,j)/nu+cy)/4; 
      u(i,j)=p(i,j)/nu;
    end
  end
  w=[p(2:end-1,2:end-1),y(2:end-1,2:end-1),u(2:end-1,2:end-1)];
  w=w(:);
  res(k)=norm(rhs-Ahat*w);
end