function [x, y, s, obj] = pdpc(A, b, c, x, y, s, eta)
%PDPC   Primal-dual predictor-corrector method
%
%       [x, y, s, obj] = pdpc(A, b, c, x, y, s, eta)
%
%       minimize        c'*x
%       subject to      Ax = b
%
%       A: m * n, b: m * 1, c: n * 1
%       x, y, s: initial values for the primal and dual
%                problems.
%       default eta is 0.99999 (0 < eta < 1)

%   Ling Li, Caltech
%   05/20/2001

[m, n] = size(A);

if nargin < 7,
    eta = 0.99999;
end

ite = 0;
while 1,
    mu = x' * s / n;
    if mu < 1e-15, break; end

    D = diag(x ./ s);   %% inv(S) * X;
    iADAT = inv(A * D * A');
    rP = b - A * x;
    rD = c - A' * y - s;

    %% Predictor step
    v = - x .* s;       %% - X S e

    iSv = v ./ s;       %% inv(S) * v;
    dy = iADAT * (rP - A * (iSv - D * rD));
    ds = rD - A' * dy;
    dx = iSv - D * ds;

    aprim = 1;
    neg_idx = find(dx < 0);
    if ~isempty(neg_idx),
        aprim = min([aprim; -x(neg_idx) ./ dx(neg_idx)]);
    end

    adual = 1;
    neg_idx = find(ds < 0);
    if ~isempty(neg_idx),
        adual = min([adual; -s(neg_idx) ./ ds(neg_idx)]);
    end

    mu_aff = (x + aprim*dx)' * (s + adual*ds) / n;
    theta = (mu_aff / mu) ^ 3;

    %% Corrector step
    v = v + theta * mu - dx.*ds;%% dX dS e

    iSv = v ./ s;
    dy = iADAT * (rP - A * (iSv - D * rD));
    ds = rD - A' * dy;
    dx = iSv - D * ds;

    %% Update
    aprim = 1;
    neg_idx = find(dx < 0);
    if ~isempty(neg_idx),
        aprim = min([aprim; - eta * x(neg_idx) ./ dx(neg_idx)]);
    end

    adual = 1;
    neg_idx = find(ds < 0);
    if ~isempty(neg_idx),
        adual = min([adual; - eta * s(neg_idx) ./ ds(neg_idx)]);
    end

    ite = ite + 1;
    y = y + adual * dy;
    x = x + aprim * dx;
    s = s + adual * ds;
    disp(['ite = ' int2str(ite) ', n mu = ' num2str(x' * s)]);
end

obj = c' * x;

disp(['x = [' num2str(x') ']''']);
disp(['s = [' num2str(s') ']''']);
disp(['y = [' num2str(y') ']''']);
disp(['c''*x = ' num2str(obj) ', b''*y = ' num2str(b'*y)]);
