function Q = gramschmidt(B)
% Gram-Schmidt orthonormalization of column vectors of matrix \texttt{B}.
% The columns of \texttt{B} are supposed to be \emph{linearly independent}!
m = size(B,2);
% Build output matrix with orthonormal columns by successively adding columns
% to the matrix \texttt{N} from the right.

% First normalized vector is just the normalized first column of \texttt{B}
Q = B(:,1)/norm(B(:,1));
for l=2:m
  % This is the \textbf{trick}: use the orthonormal system that has already been built for
  % \emph{efficient} implementation of orthogonal projection. First compute the
  % coefficient vector \texttt{c} of the orthogonal projection of the \Blue{$l$}-th
  % column of \texttt{B} onto the space spanned by the orthonormal vectors
  % constructed so far. This space agrees with the span of the first \Blue{$l-1$} columns
  % of \texttt{B}.
  c = (B(:,l)'*Q)',
  % The difference of the $l$-th column vector of \texttt{B} and its orthogonal
  % projection given by \texttt{N*c} determines the direction \texttt{n} of the \Blue{$l$}-th
  % orthonormal vector. Note that \texttt{Q*c} computes the orthogonal
  % projection of the \Blue{$l$}-th column of \texttt{B} onto the space spanned
  % by the first \Blue{$l-1$} columns of \texttt{Q}. 
  n = B(:,l)-Q*c; 
  % It remains to normalize it.
  Q = [Q, n/norm(n)];
end

