function [r] = two_loop(vec,gamma_qn,stored_vecs)
%multiplies "vec" with the inverse of the L-BFGS matrix
%assumes initial L-BFGS inverse is gamma_qn*I, i.e., B_0=(1/gamma_qn)I
%an implementation based on the two-loop recursion published in
%NOCEDAL, J. 1980. Updating quasi-Newton matrices with limited storage.
%Math. Comput. 35, 773--782.
global S Y RHO;
q = vec;
%allocating memory
alfa = zeros(stored_vecs,1);
for i = stored_vecs:-1:1
alfa(i) = RHO(i)*S(:,i)'*q;
q = q - alfa(i)*Y(:,i);
end
r = gamma_qn*q;
for i = 1:stored_vecs
beta = RHO(i)*Y(:,i)'*r;
r = r + (alfa(i)-beta)*S(:,i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%