function [x] = leastsq(A,b)
% returns a least squares solution of A x = b

[m,n] = size(A);
if (m == n)
% if we want the least norm least squares solution
% then just use x = pinv(A)*b;
  % Otherwise,. do it via QR decomposition
  [Q,R,E]=qr(full(A));
  % Determine r = effective rank
  tol = eps*norm(A,'fro');
  r = sum(abs(diag(R)) > tol);
  r = r(1); % fix for case where R is vector.
  x = Q'*b;
  y = [R(1:r,1:r)\x(1:r) ; zeros(length(b)-r,1)];
  x = E*y;
else
  x = A\b;
end;

return;

