%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This is a demo driver for gen_recursionSolve_qN.m.  
% 
% For details, see
%
% "Solving limited-memory BFGS systems with 
%  generalized diagonal upates"
% by Jennifer Erway and Roummel Marcia
%
% which is available at
%
% http://www.wfu.edu/~erwayjb/publications.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

%set parameters
M     =    5;   %number of L-BFGS updates
n     = 1000;   %dimension of problem
epsil = 1e-6;   %lower bound on gamma*sigma

%generate a positive-definite diagonal matrix
D = rand(n,1);

%form the matrices S and Y from quasi-Newton method
S = randn(n,M);
Y = randn(n,M);

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

%ensure min(D)>=sigma
%set 1/gamma and ensure that gamma*sigma > epsil
minD = min(D);
gamma_inv = norm(Y(:,M))^2/(S(:,M)'*Y(:,M));
gamma_inv = min(gamma_inv, minD/epsil);

%right-hand side
z = randn(n,1);

%unroll the L-BFGS formula
[a,b] = unrolling(S,Y,gamma_inv);

tic;
%solve (B + D)x=z
[x] = gen_recursionSolve_qN(a,b,gamma_inv,D,z);
time_recur = toc;

%solve using direct method
B = gamma_inv*eye(n);
for i=1:M
    B = B - a(:,i)*a(:,i)'+ b(:,i)*b(:,i)';
end
tic;
t = (B+diag(D))\z;
time_direct = toc;

%test recursion solution quality by computing the
%residual relative error ||(B + D)*x - z||/||z||
resid = norm( B*x+D.*x - z)/norm(z);
fprintf('\n');
fprintf('Residual relative error: %8.2e\n\n',     resid);
fprintf('Recursion solve time:    %8.6f sec\n',   time_recur);
fprintf('Direct solve time:       %8.6f sec\n\n', time_direct);

