function [temp,nft] = Amult(p,Hfit,TF,nx,ny)
%AMULT   Matrix-vector product.
%
%   [TEMP]=AMULT(P,HFIT,TF,NX,NY) evaluates the matrix-vector
%   product TEMP of the (NX*NY)-by-1 vector P with the matrix
%   A'*diag(HFIT(:)*A:
%      TEMP = (A'*diag(HFIT(:)*A)*P
%   where A is the blurring operator. In denoising applications, 
%   A is the identity. 
%   For Poisson-corrupted data, HFIT is a NX-by-NY array such that 
%   A'*diag(HFIT(:))*A is the hessian of the Kullback-Leibler divergence. 
%   For Gaussian-corrupted data, HFIT=1. 
%   For deblurring applications, the NX-by-NY array TF is the fft2 
%   of the psf image. For denoising applications, TF=[]. 
%   Note that TEMP is a (NX*NY)-by-1 vector.
%
%   [TEMP,NFT] = AMULT(P,HFIT,TF,NX,NY) also returns the number NFT of
%   performed ffts.
%
%   See also KLTKB_OBJ, KLTVB_OBJ, KLTVN_OBJ, LSTKB_OBJ, LSTVB_OBJ, 
%   LSTVN_OBJ, PCGSOLVER.

if isempty(TF) % denoising 
    temp = Hfit(:).*p;
    nft = 0;
else % deblurring
    temp = real(ifft2(TF.*fft2(reshape(p,nx,ny))));
    temp = Hfit.*temp;
    temp = real(ifft2(conj(TF).*fft2(temp)));
    temp = temp(:);
    nft = 4;
end

