disp ('constructing sparse matrix')
n = 300; % size of the grid
m = 90000; % number of points
X = (n-1)*rand (m, 1); Y = (n-1)*rand (m, 1);
IX = ceil (X); JY = ceil (Y);

A = sparse(m, n^2);
A += sparse (1:m, sub2ind ([n, n], IX  , JY  ), (IX+1-X).*(JY+1-Y), m, n^2);
A += sparse (1:m, sub2ind ([n, n], IX+1, JY  ), (X - IX).*(JY+1-Y), m, n^2);
A += sparse (1:m, sub2ind ([n, n], IX  , JY+1), (IX+1-X).*(Y - JY), m, n^2);
A += sparse (1:m, sub2ind ([n, n], IX+1, JY+1), (X - IX).*(Y - JY), m, n^2);

printf ("A'*v matrix: ")
nv = 100;
v = ones (m, nv);
tic; 
u = A'*v;
printf ("%f s\n", toc)

printf ("A'*v 100 vectors: ")
v = ones (m, 1);
tic;
for i=1:nv
  u = A'*v;
endfor
printf ("%f s\n", toc)

printf ("(v'*A)' 100 vectors: ")
v = ones (m, 1);
tic;
for i=1:nv
  u = (v'*A)';
endfor
printf ("%f s\n", toc)

printf ("v*A' matrix: ")
v = ones (nv, m);
tic; 
u = v*A';
printf ("%f s\n", toc)

printf ("v*A' 100 vectors: ")
v = ones (1, m);
tic;
for i=1:nv
  u = v*A';
endfor
printf ("%f s\n", toc)