function [DEE] = driver_eig_finder_convex()
%
% Driver for eig_finder_convex.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% "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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

phi_convex = 0.5;  %parameter for Broyden convex class
                 %phi_convex=0 (BFGS), phi_convex=1 (DFP)

%implicitly form quasiNewton matrix with random data
m         = 100;
k         = 5;
S         = randn(m,k); Y=randn(m,k);
gamma_inv = 3;

%ensure B is positive definite
for i = 1:k
  if S(:,i)'*Y(:,i)<0
    S(:,i) = -S(:,i);
  end
end



%find eigenvalues of B
flag = 0;
R1 = 0;
[dEE,R1] = eig_finder_convex(S,Y,gamma_inv,R1,phi_convex,flag);


%check
B       = Bfull_convex(S,Y,gamma_inv,phi_convex);
B       = (B+B')./2;
sortedD = sort(dEE);
sortedB = sort(real(eig(B)));  %knock out rounding errors that are complex
fprintf('\nphi: %2.2e', phi_convex);
fprintf('\nn: %8d', size(S,1));
fprintf('\nTesting eigenvalues: %41.5e\n', norm(sortedD-sortedB,Inf)/norm(sortedB,Inf));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%add a new pair to existing S and Y

s_new   = randn(m,1);
y_new   = randn(m,1);

if s_new'*y_new<0
  s_new = -s_new;
end
k       = k+1;
S       = [S s_new];
Y       = [Y y_new];


[R1]    = update_QR_add_pair_convex(S,Y,gamma_inv,R1);
flag = 1;  %R1 is known
[dEE]   = eig_finder_convex(S,Y,gamma_inv,R1,phi_convex,flag);

%test accuracy of method
B       = Bfull_convex(S,Y,gamma_inv,phi_convex);
sortedD = sort(dEE);
B       = (B+B')./2;
sortedB = sort(real(eig(B))); %knock out rounding errors that are complex
fprintf('Testing eigenvalues after adding a pair: %21.5e\n', norm(sortedD-sortedB,Inf)/norm(sortedB,Inf));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%delete first pair and then add a new pair to S and Y
S       = S(:,2:k);
Y       = Y(:,2:k);

[R1]    = update_QR_delete_pair_convex(R1);

%%% form psihat
Psi = [gamma_inv.*S Y];
[m,kplus] = size(Psi);
permute_ind =1:kplus;
permute_ind=reshape(permute_ind,length(permute_ind)/2,2)';
permute_ind=permute_ind(:);
Psihat = Psi(:,permute_ind);

%add new pair
s_new   = randn(m,1);
y_new   = randn(m,1);
S       = [S s_new];
Y       = [Y y_new];

[R1]    = update_QR_add_pair_convex(S,Y,gamma_inv,R1);

flag = 1;  %already have R1
[dEE]   = eig_finder_convex(S,Y,gamma_inv,R1,phi_convex,flag);

%check accuracy of method
B       = Bfull_convex(S,Y,gamma_inv,phi_convex);
sortedD = sort(dEE);
B       = (B+B')./2;
sortedB = sort(real(eig(B)));  %knock out rounding errors that are complex
fprintf('Testing eigenvalues after deleting and adding a pair: %8.5e\n\n', norm(sortedD-sortedB,Inf)/norm(sortedB,Inf));
