function [q,r] = qr_(A)
% Function for computing QR factorization of the square matrix A
N = size(A,1);
% for each column of the source matrix...
for i = 1:N
% normalize column of A and put it into column of Q
r(i,i) = norm(A(:,i));
% what if r(i,i) == 0 ???
q(:,i) = A(:,i)./r(i,i);
% for each column to the right of the current column...
for j = i+1:N
r(i,j) = A(:,j)'*q(:,i);
%orthogonalize A(:,j) with respect to the vector Q(:,i) and put it back to matrix A
A(:,j) = A(:,j) - r(i,j)*q(:,i);
end;
end;