m=15;                          % number of mesh points in each direction
d=2;                           % dimension of the problem
nu=1e-3;                       % regularization parameter
h=1/(m+1);                     % mesh size
A=Laplacian(m,d)/h^2;          % negative 2D discrete Laplacian
Ahat=[sparse(m^d,m^d), A, -speye(m^d); ...
      A, speye(m^d), sparse(m^d,m^d); ...
      -speye(m^d), sparse(m^d,m^d), nu*speye(m^d)]; % all-at-once operator
f=-ones(m^d,1);                % discrete f=-1
yd=ones(m^d,1);                % discrete y_d=1
fhat=[f;yd;zeros(m^d,1)];      % assemble the right-hand side vector
u0=rand(3*m^d,1);              % random initial guess
kmax=min(200,m^d);             % set max number of iterations
[u1,~,~,~,RESVEC]=gmres(@(v)(v-GS_CollectiveS(v,0*fhat,Ahat,nu,eps,1)),...
                            GS_CollectiveS(0*fhat,fhat,Ahat,nu,eps,1),...
                            [],1e-6,kmax,[],[],u0);                  
[u2,res]=GS_CollectiveS(u0,fhat,Ahat,nu,1e-14,kmax);
semilogy(0:length(RESVEC)-1,RESVEC./RESVEC(1),'b'); hold on;
semilogy(0:length(res)-1,res./res(1),'r');
axis([0 50 1e-10 1e2]); yticks([1e-10 1e-6 1e-2 1e2]);
xlabel('iterations');
ylabel('norm of the residual'); 
legend('GMRES+CS','CS Stat','location','NorthEast');
set(gca,'fontsize',18,'linewidth',2);