function lbfgs(n)
% L-BFGS, Algorithm 7.5 in N&W
% n = problem dimension

m = 5;
x0 = zeros(n,1);
xE = ones(n,1); % Exact
Hk0 = eye(n);
[fk, g0] = rosenbrock(x0);

% linesearch parameters
alpha_ini = 1.0;
alpha_min = 1.0E-08;
alpha_max = 1.0E+02;
c1        = 1.0E-04;
c2        = 0.9E-01;
xtol      = 1.0E-08;
nfuneval  = 20;
quiet     = 1;

% Vector of vectors to store sk and yk
s = zeros(m,n);
y = zeros(m,n);

Nmax  = 1000;
xk = x0;
gk = g0;

for k = 1:Nmax
      
    if norm(gk) <= 1.0E-12*norm(g0)
        % converged!
        break;
    end
    
    pk = twolooprec(Hk0, gk, m, s, y);
    % strong Wolfe linesearch
    [~,~,~,~,alpha,info,~] = minpack_cvsrch(@rosenbrock,xk,       ...
                                            fk,gk,pk,alpha_ini,    ...
                                            c1,c2,xtol,          ...
                                            alpha_min,alpha_max, ...
                                            nfuneval,quiet);
    if info ~= 1
        fprintf('Warning: linesearch routine returned error %d', ...
                     info);
    end
    xk = xk + alpha*pk;
    
    [fk, gk1] = rosenbrock(xk);
    
    fprintf('%3d: f=%10.3e, |df|=%10.3e, |x-x*| = %10.3e\n', ...
            k, fk, norm(gk1), norm(xk-xE));
    
    % Update storage
    s = [s(2:m,:); zeros(1,n)]; % Remove last and move others upwards
    y = [y(2:m,:); zeros(1,n)];
    s(m,:) = alpha*pk'; % Store new values
    y(m,:) = gk1 - gk;
    
    gamma = (s(m,:)*y(m,:)') / (y(m,:)*y(m,:)');
    Hk0 = gamma*eye(n);
    gk = gk1;
end

end


function p = twolooprec(Hk0, gk, m, s, y)
% Algorithm 7.4
q = gk;
alpha = zeros(m,1);
for i=m:-1:1
    rho_inv = s(i,:)*y(i,:)';
    if (rho_inv > 0)
        alpha(i) = (s(i,:)*q) / rho_inv;
        q = q - alpha(i)*y(i,:)';
    end
end
r = Hk0*q;
for i=1:m
    rho_inv = s(i,:)*y(i,:)';
    if (rho_inv > 0)
        beta = (y(i,:)*r) / rho_inv;
        r = r + s(i,:)'*(alpha(i) - beta);
    end
end
p = -r;

end


% Backtracking algorithm
function alpha=armijo(x,p)
alpha = 1.0; rho = 0.5; c1 = 1.0E-03;
[f0, df0] = rosenbrock(x);
df0 = c1*df0'*p;
while rosenbrock(x+alpha*p) > f0+alpha*df0
    alpha = alpha*rho;
end
end