function [Q,R] = mgs(A)
% FUNCTION [Q,R] - mgs(A)
% Modified Gram Schmidt Algorithm
% find Q,R such that A = QR
%
% Input:    A -- matrix to be factorized
% Output:   Q -- unitary matrix
%           R -- upper triangular matrix

  [m,n] = size(A);
  Q = zeros(m,n); % initialization of Q
  R = zeros(n,n); % initialization of R

  for i = 1:n
    R(i,i) = norm(A(:,i));
    Q(:,i) = A(:,i)/R(i,i);
    for j = i+1:n
      R(i,j) = Q(:,i)'*A(:,j);
      A(:,j) = A(:,j) - R(i,j)*Q(:,i);
    end
  end