function [B1] = Bfull_convex(S,Y,gamma_inv,phi_convex)
%
% Computes the full quasi-Newton matrix for debugging purposes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% "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.
%---------------------------------------------------------------------------    
% 
% Inputs: S,Y: the quasi-Newton pairs that define the matrix
%         gamma_inv: the initial quasi-Newton matrix (B_0 = gamma_inv*I)
%         R1: the leading kxk upper triangular matrix of the QR decomposition
%             of Psi
%         phi_convex is the parameter for the convex class:
%              Phi_convex = 0 (BFGS), Phi_convex=1 (DPF)
%         flag = 0 if phi is of the form [B0SK YK]
%         flag = 1 if phi has the last columns are new [phi B0snew ynew]
%%
% Output: B1 is the full matrix 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n,k]=size(S);

B1 = gamma_inv*eye(n);

for i=1:k
  Bs = B1*S(:,i);
  w  = (1/(Y(:,i)'*S(:,i))).*Y(:,i) - (1/(S(:,i)'*Bs))*Bs;
  B1 = B1 - (1/(S(:,i)'*Bs))*Bs*(Bs)'+(1/(Y(:,i)'*S(:,i)))*Y(:,i)*Y(:,i)'+phi_convex*(S(:,i)'*Bs)*w*w';
end

