% INPUT   : A - MxN-matrix in sparse form
%           y - (M,1)-vector containing the data
%           x0- (N,1)-vector, the initial guess
%           w - relaxation parameter, 0<w<2
%           n - integer, the number of desired iteration rounds
% OUTPUT  : x - (N,1)-vector, the iterated reconstruction
%-----------------------------------------------------------------------------
% CALLS TO : None
% 29/2/96 Erkki Somersalo
%-----------------------------------------------------------------------------

function x = ART(A,y,x0,w,n)
 
j = 0;
x = x0;
M = length(y);
while j < n
  for ell = 1:M
    a = A(ell,:); 
    aa = full(a*a');
    if aa ~= 0
      x = x + w/aa*(y(ell) - a*x)*a';
    end
  end
  j = j+1;
end