function x = LSviaGS(A,b,itMax,x0)
% A is an m-by-n matrix and rank(A) = n.
% Assume that A is represented in sparse format.
% b is a column m-vector
% x is obtained by applying itMax steps of Gauss-Seidel to the linear system
% A'*A*x = A'*b with starting vector x = b.

% Apply the GS formula on p510 to system Fx = c:
% for i=1:n
% x(i) = (c(i) - (F(i,:)x - F(i,i)*x(i)))/(F(i,i)
% = x(i) + mu
% end
% where mu = (c(i) - F(i,:)x)/F(i,i)

% In our problem, F = A'A so F(i,:) = A(:,i)'*A
% and c = A'*b so c(i) = A(:,i)'b so
%
% for i=1:n
% mu = A(:,i)'*(b - A*x)/(A(:,i)'A(:,i));
% x(i) = x(i) + mu;
% end

% Note that r = b - A*x is the current residual. Instead of computing
% it from scratch we can update it:
%
% rNew = b - A*xNew = b - A(xOld + mu*I(:,i)) = rOld - mu*A(:,i)

[m,n] = size(A);
% Precompute the diagonal of A'*A
d = zeros(n,1);
for k=1:n
d(k) = A(:,k)'*A(:,k);
end

if nargin==3
x = zeros(n,1);
else
x = x0;
end
% Initialize the residual...
r = b - A*x;
for j=1:itMax
for i=1:n
mu = (A(:,i)'*r)/d(i);
x(i) = x(i) + mu;
r = r - mu*A(:,i);
% 2 length-m, sparse vector operations per step
end
end