function [out] = recursionSolve_qN(a,b,b0,sigma,z)
%
%
% This function directly solves systems of the form (B+\sigma I)x=z, 
% where B is an L-BFGS quasi-Newton matrix.  
%
% "Limited-Memory BFGS Systems with Diagonal Updates" 
% by Jennifer Erway and Roummel Marcia
%
% Copyright (2011): 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: a(:,k) and b(:,k) make up the kth L-bfgs update and
% B_0 = b0*I is the initial L-BFGS matrix, i.e., 
% B_k = B_0 - sum_{i=1}^k a(:,i)a(:,i)^T + sum_{i=1}^{k} b(:,i)b(:,i)^T

% initialization
k0 = size(a,2);
k  = 2*k0;
inv_c0 = 1/( b0 + sigma );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Recursion formula
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
invCz  = inv_c0*z;

for j = 1:k;
    if (mod(j,2)==1)
        c  = a(:,(j+1)/2);
    else
        c  = b(:,j/2);
    end
    p(:,j) = inv_c0*c;
    for t = 1:j-1
        p(:,j) = p(:,j) + (-1)^(t+1)*v(t)*(p(:,t)'*c)*p(:,t);
    end
    v(j) = 1/(1 + ((-1)^j)*c'*p(:,j));
    invCz = invCz + (-1)^(j+1)*v(j)*(p(:,j)'*z)*p(:,j);

end

out = invCz;



