  function [mode,x,y]=aprod(mode,m,n,x,y)
    %  user supplied function for nag routine fo4qaf
global Kop
x=sparse(x);
y=sparse(y);
if mode == 1,
  y=y+Kop*x;
  y=full(y);
  x=full(x);
elseif mode == 2
  x=x+Kop'*y;
  x=full(x);
  y=full(y);
else
  'error in APROD - wrong mode'
end
