function [xk, f, xl] = trunc_newton (obj, x0, Hessian, converg, eps, epsilon)
%TRUNC_NEWTON   Truncated Newton (ACM 113, Problem 3.4)
%
%   [xmin, fmin, xl] = trunc_newton(obj, x0, Hessian, converg, eps, epsilon)
%
%       obj:         Objective function
%                    (obj(x, n) is the n-th deriv of obj at x)
%       x0:          Initial guess of x_min
%       HessianFree: 0 - not Hessian free,
%                    1 - finite difference method
%                    2 - complex variable method
%       eps:         epsilon for Hessian free method
%
%       converg:     Expected convergence rate, usually 1--2.
%       epsilon:     epsilon for outer Newton iteration.
%
%   By default, HessianFree = 2, eps = 1e-8;
%               converg = 1.5; epsilon = 1e-12.
%
%   See also: LINESCH, TRUST_REGION.

%   Ling Li, Caltech
%   05/05/2001

if nargin < 3,
    Hessian = 2;
end
if nargin < 4,
    converg = 1.5;
end
if nargin < 5,
    %% since e_u = 1e-16, i.e. (1+e_u)-1 == 0 in the machine
    %% we set eps = sqrt(e_u).
    eps = 1e-8;
end
if nargin < 6,
    epsilon = 1e-12;
end


% Outer Newton iteration
%% Initialize
out_iter = 0;
xk = x0; xl = xk;
f = eval([obj '(xk, 0)']); df = eval([obj '(xk, 1)']);

while norm(df) > epsilon * (1+abs(f)),
    if not(Hessian),
        A = eval([obj '(xk, 2)']);
    end

    %% Inner linear CG: solve d2f p = -df
    in_iter = 0;
    xi = 0;
    ri = df; %% A*xi - b;
    pi = - ri;

    ndf = norm(df);
    switch converg
    case 1,
        eta = 0.5;
    case 2,
        eta = min(0.5, ndf);
    otherwise,
        eta = min(0.5, ndf^(converg-1));
    end

    while norm(ri) > eta * ndf,
        switch Hessian
        case 0,
            Api = A * pi;
        case 1, %% d2f * pi -> (df(x+ep) - df(x)) / e;
            dfep = eval([obj '(xk + eps.*pi, 1)']);
            Api = (dfep - df) ./ eps;
        case 2, %% d2f * pi -> Im(df(x+i*ep)) / e;
            Api = imag(eval([obj '(xk + 1i*eps.*pi, 1)'])) ./ eps;
        otherwise,
            error('Wrong HessianFree method!');
        end

        piApi = pi' * Api;
        if piApi <= 0,
            if in_iter == 0,
                xi = pi;
            end
            break;
        end

        ri2 = ri' * ri;
        ai = ri2 / piApi;
        xi = xi + ai * pi;
        ri = ri + ai .* Api;
        bi = ri' * ri / ri2;
        pi = bi * pi - ri;
        in_iter = in_iter + 1;
    end
    %% Inner linear CG: END

    pk = xi;
%    disp(['real df = ' num2str((-inv(eval([obj '(xk,2)']))*df)')]);
%    disp(['  now p = ' num2str(pk')]);

    %ak = linesch(xk, pk, 1, 0.5, f, df);
    ak = 1;
%    disp(['  alpha_' int2str(out_iter) ' = ' num2str(ak)]);
%    disp(' ');

    out_iter = out_iter + 1;
    xk = xk + ak * pk;
    xl = [xl, xk];
    f = eval([obj '(xk, 0)']);
    df = eval([obj '(xk, 1)']);
end

disp('Result: ');
disp(['  x_min = [' num2str(xk') ']'' with norm = ' num2str(norm(xk)) '.']);
disp(['  f_min = ' num2str(f)]);
disp(['  iterations = ' int2str(out_iter)]);
