function xLS = OutOfCoreLS(A,b)
% A is mxn with full column rank and in sparse format.
% Assume m>n
% b is mx1.
% xLS minimizes norm(Ax-b)

[m,n] = size(A);
% Initialization...
[Q,C] = qr(full(A(1:n+1,:)));
%b(1:n+1) = Q(:,1:n+1)'*b(1:n+1);
b(1:n+1) = Q'*b(1:n+1); % AH
% Now "read in" rows n+2 through n, updating R and b all along the way...
for k=n+2:m
    % Bring in row k...
    C(n+1,:) = A(k,:);
    % The C update process...
    for j=1:n
      % Determine c and s so [c s;-s c]'*[C(j,j);C(n+1,j)] = [*;0]
      [c,s] = Givens(C(j,j),C(n+1,j));
      C([j n+1],j:n) = [c s;-s c]'*C([j n+1],j:n);
      % Must update b accordingly...
      b([j k]) = [c s;-s c]'*b([j k]);
    end
end
xLS = C(1:n,1:n)\b(1:n);

% -2 for not applying correctly to b and/or C
% -1 for updating entire row of C each time through j loop above
% -1 for not setting C equal to R
% -2 for multiplying large matrices in performing Givens rotation