function A=GivensQR(A);
% GIVENSQR Computes the QR-decomposition using Givens reflections
%   A=GivensQR(A); computes the QR-decomposition of the matrix A using
%   Givens reflections; after the decomposition, A contains the
%   implicit QR-decomposition of A. Instead of the zeros the
%   reflections are stored following a proposal of G. W. Stewart
%   (cf. srotg in BLAS). The decomposition is computed column wise.

[m n]=size(A);
for i=1:n,
  for k=i+1:m,
    if A(k,i)==0,                             % Compute  co and si
      co=1; si=0; z=0;
    else if abs(A(k,i))>=abs(A(i,i)),
      h=A(i,i)/A(k,i);                        % cot
      si=1/sqrt(1+h*h); co=si*h;
      if co~=0, z=1/co; else z=1; end
      else
        h=A(k,i)/A(i,i);                      % tan
        co=1/sqrt(1+h*h); si=co*h; z=si;
      end;
    end;
    A(i,i)=A(i,i)*co+A(k,i)*si;
    A(k,i)=z;                                 % store co or si in A(k,i) 
    if (si~=0)&(i<n),                         % Apply reflection
      S=[co,si;si,-co];
      A(i:k-i:k,i+1:n)=S*A(i:k-i:k,i+1:n);
    end;
  end
end
