function [x_cur,r,k] = GaussSeidel(A,b,fb)

n=size(A,1);
e = spalloc(3*n,1,3*n);
x_cur  = zeros(n,1);
x_prev = ones(n,1);
k    = 1;
e(k) = norm(b);
r(k) = e(k);

% backward Gauss-Seidel
if fb=='ba',
   disp('Backward Gauss-Seidel')
   C=triu(A,0);
   R=tril(A,-1);
else
   disp('Forward Gauss-Seidel')
   C=tril(A,0);
   R=triu(A,1);
end 

while (e(k) > 1e-6),
      x_next = C\(b-R*x_cur);
      x_prev=x_cur;
      x_cur=x_next;
      k=k+1;
      e(k)=norm(x_cur-x_prev);
%      r(k)=norm(b-A*x_cur);
end
k=k-1;
disp(['Convergence ' num2str(1e-6) ' reached for ' int2str(k) ' iterations.'])
