function [Q,R] = mgs(A);

    p = size(A,2);
    Q = A;
    R = zeros(p,p);
    
    for m=1:p,
                
        for j=1:m-1,
            
            R(j,m) = Q(:,j)'*Q(:,m);
            Q(:,m) = Q(:,m) - R(j,m)*Q(:,j);
            
        end
        
        R(m,m) = norm(Q(:,m));
        Q(:,m) = Q(:,m)/R(m,m);
    end
        