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));