  function [xB,rhoxB,r] = LSviaQRP(A,b,tol)
% A is an mxn matrix (m>=n) and b is a column
% Uses the QR-with-column-pivoting factorization to solve the 
% (nearly) rank-deficient least squares problem min || Ax - b ||.
% tol is a small positive number that "defines" zero. In particular,
% if a floating point number f satisfies |f| <= tol, then it is treated
% as if it were zero. 
% xB is the minimizing x and rhoxB = || A*xB - b ||.
% r is an estimate of A's rank given the tolerance tol.

[m,n] = size(A);
[Q,R,piv] = qr(A,0);
% This means that AP = QR where
%      Q is mxn with orthonormal columns
%      R is nxn and upper triangular
%      P = I(:,piv)
% Look for the first negligble diagonal entry in R
k = 0;
while k<n && abs(R(k+1,k+1)) > tol
    k = k+1;
end
% r is an estimate of A's rank where
r = k;
% Note that if [Q,R,P] = qr(A), then 
%  || Ax - b ||^2 =  || (Q'AP)(P'x) - Q'b ||^2 
%                 =  || Rz - d ||^2
% where z = P'x and d = Q'b. Since the upper triangular matrix has the
% property that R(r+1:m) = 0 it follows that
%
%  || Ax - b ||^2 = || R(1:r,:)z - d(1:r)||^2 + || d(r+1:m) ||^2
%
% To obtain the "basic" solution xB we set z(r+1:n) = 0 and solve 
% the nonsingular triangular system R(1:r,1:r)z(1:r) = d(1:r).
% We can do all this even if we have the thin qr(A,0) factorization...
z = [R(1:r,1:r)\(Q(:,1:r)'*b); zeros(n-r,1)];
xB = zeros(n,1);
xB(piv) = z;
rhoxB = norm(A*xB - b);



    

