clc

n = 6;
m = 10;

A = hilb(m);
A = pascal(m);
A = A(1:m,1:n);

% CGS
Q = zeros(m,n);
R = zeros(n);
for k = 1:n
    Q(:,k) = A(:,k);
    for j = 1:k-1
        R(j,k) = Q(:,j)'*A(:,k);
        Q(:,k) = Q(:,k) - R(j,k)*Q(:,j);
    end
    R(k,k) = sqrt(Q(:,k)'*Q(:,k));
    Q(:,k) = Q(:,k)/R(k,k);
end

Q'*Q
norm(Q'*Q-eye(n), inf)

% MGS
Q = A;
R = zeros(n);
for k = 1:n
    R(k,k) = sqrt(Q(:,k)'*Q(:,k));
    Q(:,k) = Q(:,k)/R(k,k);
    for j = k+1:n
        R(k,j) = Q(:,k)'*Q(:,j);
        Q(:,j) = Q(:,j) - R(k,j)*Q(:,k);
    end
end

Q'*Q
norm(Q'*Q-eye(n), inf)

