function [ p ] = LBFGS_two_loop( m, grad, s, y ,seq)
% two loop L-BFGS 
% Input :

% m    - number of former vectors used
% grad - gradient of the current point
% s    - n by M matrix, record the difference of points
% y    - n by M matrix, record the difference of gradients

% Output:

% p    - descent direction

 
q = grad;
alpha_rec = zeros(m,1);
[~,seqsort] = sort(seq,'ascend');

for i = m : -1 : 1
    alpha = s(:,seqsort(i))'*q/(y(:,seqsort(i))'*s(:,seqsort(i)));
    q = q - alpha*y(:,seqsort(i));
    alpha_rec(seqsort(i)) = alpha;
end
H0 = y(:,seqsort(m))'*s(:,seqsort(m))/(y(:,seqsort(m))'*y(:,seqsort(m)));
z = H0*q;

for i = 1 : m
    beta = y(:,seqsort(i))'*z/(y(:,seqsort(i))'*s(:,seqsort(i)));
    z = z + s(:,seqsort(i))*(alpha_rec(seqsort(i))-beta);
end
 
p = -z;

end

