function [R1] = update_QR_add_pair(S,Y,gamma_inv,R1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Updates the R part of the QR decomposition of Psi
% after a new pair has been computed and added to S, Y, and Psi.
%
% "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)
%         Psi is a mx(k+1) matrix
%         R1, the leading kxk upper triangular matrix of the QR decomposition
%             of Psi before the newest column was added to Psi
%
% Output: R1, the (k+1)x(k+1) upper triangular matrix of the QR
%         decomposition of Psi after the addition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
feedback = 1;  %flag

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


if rank(R1)==size(R1,2)  %compute numerical rank to MATLAB tolerances


  %add kplus-1 column to R1
  u_1 = R1'\(Psihat(:,1:kplus-2)'*Psihat(:,kplus-1));
  eta = sqrt((norm(Psihat(:,kplus-1)))^2 - (norm(u_1))^2);  
  R1  = [R1 u_1; zeros(1,kplus-2) eta];


  %add kplus column to R1
  u_1 = R1'\(Psihat(:,1:kplus-1)'*Psihat(:,kplus));
  eta = sqrt((norm(Psihat(:,kplus)))^2 - (norm(u_1))^2);  
  R1  = [R1 u_1; zeros(1,kplus-1) eta];


else
  if feedback
    fprintf('\nR is not full rank... redoing full QR decomposition\n');
  end
  R   = qr(Psihat);
  R1  = R(1:kplus,1:kplus);
  fprintf('error check rTr-PsiTPsi: %8.2e', norm(Psihat'*Psihat-R1'*R1));
end

