function [x,J,it]=simplex(c,A,b,x0,J0)

 eps=1e-12;
 x=x0;
 J=J0;
 n=max(size(x));
 v1n=[1:n]';
 for it=0:inf
  Ji=v1n(J);
  Jci=v1n(~J);
  % A=A(:,J)*D;
  D=A(:,J)\A;
  cr=c(~J)-D(:,~J)'*c(J);
  tr=min(cr);
  if tr>-eps
   return;
  end
  r=Jci(min(find(cr==tr)));
  dr=D(:,r);
  if dr<=0
   fprintf(1,'Problem nach unten unbeschränkt\n!');
   return;
  end
  I=(dr>0);
  Ii=Ji(I);
  delta=min(x(Ii)./dr(I));
  nu=Ii(min(find(x(Ii)./dr(I)==delta)));
  x(r)=delta;
  x(J)=x(J)-delta*dr;
  J(nu)=false;
  J(r)=true;
 end