   function [xLS,res] = LSq(A,b)
%
%  Pre:
%    A     m-by-n and rank(A) = n.
%    b     m-by-1.
%  Post:
%    xLS   n-by-1, minimizes the 2-norm of Ax-b.
%    res   norm(A*xLS-b)
%

   [m,n] = size(A);
   for j=1:n
      for i=m:-1:j+1
         %Zero A(i,j)
         [c,s] = Rotate(A(i-1,j),A(i,j));
         A(i-1:i,j:n) = [c s; -s c]*A(i-1:i,j:n);
		 b(i-1:i) = [c s; -s c]*b(i-1:i);
      end
   end
   xLS = UTriSol(A(1:n,1:n),b(1:n));
   if m==n
      res = 0;
   else
      res = norm(b(n+1:m));
   end
