function [ Q, R ] = gs_factor ( A )
% [ Q, R ] = gs_factor ( A )
% Use modified Gram-Schmidt to construct Q and R factors of A, where
% Q is an orthogonal matrix
% R is an upper-triangular matrix

% $Id: gs_factor.m,v 1.2 2017/02/08 15:19:40 mike Exp $
% M. M. Sussman

for k=1:size(A,2)
  y=A(:,k);
  for ell=1:k-1
    R(ell,k)=dot(Q(:,ell),y);
    y=y-R(ell,k)*Q(:,ell);
  end
  R(k,k)=norm(y);
  if R(k,k) > 0
    Q(:,k)=y/R(k,k);
  else
    Q(:,k)=0;
  end
end
