function A = makematrix(ndim);
%-------------------------------------------------------------------
%
%   Call:  A = makematrix(ndim);
%
%   Given parameter "ndim",  this function returns
%   a square symmetric matrix A of the size N x N,
%   where N = ndim^2.
%   Matrix A represents a blurring operator, in the following sense.
%   If "f" is a square image of size ndim x ndim, and "X"
%   is reordering of "f" into a column vector of size ndim^2 x 1
%   then  Y = A X represents a blurred image (after reshaping
%   "Y" into an image of size ndim x ndim)
%
%-------------------------------------------------------------------


ndimx = ndim;
ndimy = ndim;

iconvsize = 3;

lengvec = ndimy*ndimx;

A =zeros(lengvec,lengvec,'single');


for iy = 1:ndimy
    for ix = 1:ndimx
        iveccoor = (ix-1)*ndimy+iy;

        for newy = (max(iy-iconvsize,1)):(min(iy+iconvsize,ndimy))
            for newx = max(ix-iconvsize,1):(min(ix+iconvsize,ndimx))
                ioldcoor = (newx-1)*ndimy+newy;
                dx = ix - newx;
                dy = iy - newy;
                weight = 1./(1.+dx*dx+dy*dy);

                A(iveccoor,ioldcoor) = weight/sqrt(ndim);
            end
        end
    end
end


A = A*A*A;


return

