function [A,alpha] = QRDecomposition( A )
% QR - decomposition
% argument "A" is the matrix we want to decompose

% sign, being either +1 or -1, but not 0
function sig = pmsign (val)
    if (val < 0)
        sig = -1;
    else
        sig = +1;
    end
end

% auxiliary function: apply Q
function x = ApplyQ( k, alpha, v, x )
    s = 0;
    n = size(v);
    for l = k:n
        s = s + v(l)*x(l);
    end
    s = s / (alpha*v(k));
    for l=k:n
        x(l) = x(l) + s*v(l);
    end
end

% main program
% get the length of the vector
[m,n] = size( A );

for k = 1:(n-1)
    norm = 0;
    for i =k:n
        norm = norm + A(i,k)*A(i,k);  
    end
    norm = sqrt(norm);
    alpha(k) = -pmsign(A(k,k)) * norm;
    A(k,k) = A(k,k) - alpha(k);
    for j = k+1:n
       A(:,j) = ApplyQ( k, alpha(k), A(:,k), A(:,j) );
    end
end

end
