rand("state", 1)
A = rand(9, 6);

[w k] = size(A);
for j = 1:k
  y(:,j) = A(:,j) ;

  for jj = 1:(j-1)
    y(:,j) = y(:,j) - dot(A(:,j),y(:,jj)) * y(:,jj)/ sum(y(:,jj).^2);
  end

  Q(:,j) = y(:,j) / sqrt(sum(y(:,j).^2));
end

for i = 1:k
  for j = i:k
    R(i,j) = Q(:,i)' * A(:,j);
  end
end
