function [x,iters,cgft] = ...
    PCGsolver(b,TF,Hfit,Hreg,eta,Amult,tol,maxiters,precond)
%PCGSOLVER   Preconditioned Conjugate-Gradient method.
%   PCGSOLVER solves the linear system 
%             (A'*diag(HFIT(:))*A+ETA*HREG) X = B
%   whose coefficient matrix is the hessian of the objective function
%   J=J0+ETA*JR where J0 is the data-fitting term and JR is the 
%   regularization term. The matrix A'*diag(HFIT(:))*A is the hessian of J0  
%   and HREG is the hessian of JR. 
%   For deblurring applications, A is the blurring operator while, 
%   for denoising applications, A is the identity.  
%   The initial iterate is set to zero.
%
%   [X]=PCGSOLVER(B,TF,HFIT,HREG,ETA,AMULT,TOL,MAXITERS,PRECOND)
%   returns the (NX*NY)-by-1 solution X. The NX-by-NY array B is the 
%   right hand side. The NX-by-NY array TF is the fft of the psf image. 
%   HFIT is a NX-by-NY array such that A'*diag(HFIT(:))*A is the 
%   hessian of the data-fitting term J0. 
%   The (NX*NY)-by-(NX*NY) array HREG is the (possibly approximate) hessian
%   of the regularization term JR. ETA is the regularization parameter.
%   AMULT is a function handle. AMULT(P,HFIT,TF,NX,NY) accepts the input 
%   arrays P, HFIT and TF and returns the matrix-vector product 
%   (A'*diag(HFIT(:))*A)*P.
%   TOL is the relative residual tolerance. MAXITERS is the maximum
%   number of iterations. PRECOND specifies if the diagonal preconditioner
%   is used. If PRECOND is 1, PCGSOLVER uses the diagonal preconditioner
%   else, if precond is 0, a preconditioner is not applied. 
%
%   [X,ITERS]=PCGSOLVER(B,TF,HFIT,HREG,ETA,AMULT,TOL,MAXITERS,PRECOND)
%   also returns the number ITERS of performed iterations.
%
%   [X,ITERS,CGFT]=PCGSOLVER(B,TF,HFIT,HREG,ETA,AMULT,TOL,MAXITERS,PRECOND)
%   also returns the number CGFT of performed FFTs.
%
%   See also AMULT.

[nx,ny] = size(b); N = nx*ny; 
b = reshape(b,N,1);
tol = tol*norm(b); % residual relative tolerance 
cgft = 0;

% set initial iterate
x = zeros(size(b));
r = b;
if precond 
    if isempty(TF)
        dd1 = Hfit(:);
        dd2 = diag(Hreg);
        denom = dd1+eta*dd2;
        z = r./denom;
    else
        dd1 = mean(Hfit(:));
        dd2 = mean(diag(Hreg));
        denom = dd1*abs(TF).^2+eta*dd2;
        z = (real(ifft2(fft2(reshape(r,nx,ny))./denom))) ; z = z(:);
        cgft = cgft+2;
    end    
else
    z = r;
end
rho = z'*r; p = z;
iters = 1;
eflag = 1;

while eflag
    % computes the matrix vector product of p with A'*diag(Hfit(:))*A
    [temp,nft] = feval(Amult,p,Hfit,TF,nx,ny); cgft = cgft+nft;  
    Ap = temp+eta*Hreg*p;
    w = Ap;
    alpha = p'*w;
    alpha = rho/alpha;
    % update CG iterate and residual
    x = x+alpha*p;
    r = r-alpha*w;
    rhoold = rho;
    if precond 
        if isempty(TF)
            z = r./denom;
        else
            z = (real(ifft2(fft2(reshape(r,nx,ny))./denom))) ; z = z(:);
            cgft = cgft+2;
        end  
    else
        z = r;
    end
    rho = z'*r;
    beta = rho/rhoold;
    p = z+beta*p;
    iters = iters+1; 
    % check stopping criteria
    eflag = (norm(r) > tol) & (iters <= maxiters);   
end
iters = iters-1;