% limited memory Broyden method

function [dir, tau0, cache] = Direction_lbroyden(prob, opt, it, restart, sk, yk, v, cache)

sk = sk(:);
yk = yk(:);

[m, n] = size(v);
v = v(:);

if it == 1 || restart
    dir = -v; % use steepest descent direction initially
    cache.S = []; % stores vectors sk
    cache.Y = []; % stores vectors yk
    cache.W = []; % stores columns of H0*Y-S
    cache.StY = []; % stores inner products <sk,yk>
    cache.M = []; % these two are m-by-m where m is the memory of the method
    cache.LBroyden_mem = 0;
else
    % damping
    if opt.modBroyden == 2 % enforces positive curvature along sk
        sig = 0.1;
        prev_v = cache.prev_v;
        prev_tau = cache.prev_tau;
        sty = sk'*yk;
        stv = sk'*prev_v;
        if sty < sig*prev_tau*abs(stv)
            theta = (1+sign(stv)*sig)*prev_tau*stv/(prev_tau*stv + sty);
            yk = theta*yk - (1-theta)*prev_tau*prev_v;
        end
    end
    delta = 1; % diagonal of H0
    wk = delta*yk - sk;
    if cache.LBroyden_mem == opt.memory, idx0 = 2;
    else idx0 = 1; cache.LBroyden_mem = cache.LBroyden_mem+1; end
    S0 = cache.S(:,idx0:end);
    Y0 = cache.Y(:,idx0:end);
    W0 = cache.W(:,idx0:end);
    StY0 = cache.StY(idx0:end,idx0:end);
    M0 = cache.M(idx0:end,idx0:end);
    % update matrices S, Y, W, StY, M
    cache.S = [S0, sk];
    cache.Y = [Y0, yk];
    cache.W = [W0, wk];
    if isempty(Y0), cache.StY = sk'*yk;
    else cache.StY = [[StY0; sk'*Y0], cache.S'*yk]; end
    if isempty(S0), cache.M = 0;
    else cache.M = [[M0; S0(:,end)'*S0], zeros(cache.LBroyden_mem, 1)]; end
    K = delta * cache.StY - cache.M;
    % compute direction
    dir = delta*(cache.W * (K\(cache.S'*v)) - v);
end

tau0 = 1.0;
dir = reshape(dir, m, n);

end