function Q = gramschmidt(A)
% Gram-Schmidt orthogonalization of column vectors
% Arguments: Matrix A passes vectors in its columns
% Return values: Matrix Q contains the orthornormal basis in its columns
[n,k] = size(A); % Get number k of vectors and dimension n of space
Q = A(:,1)/norm(A(:,1)); % First basis vector
for j=2:k
  q = A(:,j) - Q*(Q'*A(:,j)); % Orthogonal projection; loop-free implementation
  nq = norm(q);               % Check premature termination 
  if (nq < (1E-9)*norm(A(:,j))), break; end % Safe check for == 0 
  Q = [Q,q/nq]; % Add new basis vector as another column of Q
end
