function R=RMatrix(m)
% RMATRIX construct restriction matrix of size M^2 x m^2 
%   where m must have the form m=2^l-1 and hence M is M=2^(l-1)-1
%   R=RMatrix(m) constructs a restriction matrix of size M^2 x m^2
%   to be used in a multigrid code.

e=ones(m^2,1);                      % diagonal entries
e1=repmat([ones(m-1,1);0],m,1);
e2=ones(m^2,1);
B=zeros(m^2,9);
B(1:length(e1),1)=e1/16;
B(1:length(e2),2)=e2/8;        
B(1:length(e1),3)=flipud(e1/16);  
B(1:length(e1),4)=e1/8;
B(1:length(e),5)=e/4;
B(1:length(e1),6)=flipud(e1/8);     % flipud for spdiags convention 
B(1:length(e1),7)=e1/16;
B(1:length(e2),8)=flipud(e2/8);
B(1:length(e1),9)=flipud(e1/16);

R=spdiags(B,[-m-1 -m -m+1 -1 0 1 m-1 m m+1],m^2,m^2); 
G=numgrid('S',m+2);                 % extract even gridpoints
id=G(3:2:m,3:2:m);
R=R(id(:),:);

