function C=CovarianceBjoerck(R)
% COVARIANCEBJOERCK computes the covariance matrix by Bjoerck's method
%   C=CovarianceBjoerck(R) computes the covariance matrix from the
%   R factor of A=QR by solving the system R C=R^-T

n=min(size(R)); 
C=zeros(n,n);
for k=n:-1:1,                  
  C(k,k)=(1/R(k,k)-R(k,k+1:n)*C(k,k+1:n)')/R(k,k);
  for i=k-1:-1:1,
    C(i,k)=-(R(i,i+1:k)*C(i+1:k,k)+R(i,k+1:n)*C(k,k+1:n)')/R(i,i);
  end
end