1 LPcode of Matlab 2006-05-11 1 ACcode of Matlab 2006-05-12

pdnewton.m

if (iter == 0)

    [m,n]  = size(A);

    x=ones(n,1);

    y=0*ones(m,1);

    c=2*ones(n+1,1)+rand(n+1,1);

    b=A*x;

    A=[A -b];

    x(n+1)=1;

    n=n+1;

    s=c-A'*y;

    x= (n/(c'*x))*x;

    rk     = ones(n,1)-x.*s;

    norm(rk)

end;

Ak     = A*diag(x./s);

Mk     = Ak*A';

dy     = -Mk\(A*(rk./s));

ds     = -A'*dy;

dx     = rk./s - (x.*ds)./s;

y      = y+dy;

s      = s+ds;

x      = x+dx;

rk     = ones(n,1)-x.*s;

norm(rk)

iter   = iter+1

===============================================================================

pdpotentiallp.m

if (iter == 0)

    [m,n]  = size(A);

    x=ones(n,1);

    y=0*ones(m,1);

    c=2*ones(n,1)+rand(n,1);

    b=A*x;

    s=c-A'*y;

    potential=(n+rho)*log(x'*s)-ones(n,1)'*log(x.*s)

end;

mu     = x'*s/n;

gamma  = n/(n+rho);

rk     = gamma*mu*ones(n,1)-x.*s;

Ak     = A*diag(x./s);

Mk     = Ak*A';

dy     = -Mk\(A*(rk./s));

ds     = -A'*dy;

dx     = rk./s - (x.*ds)./s;

theta  = min(x.*s)/((x.*rk)'*(rk./s));

theta  = alpha*sqrt(theta);

y      = y+theta*dy;

s      = s+theta*ds;

x      = x+theta*dx;

potential=(n+rho)*log(x'*s)-ones(n,1)'*log(x.*s)

iter   = iter+1