%-----------------------------------------------------------------------
% function [Q,R] = mgs(A)
% computes QR factorization via modified Gram-Schmidt
%-----------------------------------------------------------------------
function [Q,R,a,b] = mgs(A)
flops(0);tic;
% determine size
[m,n] = size(A);
% Initialize R matrix
d=zeros(n,1);
R = zeros(n,n);
for i=1:n-1
% use 1st of remaining non-orthogonal set as the q_i
% could also avoid dividing by norm, since we can
% simply divide by norm^2 after projecting (numerically
% harmless except for special cases or in symbolic math).
   R(i,i)=norm(A(:,i));
%% nn = A(:,i)'*A(:,i);
   A(:,i) = A(:,i)/R(i,i);
   for k=i+1:n
      R(i,k)=dot(A(:,i),A(:,k));
      A(:,k)=A(:,k)-A(:,i)*R(i,k);
   end
%  R(i,i+1:n)=A(:,i)'*A(:,i+1:n);
%  A(:,i+1:n)=A(:,i+1:n)-A(:,i)*R(i,i+1:n);
%% A(:,i+1:n)=A(:,i+1:n)-A(:,i)*R(i,i+1:n)/nn;
%% R(i,i) = sqrt(nn);
%% R(i,i+1:n)=R(i,i+1:n)/R(i,i);
end
R(n,n)=norm(A(:,n));
A(:,n)=A(:,n)/R(n,n);
Q=A;
a=flops; b=toc;
%-----------------------------------------------------------------------
