function [DEE] = driver_eig_finder()
%
% Driver for eig_finder.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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

%set up for eig_finder
Psi   = Y-gamma_inv.*S;
if rank(Psi)<size(Psi,2)
  ferror('Psi is not full rank... rerun');
end

%full QR
[Q,R] = qr(Psi);
R1    = R(1:k,1:k);

%find eigenvalues of B
[dEE] = eig_finder(S,Y,R1,gamma_inv);

%check accuracy of method
B       = Bfull(S,Y,gamma_inv);
B       = (B+B')./2;
sortedD = sort(dEE);
sortedB = sort(real(eig(B)));  %knock out rounding errors that are complex
fprintf('\nUsing an SR1 matrix');
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);
k       = k+1;
S       = [S s_new];
Y       = [Y y_new];
Psi     = [Psi y_new-gamma_inv*s_new];

if rank(Psi)<size(Psi,2)
  ferror('Psi is not full rank... rerun');
end


[R1]    = update_QR_add_pair(Psi,R1);
[dEE]   = eig_finder(S,Y,R1,gamma_inv);

%check
B       = Bfull(S,Y,gamma_inv);
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);
Psi     = Psi(:,2:k);
[R1]    = update_QR_delete_pair(Psi,R1);

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

if rank(Psi)<size(Psi,2)
  ferror('Psi is not full rank... rerun');
end

[R1]    = update_QR_add_pair(Psi,R1);
[dEE]   = eig_finder(S,Y,R1,gamma_inv);

%check accuracy of method
B       = Bfull(S,Y,gamma_inv);
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));
