wolffd@0: function [z, w, status, z0] = lcprog(M, q, varargin)
wolffd@0: %LCPROG - Linear complementary programming with Lemke's algorithm.
wolffd@0: %
wolffd@0: %	[z, w] = lcprog(M, q);
wolffd@0: %	[z, w, status, z0] = lcprog(M, q, varargin);
wolffd@0: %
wolffd@0: % Inputs:
wolffd@0: %	M - Matrix (n x n matrix)
wolffd@0: %	q - Vector (n vector)
wolffd@0: %
wolffd@0: % Optional Inputs:
wolffd@0: %	Optional name-value pairs are:
wolffd@0: %	'maxiter'	Maximum number of iterations before termination of algorithm (default is (1+n)^3).
wolffd@0: %	'tiebreak'	Method to break ties for pivot selection (default is 'first'). Options are:
wolffd@0: %				'first'		Select pivot with lowest index.
wolffd@0: %				'last'		Select pivot with highest index.
wolffd@0: %				'random'	Select a pivot at random.
wolffd@0: %	'tolpiv'	Pivot tolerance below which a number is considered to be zero (default is 1.0e-9).
wolffd@0: %
wolffd@0: % Outputs:
wolffd@0: %	z - Solution vector (n vector)
wolffd@0: %	w - Complementary solution vector (n vector)
wolffd@0: %	status - Status of optimization with values (integer)
wolffd@0: %		 6		Ray termination
wolffd@0: %		 1		Convergence to a solution (normal termination)
wolffd@0: %		 0		Termination prior to convergence
wolffd@0: %		-2		Infeasible problem
wolffd@0: %		-3		Bad pivot to an infeasible solution
wolffd@0: %	z0 - Artificial variable which is 0 upon normal termination (scalar)
wolffd@0: %
wolffd@0: % Comments:
wolffd@0: %	Given a matrix M and a vector q, the linear complementary programming (LCP) problem seeks an
wolffd@0: %	orthogonal pair of non-negative vectors w and z such that
wolffd@0: %		w = M * z + q
wolffd@0: %	with
wolffd@0: %		w' * z = 0
wolffd@0: %		w >= 0
wolffd@0: %		z >= 0 .
wolffd@0: %
wolffd@0: %	The algorithm starts with the problem
wolffd@0: %		w = M * z + q + z0 * 1
wolffd@0: %	with an artificial variable z0 > 0 such that q + z0*1 >= 0. The algorithm goes through a series
wolffd@0: %	of pivot steps that seeks to drive z0 to 0. The algorithm can stop in two ways with either
wolffd@0: %	normal termination or ray termination.
wolffd@0: %
wolffd@0: %	With normal termination, the artificial variable z0 is zero so that the original problem is
wolffd@0: %	solved.
wolffd@0: %
wolffd@0: %	With ray termination, the artificial variable z0 is non-negative with "ray" solution in the form
wolffd@0: %		w = wR + lambda * dw
wolffd@0: %		z = zR + lambda * dz
wolffd@0: %		z0 = z0R + lambda * dz0
wolffd@0: %	that satisfies
wolffd@0: %		w = M * z + q + z0 * 1
wolffd@0: %		w' * z = 0
wolffd@0: %	for any lambda >= 0. If ray termination occurs, the returned variables w, z, and z0 are the
wolffd@0: %	points on the ray with lambda = 0.
wolffd@0: %
wolffd@0: %	If tolpiv is too small, the algorithm may identify a pivot that is effectively zero. If this
wolffd@0: %	situation occurs, the status code -3 is triggered. Try a larger pivot tolerance 'tolpiv' although
wolffd@0: %	a pivot that is too large can trigger false convergence due to failure to find a pivot.
wolffd@0: %
wolffd@0: 
wolffd@0: %
wolffd@0: % Copyright 2007 The MathWorks, Inc.
wolffd@0: % $Revision: 1.1.6.4 $   $ Date: $
wolffd@0: 
wolffd@0: % initialization
wolffd@0: 
wolffd@0: if nargin < 2 || isempty(M) || isempty(q)
wolffd@0:    error('finance:lcprog:MissingInputArgument', ...
wolffd@0:       'One or more missing required input arguments M, q.');
wolffd@0: end
wolffd@0: 
wolffd@0: q = q(:);
wolffd@0: n = size(q,1);
wolffd@0: 
wolffd@0: % process name-value pairs (if any)
wolffd@0: 
wolffd@0: if mod(nargin-2,2) ~= 0
wolffd@0:    error('finance:lcprog:InvalidNameValuePair', ...
wolffd@0:       'Invalid name-value pairs with either a missing name or value.');
wolffd@0: end
wolffd@0: 
wolffd@0: names = { 'maxiter', 'tiebreak', 'tolpiv' };				% names
wolffd@0: values = { (1 + n)^3, 'first', 1.0e-9 };					% default values
wolffd@0: try
wolffd@0:    [maxiter, tiebreak, tolpiv] = parsepvpairs(names, values, varargin{:});
wolffd@0: 
wolffd@0: catch E
wolffd@0:    E.throw
wolffd@0: end
wolffd@0: 
wolffd@0: % check optional arguments
wolffd@0: 
wolffd@0: if maxiter <= 0
wolffd@0:    error('finance:lcprog:NonPositiveInteger', ...
wolffd@0:       'Maximum number of iterations (maxiter) must be a positive integer.');
wolffd@0: end
wolffd@0: 
wolffd@0: if ~isempty(strmatch(tiebreak,{'default','first'}))
wolffd@0:    tb = 1;
wolffd@0: elseif ~isempty(strmatch(tiebreak,'last'))
wolffd@0:    tb = 2;
wolffd@0: elseif ~isempty(strmatch(tiebreak,'random'))
wolffd@0:    tb = 3;
wolffd@0: else
wolffd@0:    tb = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: if tolpiv < 2*eps
wolffd@0:    error('finance:lcprog:InvalidTolerance', ...
wolffd@0:       'Unrealistically small tolerance (tolpiv) specified.');
wolffd@0: end
wolffd@0: 
wolffd@0: % trivial solution
wolffd@0: 
wolffd@0: if all(q >= 0)
wolffd@0:    z = zeros(n,1);
wolffd@0:    w = q;
wolffd@0:    status = 1;
wolffd@0:    z0 = 0;
wolffd@0:    return
wolffd@0: end
wolffd@0: 
wolffd@0: % set up general problem
wolffd@0: 
wolffd@0: w = zeros(n,1);
wolffd@0: z = zeros(n,1);
wolffd@0: u = ones(n,1);
wolffd@0: z0 = max(-q);
wolffd@0: 
wolffd@0: tableau = [ eye(n), -M, -u, q];
wolffd@0: basic = (1:n)';
wolffd@0: 
wolffd@0: % initial pivot and initial basis
wolffd@0: 
wolffd@0: r = find(q == min(q),1);			% r are pivot rows
wolffd@0: cx = r;								% w(r) leaves the basis
wolffd@0: c = 2*n + 1;						% initial pivot column
wolffd@0: basic(r) = c;						% z0 enters the basis
wolffd@0: 
wolffd@0: % pivot step
wolffd@0: 
wolffd@0: pr = (1/tableau(r,c)) * tableau(r,:);
wolffd@0: pc = tableau(:,c);
wolffd@0: tableau = tableau - pc * pr;
wolffd@0: tableau(r,:) = pr;
wolffd@0: 
wolffd@0: % main loop
wolffd@0: 
wolffd@0: for iter = 2:maxiter
wolffd@0: 
wolffd@0:    % get new basic variable
wolffd@0: 
wolffd@0:    if cx <= n
wolffd@0:       c = n + cx;
wolffd@0:    elseif cx <= 2*n
wolffd@0:       c = cx - n;
wolffd@0:    end
wolffd@0: 
wolffd@0:    % do ratio test
wolffd@0: 
wolffd@0:    d = tableau(:,c);
wolffd@0:    mind = NaN(n,1);
wolffd@0:    for i = 1:n
wolffd@0:       if d(i) > tolpiv
wolffd@0:          mind(i) = tableau(i,end)/d(i);
wolffd@0:       end
wolffd@0:    end
wolffd@0:    iind = find(mind == nanmin(mind));
wolffd@0: 
wolffd@0:    % ray termination
wolffd@0: 
wolffd@0:    if isempty(iind)						% cannot find a new basic variable
wolffd@0:       for i = 1:n
wolffd@0:          if basic(i) <= n
wolffd@0:             ii = basic(i);
wolffd@0:             w(ii) = tableau(i,end);
wolffd@0:             z(ii) = 0;
wolffd@0:          elseif basic(i) <= 2*n
wolffd@0:             ii = basic(i) - n;
wolffd@0:             w(ii) = 0;
wolffd@0:             z(ii) = tableau(i,end);
wolffd@0:          else
wolffd@0:             z0 = tableau(i,end);
wolffd@0:          end
wolffd@0:       end
wolffd@0:       if lcprealitycheck(M, q, w, z)
wolffd@0:          status = -3;
wolffd@0:          z = NaN(n,1);
wolffd@0:          w = NaN(n,1);
wolffd@0:          z0 = NaN;
wolffd@0:       else
wolffd@0:          status = 6;
wolffd@0:       end
wolffd@0:       return
wolffd@0:    end
wolffd@0: 
wolffd@0:    % identify next basic variable to be replaced using ratio test
wolffd@0: 
wolffd@0:    if tb == 3
wolffd@0:       ii = ceil(rand()*numel(iind));
wolffd@0:       if ii < 1
wolffd@0:          ii = 1;
wolffd@0:       end
wolffd@0:       r = iind(ii);
wolffd@0:    elseif tb == 2
wolffd@0:       r = iind(end);
wolffd@0:    elseif tb == 1
wolffd@0:       r = iind(1);
wolffd@0:    end
wolffd@0:    for i = 1:numel(iind)					% always select z0 if it is in the index set
wolffd@0:       if basic(iind(i)) == (2*n + 1)
wolffd@0:          r = iind(i);
wolffd@0:       end
wolffd@0:    end
wolffd@0: 
wolffd@0:    % new basis
wolffd@0: 
wolffd@0:    cx = basic(r);							% get variable leaving the basis
wolffd@0:    basic(r) = c;							% move new basic variable into the basis
wolffd@0: 
wolffd@0:    % pivot step
wolffd@0: 
wolffd@0:    pr = (1/tableau(r,c)) * tableau(r,:);
wolffd@0:    pc = tableau(:,c);
wolffd@0:    tableau = tableau - pc * pr;
wolffd@0:    tableau(r,:) = pr;
wolffd@0: 
wolffd@0:    if ~isfinite(tableau(1,end))			% trap overflow/underflow
wolffd@0:       warning('finance:lcprog:InfeasibleProblem', ...
wolffd@0:          'Tableau overflow/underflow due to bad pivot.');
wolffd@0:       status = -3;
wolffd@0:       z = NaN(n,1);
wolffd@0:       w = NaN(n,1);
wolffd@0:       z0 = NaN;
wolffd@0:       return
wolffd@0:    end
wolffd@0: 
wolffd@0:    % convergence test
wolffd@0: 
wolffd@0:    if all(basic <= 2*n)					% test if z0 left the basis -> convergence
wolffd@0:       for i = 1:n
wolffd@0:          if basic(i) <= n
wolffd@0:             ii = basic(i);
wolffd@0:             w(ii) = tableau(i,end);
wolffd@0:             z(ii) = 0;
wolffd@0:          else
wolffd@0:             ii = basic(i) - n;
wolffd@0:             w(ii) = 0;
wolffd@0:             z(ii) = tableau(i,end);
wolffd@0:          end
wolffd@0:       end
wolffd@0:       z0 = 0;
wolffd@0:       if lcprealitycheck(M, q, w, z)
wolffd@0:          status = -3;
wolffd@0:          z = NaN(n,1);
wolffd@0:          w = NaN(n,1);
wolffd@0:          z0 = NaN;
wolffd@0:       else
wolffd@0:          status = 1;
wolffd@0:       end
wolffd@0:       return
wolffd@0:    end
wolffd@0: end
wolffd@0: 
wolffd@0: % too many iterations
wolffd@0: 
wolffd@0: for i = 1:n
wolffd@0:    if basic(i) <= n
wolffd@0:       ii = basic(i);
wolffd@0:       w(ii) = tableau(i,end);
wolffd@0:       z(ii) = 0;
wolffd@0:    elseif basic(i) <= 2*n
wolffd@0:       ii = basic(i) - n;
wolffd@0:       w(ii) = 0;
wolffd@0:       z(ii) = tableau(i,end);
wolffd@0:    else
wolffd@0:       z0 = tableau(i,end);
wolffd@0:    end
wolffd@0: end
wolffd@0: if lcprealitycheck(M, q, w, z)
wolffd@0:    status = -3;
wolffd@0:    z = NaN(n,1);
wolffd@0:    w = NaN(n,1);
wolffd@0:    z0 = NaN;
wolffd@0: else
wolffd@0:    status = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function notok = lcprealitycheck(M, q, w, z)
wolffd@0: %LCPREALITYCHECK - Make sure candidate solution z and w is complementary basic feasible.
wolffd@0: 
wolffd@0: u = w - M*z - q;
wolffd@0: 
wolffd@0: if (max(u) - min(u)) > 1
wolffd@0:    warning('finance:lcprog:InfeasibleProblem', ...
wolffd@0:       'Candidate solution is infeasible due to a bad pivot.');
wolffd@0:    notok = true;
wolffd@0: else
wolffd@0:    notok = false;
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % [EOF]