function [Ak,bk,ck] = IrkaN(A,b,c,k, tol);
   n = length(b);
   I = eye(n);
   itmax = 100; 
   [V,H,f] = ArnoldiC(inv(A),k,b);
   W = V;
   Ak = V'*A*V;
   Ak_eig = eig(Ak);
   Ak_eigOld = Ak_eig;
   Dist = 1;
   itcount = 1;
   while (Dist > tol & itcount < itmax),
      for j = 1:k,
          sigmaj = -Ak_eig(j);
          v = (A - sigmaj*I)\b;
          V(:,j) = v/norm(v);
          w = (A' - sigmaj*I)\c';
          W(:,j) = w/norm(w); 
      end
      [Z,S,Y] = svd(W'*V);
      S = sqrt(S);
      W = W*Z/S;  V = V*Y/S;
      Ak = W'*A*V; bk = W'*b; ck = c*V;
      Ak_eig = eig(Ak);
      [s,ii] = sort(real(Ak_eig));
      [s,jj] = sort(real(Ak_eigOld));
      SS = [abs(Ak_eig(ii)) abs(Ak_eigOld(jj))];
      Dist = norm(abs(Ak_eig(ii)) - abs(Ak_eigOld(jj)))  
      Ak_eigOld = Ak_eig;
      itcount = itcount + 1;
  end
  itcount = itcount
  SS = SS