function [x,y,s]=ipmpk(A,b,c,epsilon)
% Praediktor-Korrektor-Verfahren von Mehrotra
% 
% Eingabe: A,b,c......Daten des gegebenen linearen Programms
%          epsilon....gewuenschte Loesungsgenauigkeit

% Schrittweitensteuerung:
eta=0.9995;

[m n]=size(A);
bnorm=norm(b);
cnorm=norm(c);
At=A';

% Startpunkt:
x=ones(n,1);
s=x;
y=ones(m,1);

mu=x'*s/n;

while norm(A*x-b)/max(1,bnorm)+norm(At*y+s-c)/max(1,cnorm)+abs(c'*x-b'*y)/max([1,abs(c'*x),abs(b'*y)])>epsilon
    
    % Praediktor-Schritt
    rhs1=c-At*y-s;
    rhs2=b-A*x;
    rhs3=-x.*s;
    D=diag(x./s);
    dy=(A*D*At)\(rhs2-A*(rhs3./s)+A*D*rhs1);
    ds=rhs1-At*dy;
    dx=(rhs3-x.*ds)./s;
    
    Ixp=find(dx<0);
    Isp=find(ds<0);
    tpP=min([1;-x(Ixp)./dx(Ixp)]);
    tdP=min([1;-s(Isp)./ds(Isp)]);

    muP=(x+tpP*dx)'*(s+tdP*ds)/n;
    sigma=(muP/mu)^3;
    
    % Korrektor-Schritt
    rhs3=rhs3-dx.*ds+sigma*mu;
    
    dyK=(A*D*At)\(rhs2-A*(rhs3./s)+A*D*rhs1);
    dsK=rhs1-At*dyK;
    dxK=(rhs3-x.*dsK)./s;
    
    IxKp=find(dxK<0);
    IsKp=find(dsK<0);
    tpK=min([1;-x(IxKp)./dxK(IxKp)]);
    tdK=min([1;-s(IsKp)./dsK(IsKp)]);
    
    tp=min(1,eta*tpK);
    td=min(1,eta*tdK);
    
    % Update-Schritt
    x=x+tp*dxK;
    y=y+td*dyK;
    s=s+td*dsK;
    
    mu=x'*s/n;
end