[SOLVED] CS function [x,y,s] = ipdriver(c,A,b,tol,maxit,printlevel)

function [x,y,s] = ipdriver(c,A,b,tol,maxit,printlevel)
%IPDRIVERPrimal-dual interior-point method.
%
%This is the driver function of a very simple IPM for solving the
%linear programming problem
%
%min cx subject to A*x=b, x>=0,(1)
%
%where A has full row rank. In each iteration NEWTONSOLVE[R] is called
%for solving the linear system of step equations.
%
%[x,y,s] = ipdriver(c,A,b) finds the primal-dual solution x, (y,s) of
%the linear program in form (1). y are the Lagrange multipliers to the
%equality constraints and s>=0 the Lagrange multipliers to the
%nonnegativity constraints on x.
%
%[x,y,s] = ipdriver(c,A,b,tol) specifies the termination tolerance.
%The method terminates when the relative primal residual, the relative
%dual residual and the duality measure are all below tol.
%Default: 1e-7.
%
%[x,y,s] = ipdriver(c,A,b,tol,maxit) specifies the maximum number of
%iterations. Default: 100.
%
%[x,y,s] = ipdriver(c,A,b,tol,maxit,printlevel) sets the printlevel.
%0: turn off iteration output
%1: print primal and dual residual and duality measure
%2: print centering parameter and step length
%3: print residuals in the solution of the step equations
%Default: 1
%
%See ipexample.m for an illustrative example.
%
%Lukas Schork, 2015

[m n] = size(A);

% Make sure that b and c are column vectors of dimension m and n.
if (size(b,2) > 1) b = b; end
if (size(c,2) > 1) c = c; end
if (~isequal(size(c),[n,1]) || ~isequal(size(b),[m,1]))
error(problem dimension incorrect);
end

% Make sure that A is sparse and b, c are full.
if (~issparse(A)) A = sparse(A); end
if (issparse(b)) b = full(b); end
if (issparse(c)) c = full(c); end

% Set default values for missing parameters.
if (nargin < 4 || isempty(tol)) tol = 1e-7; endif (nargin < 5 || isempty(maxit)) maxit = 100; endif (nargin < 6 || isempty(printlevel)) printlevel = 1; endpl = printlevel;% ********************************************************************% Initialization% ==============%% Choose initial point (x,y,s). We make the very simple choice%% x(i)=1, y(j)=0, s(i)=1 for 1<=i<=n, 1<=j<=m.%% This point will be primal and dual infeasible, but it is perfectly% centered.% ********************************************************************x = ones(n,1);y = zeros(m,1);s = ones(n,1);% Set fixed parameters, see below.sigma1 = 0.5;sigmamin = 0.05;sigmamax = 0.95;% Initialize persistent variables.iter = 0;alpha_x = 0;alpha_s = 0;% ********************************************************************% IPM Main-Loop% =============%% In each step of the interior-point method we solve the Newton system%% [ A 00 ] (dx) (b-A*x)=: (r)% [ 0 A’ I ] (dy) = (c-A’*y-s) =: (p)% [ S 0X ] (ds) (simga*mu*e-X*S*e) =: (q),%% where mu=dot(x,s)/n and sigmain[sigmamin,sigmamax]. The next% iterate is given by%% x = x+alpha_x*dx, y = y+alpha_s*dy, s = s+alpha_s*ds,%% where alpha_x and alpha_s are the primal and dual step length.% ********************************************************************header(pl);while (iter < maxit)% Check for termination. We have found a sufficiently accurate% solution if the relative primal residual, the relative dual% residual and the duality measure are all below tol.r = b-A*x;p = c-A’*y-s;mu = dot(x,s)/n;if (norm(r)/norm(b) < tol && norm(p)/norm(c) < tol && mu < tol)fprintf(‘optimal solution found
‘);break;enditer = iter+1;% Choose the centering parameter sigma. In the first iteration we% use the default value sigma1. In all later iterations we choose% simga depending on the step lengths in the previous iteration.if (iter == 1)sigma = sigma1;elsesigma = max(1-alpha_x,1-alpha_s)^5;endsigma = min(sigma,sigmamax);sigma = max(sigma,sigmamin);q = sigma*mu*ones(n,1) – x.*s;% Solve the Newton system and calculate residuals.NS = newtonsolver(A,x,s);[dx,dy,ds] = newtonsolve(NS,r,p,q);res1 = norm(A*dx-r);res2 = norm(A’*dy+ds-p);res3 = norm(x.*ds+s.*dx-q);% Determine primal and dual step length. Calculate “step to the% boundary” alphamax_x and alphamax_s. Then choose 0 < rho < 1% heuristically and set step length = rho * step to the boundary.idx = dx < 0;ids = ds < 0;alphamax_x = min([1;-x(idx)./dx(idx)]);alphamax_s = min([1;-s(ids)./ds(ids)]);rho = 0.95;alpha_x = rho*alphamax_x;alpha_s = rho*alphamax_s;% Make step.x = x+alpha_x*dx; y = y+alpha_s*dy; s = s+alpha_s*ds;% Print iteration output.xinf = norm(A*x-b);sinf = norm(A’*y+s-c);mu = dot(x,s)/n;output(pl,iter,xinf,sinf,mu,sigma,alpha_x,alpha_s,res1,res2,res3);end % while (iter < maxiter)% The IPM has terminated either because the solution accuracy is% reached or the maximum number of iterations is exceeded. Print% result.fprintf(‘iterations: %4d
‘, iter);fprintf(‘primal feasibility: %8.2e
‘, norm(A*x-b));fprintf(‘dual feasibility: %8.2e
‘, norm(A’*y+s-c));fprintf(‘complementarity: %8.2e
‘, dot(x,s)/n);end% ======================================================================% header% ======================================================================function header(pl)if (pl >= 1)
fprintf( );
fprintf(%4s, iter);
fprintf(%8s, pr feas);
fprintf(%8s, dl feas);
fprintf(%8s, mu);
end
if (pl >= 2)
fprintf();
fprintf(%8s, sigma);
fprintf(%8s, alpha_x);
fprintf(%8s, alpha_s);
end
if (pl >= 3)
fprintf();
fprintf(%8s, res1);
fprintf(%8s, res2);
fprintf(%8s, res3);
end
if (pl >= 1)
fprintf(
============================);
end
if (pl >= 2)
fprintf(========================);
end
if (pl >= 3)
fprintf(========================);
end
if (pl >= 1) fprintf(
); end
end

% ======================================================================
% output
% ======================================================================

function output(pl,it,xinf,sinf,mu,sigma,alpha_x,alpha_s,res1,res2,res3)
if (pl >= 1)
fprintf( );
fprintf(%4d, it);
fprintf(%8.2e, xinf);
fprintf(%8.2e, sinf);
fprintf(%8.2e, mu);
end
if (pl >= 2)
fprintf();
fprintf(%8.2e, sigma);
fprintf(%8.2e, alpha_x);
fprintf(%8.2e, alpha_s);
end
if (pl >= 3)
fprintf();
fprintf(%8.2e, res1);
fprintf(%8.2e, res2);
fprintf(%8.2e, res3);
end
if (pl >= 1) fprintf(
); end
end