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
Mt = blkdiag(A*A,speye(m^d),nu*speye(m^d)); % block-diag precond Mtilde
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
[u,~,~,~,rv]=gmres(Ahat,fhat,[],1e-6,kmax,Mt,[],u0);
semilogy(0:length(rv)-1,rv./rv(1),'b-');
axis([0 50 1e-10 1]); yticks([1e-10 1e-8 1e-6 1e-4 1e-2 1]);
xlabel('iterations');
ylabel('norms of the residuals'); 
set(gca,'fontsize',18,'linewidth',2);

