DN1

function [X,ctr_j,ctr_g,ctr_s] = iter(A,b,x0,nat,omega)
  %
  % if X==x0 and ctr_j=ctr_g=ctr_s=0
  %   failed
  %

  N = length(A);
  j_ctr = 0;
  g_ctr = 0;
  s_ctr = 0;
  X = x0;

  % pd detection
  B = abs(A)';

  % required, not sufifcient
  d = 2*max(B) < sum(B);
  if(find(d))
    % cannot build p-d
    return;
  end

  % sufficient: every max alament in its column
  [_, perm] = max(B);
  if length(unique(perm)) ~= N
    % cannot build p-d
    return;
  end

  % tests done -> permutate matrix A for it to be diagonally dominant
  A = A(:, perm);
  Ar = A - eye(N).*A;
  Ad = diag(1./diag(A));

  % for jacobi and g-s
  T = -Ad*Ar;
  C = Ad*b;

  % jacobi
  ctr_j = 0;
  x = x0;
  while(max(abs(A*x-b)) > nat)
    x = T*x + C;
    ctr_j = ctr_j+1;
  end
  X(perm) = x;

  % GS
  ctr_g = 0;
  x = x0;
  while(max(abs(A*x-b)) > nat)
    for j = [1:N]
      x(j) = T(j,:)*x + C(j);
    end
    ctr_g = ctr_g+1;
  end

  % sor
  ctr_s = 0;
  x = x0;
  while(max(abs(A*x-b)) > nat)
    for j = [1:N]
      resid = Ad(j,j) * (b(j) - A(j,:)*x);
      x(j) = x(j) + omega*resid;
    end
    ctr_s = ctr_s+1;
  end
end