function [dEE,R1] = eig_finder_convex(S,Y,gamma_inv,R1,phi_convex,flag)
%
% Computes the eigenvalues of the quasi-Newton convex class matrix given S, Y, Psi, 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 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
%         phi_convex is the parameter for the convex class:
%              Phi_convex = 0 (BFGS), Phi_convex=1 (DPF)
%         flag = 0 if phi is of the form [B0SK YK]
%         flag = 1 if phi has the last columns are new [phi B0snew ynew]
%%
% Output: dEE: the eigenvalues of the quasi-Newton matrix (in a vector)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

%check Psi has full column rank 
Psi   = [gamma_inv.*S Y];
if rank(Psi)<size(Psi,2)
  ferror('Psi is not full rank... rerun');
end


sty = S(:,1)'*Y(:,1);
sBs = gamma_inv.*S(:,1)'*S(:,1);
M11 = -(1-phi_convex)/(sBs);
M12 = -phi_convex/(sty);
M22 = (1/sty)*(1 + ((phi_convex*sBs)/(sty)));
M = [M11 M12; M12 M22]; 
for j=1:k-1
  sty = S(:,j+1)'*Y(:,j+1);
  PsiTJminusS = [Psi(:,1:j) Psi(:,k+1:k+j)]'*S(:,j+1);
  pj = M*(PsiTJminusS);
  sBs = gamma_inv.*(S(:,j+1)'*S(:,j+1)) + (PsiTJminusS'*pj);
  alphaj = -(1-phi_convex)/(sBs);
  betaj = -phi_convex/(sty);
  deltaj = (1+(phi_convex*sBs)/(sty))/(sty);
  M =  [ M+alphaj*pj*pj'   alphaj*pj  betaj*pj; alphaj*pj'   alphaj ...
	 betaj;  betaj*pj'  betaj deltaj];
  % apply permutations
  % move 2j-2 column to (j)th col, shift rest
  M = M(:, [1:j 2*j+1 (j+1):2*j 2*j+2]);
  % move 2j-2 row to jth row and shift rest
  M = M([1:j 2*j+1 (j+1):(2*j) 2*j+2],:);
end


%resymmetrize
M     = (M+M')./2;  

%create the perfect shuffle permutation
permute_ind =1:2*k;
permute_ind=reshape(permute_ind,length(permute_ind)/2,2)';
permute_ind=permute_ind(:);


if flag == 0

  %permute Psi and M using perfect shuffle to get Psihat, PMP
  Psihat=Psi(:,permute_ind);
  PMP  = M(:,permute_ind);
  PMP =  PMP(permute_ind,:);
  [Q,R] = qr(Psihat);
  R1    = R(1:2*k,1:2*k);


elseif flag == 1  %have r
  %permute Psi and M using perfect shuffle to get Psihat, PMP
  Psihat=Psi(:,permute_ind);
  PMP  = M(:,permute_ind);
  PMP =  PMP(permute_ind,:);  

else 
  ferror('\nflag must either be set to 0 or 1\n\n');
end
  

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

