function [d,R,y]=FastGivens(A,b)
% FASTGIVENS reduce linear system to upper triangular form. 
%   [d,R,y]=FastGivens(A,b) reduces Ax=B using fast Givens reflections
%   to Rx=y. d contains the scaling factors. If [q,r]=qr(A) then
%   abs(diag(d)^(-0.5)*R)=abs(r)

[m n]=size(A); d=ones(m,1);
for i=1:n,
  for k=i+1:m,
    if A(k,i)~=0,
      alpha=-A(i,i)/A(k,i); beta=-alpha*d(k)/d(i);
      gamma=-alpha*beta;
      if gamma<=1,
        A(i,i)=A(i,i)*beta+A(k,i);
        if i<n,
          h=A(i,i+1:n)*beta+A(k,i+1:n);
          A(k,i+1:n)=A(i,i+1:n)+A(k,i+1:n)*alpha;
          A(i,i+1:n)=h;
        end;
        h=b(i)*beta+b(k); b(k)=b(i)+b(k)*alpha; b(i)=h;
        h=d(i); d(i)=d(k); d(k)=h;               % swap scaling factors
      else
        alpha=-1/alpha; beta=1/beta; gamma=1/gamma;
        A(i,i)=A(i,i)+beta*A(k,i);
        if i<n,
          h=A(i,i+1:n)+A(k,i+1:n)*beta;
          A(k,i+1:n)=A(i,i+1:n)*alpha -A(k,i+1:n);
          A(i,i+1:n)=h;
        end;
        h=b(i)+b(k)*beta; b(k)=b(i)*alpha-b(k); b(i)=h;
      end;                                  
      d(i)=d(i)*(1+gamma); d(k)=d(k)*(1+gamma);  % update for both cases
    end;                                    
  end                                       
end                                         
R=triu(A); y=b;
