conjgrad(x0, b, operator, adjoint, res_limit, max_steps)

function [x, Res] = conjgrad(x0, b, operator, adjoint, res_limit, max_steps)
% Uses the method of conjugate gradients to minimze the following
% expression:
% f(x) = (A*x-b)'*(A*x-b)
%
% See http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
%  and http://en.wikipedia.org/wiki/Conjugate_gradient_method
%
% Inputs:
%  x0: initial value, must be of correct size
%  b:  as in expression above
%  operator: handle to a function taking one parameter x. It computes A*x.
%  adjoint:  handle to a function taking one parameter y. It computes A'*y.
%  res_limit: relative residual limit sqrt(f(x)/(b'*b)) stopping criterion
%  max_steps: maximum number of steps stopping criterion
%
% Output:
%  x:  the optimized input vector
%  Res: vector of relative residuals: sqrt(f(x)/(b'*b))
%
% (c) 2009 Peter Divos - pdivos at gmai dot com

x = x0;

Res = [];

dX = operator(x)-b;

res = sqrt(mydot(dX,dX)/mydot(b,b));
Res = [Res, res];

Dx = adjoint(dX);
Lx = Dx;
LX = operator(Lx);
alpha = -mydot(dX,LX)/mydot(LX,LX);
x = x + alpha * Lx;

while(1)
    Lx_prev = Lx;
    Dx_prev = Dx;

    dX = operator(x)-b;
    res = sqrt(mydot(dX,dX)/mydot(b,b));
    Res = [Res, res];
    if(res<res_limit || numel(Res)>max_steps)
        break;
    end
    Dx = adjoint(dX);

    beta = mydot(Dx,Dx)/mydot(Dx_prev,Dx_prev);
    Lx = Dx + beta * Lx_prev;

    LX = operator(Lx);
    alpha = -mydot(dX,LX)/mydot(LX,LX);
    x = x + alpha * Lx;
end

    function d = mydot(x1,x2)
        % my dot product for any number of dimensions
        if(size(x1)~=size(x2))
            error('x1 and x2 not of same size');
        end
        d = x1.*conj(x2);
        for ii = 1:numel(size(x1))
            d = sum(d);
        end
    end
end