function [x_out, fval, it_out, ll_out] = MPW_ZPK(D, c, A, b, eps, M)
    MAX_IT = 100;
    delta = 0.2;
    beta = 0.999;
    it = 0;
    n = length(c);
    m = length(b);

    ex = ones(n,1);
    ey = ones(m,1);

    x_0 = ones(n,1);
    y_0 = ones(m,1);
    w0 = ones(m,1);
    z0 = ones(n,1);

    x_k = x_0;
    y_k = y_0;
    w_k = w0;
    z_k = z0;

    while(it < MAX_IT)
        it = it + 1;

        if abs(max(x_k))>M || abs(max(y_k)) > M
            x_out = -1;
            fval_out = -1;
            it_out = it;
            break
        end

        rho = b - A * x_k + w_k;
        sigma = c - A' * y_k - z_k + D * x_k;
        gamma = z_k' * x_k + y_k' * w_k;

        if sum(abs(rho))<eps || sum(abs(sigma)) < eps || abs(gamma) < eps
            x_out = x_k;
            ll_out = [y_k; z_k];
            fval_out = 1/2* x_k' * D * x_k + c'*x_k;
            it_out = it;
            break
        end

        r_k = delta * gamma / (n + m);

        X_inv = inv(diag(x_k));
        Y_inv = inv(diag(y_k));

        X = diag(x_k);
        Z = diag(z_k);
        Y = diag(y_k);
        W = diag(w_k);

        matrix = [-(X_inv*Z + D) A'; A Y_inv * W];
        vector = [c - A' * y_k - r_k * X_inv * ex + D * x_k; ...
             b - A * x_k + r_k * Y_inv * ey];
        result = matrix \ vector;

        delta_x = result(1 : n);
        delta_y = result(n+1 : n+m);
        delta_z = X_inv*(r_k*ex - X*Z*ex - Z*delta_x);
        delta_w = Y_inv*(r_k*ey - Y*W*ey - W*delta_y);

        indices = delta_w < 0;
        values =  -beta * w_k ./ delta_w;
        alpha1 = min([1; values(indices)]);

        indices = delta_y < 0;
        values = -beta * y_k ./ delta_y;
        alpha2 = min([1; values(indices)]);

        indices = delta_x < 0;
        values =  -beta * x_k ./ delta_x;
        alpha3 = min([1; values(indices)]);

        indices = delta_z < 0;
        values =  -beta * z_k ./ delta_z;
        alpha4 = min([1; values(indices)]);

        alpha = min([1 alpha1 alpha2 alpha3 alpha4]);

        x_k = x_k + alpha * delta_x;
        z_k = z_k + alpha * delta_z;
        y_k = y_k + alpha * delta_y;
        w_k = w_k + alpha * delta_w;

        delta = 0.2*delta;
    end

end