function [U,R]=gramschmidt(A,test)
% [U,R]=gramschmidt(A,test)
% computes the orthogonal decomposition 
%
% A=Q*R   , Q'*Q=I and of the same shape as A
% R upper triangular , with A n times m and m<=n
% using the modified gram schmidt orthogonalization.
% if test==1 then norm(U'*U-I) and norm(A-U*R) will
% be displayed
[n,m]=size(A);
U=A;
R=zeros(m,m);
for i=1:m
    h=norm(U(:,i));
    if h==0
        error('gram schmidt: matrix ist rangdefizient');
    end
    U(:,i)=U(:,i)/h;
    R(i,i)=h;
    for j=i+1:m
        R(i,j)=U(:,i)'*U(:,j);
        U(:,j)=U(:,j)-U(:,i)*R(i,j);
    end
end
if test==1 
disp('gramschmidt test:orthogonality');
norm(U'*U-eye(m))
disp('gramschmidt test: factorization');
norm(A-U*R)
end