function [xk,lk,f] = aug_lagran (obj,x0,l0,a0,Hessian,converg,eps,epsilon)
%AUG_LAGRAN     Augmented Lagrangian (ACM 113, Final Problem 4)
%
%   [xmin, lmin, fmin] =
%           aug_lagran(obj, x0, lambda0, a, 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 (xmin)
%       l0:          Initial guess of lambda_min (lmin)
%       a:           value of a, or [a0, inc_rate];
%       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: TRUNC_NEWTON.

%   Ling Li, Caltech
%   06/06/2001

global l a

% Default parameters
if nargin < 5, Hessian = 2;     end
if nargin < 6, converg = 1.5;   end
if nargin < 7, eps = 1e-8;      end
if nargin < 8, epsilon = 1e-12; end


xk = x0; nx = length(x0);
lk = l0; nl = length(l0);
ak = a0(1);
ite = 0;

% store xk and lk to calc converge rate
xl = []; ll = [];
d = Inf;

while d > epsilon,
    xl = [xl, xk]; ll = [ll, lk];
    l = lk; a = ak;
    [xn, f] = trunc_newton(obj, xk, Hessian, converg, eps, epsilon);
    c = eval([obj '(xn, -1)']);

    d = norm(xn - xk) + ak*norm(c);
    xk = xn;
    lk = lk - ak * c;

    ite = ite + 1;
    if length(a0) > 1,
        ak = ak * a0(2);
    end
end

% Calculate distance to the final point.
xl = xl - repmat(xk, 1, ite); xd = sqrt(sum(xl.*xl));
ll = ll - repmat(lk, 1, ite); ld = sqrt(sum(ll.*ll));

semilogy(0:ite-1, [xd; ld], 'o-');
legend('error of x', 'error of \lambda');
xlabel('iteration #');
ylabel('error');
