function [dEE] = eig_finder(S,Y,R1,gamma_inv)
%
% Computes the eigenvalues of the SR1 matrix given S, Y, R1, gamma_qn
%
% "On efficiently computing the eigenvalues of limited-memory
% quasi-Newton matrices"
% by Jennifer Erway and Roummel Marcia
%
% Copyright (2015): Jennifer Erway and Roummel Marcia
%
% The technical report and software are available at 
% www.wfu.edu/~erwayjb/publications.html
% www.wfu.edu/~erwayjb/software.html
%
% 
% This code is distributed under the terms of the GNU General Public
% License
% 2.0.
%
% Permission to use, copy, modify, and distribute this software for
% any purpose without fee is hereby granted, provided that this entire 
% notice is included in all copies of any software which is or includes
% a copy or modification of this software and in all copies of the
% supporting documentation for such software.
% This software is being provided "as is", without any express or
% implied warranty.  In particular, the authors do not make any
% representation or warranty of any kind concerning the merchantability
% of this software or its fitness for any particular purpose.
%---------------------------------------------------------------------------    
% 
% Inputs: S,Y: the quasi-Newton pairs that define the SR1 matrix
%         gamma_inv: the initial quasi-Newton matrix (B_0 = gamma_inv*I)
%         R1: the leading kxk upper triangular matrix of the QR decomposition
%             of Psi
%
% Output: dEE: the eigenvalues of the SR1 matrix (in a vector)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%constants
[m,k] = size(S);
dEE = zeros(m,1);

%decomposition of M
STY   = S'*Y;
M     = inv(tril(STY)+(tril(STY,-1))'-gamma_inv*S'*S); 
M     = (M+M')./2;  

%eigs computation for smaller kxk block
R1MR1  = R1*M*R1';
R1MR1  = (R1MR1+(R1MR1'))./2;
[v,d1] = eig(R1MR1);
d1     = real(diag(d1));  %store as a vector, avoid complex problems

dEE(1:k)   = d1+gamma_inv;
dEE(k+1:m) = gamma_inv;

