% PATHMCP - mixed complementarity problems by the PATH method.
% Package written in C-language.  Sparsity is preserved.
%
% Calling syntax: z = pathmcp(z,l,u,cpfunjac)
%
% PATH solves an MCP of the form
%
%       l <= z <= u
%       where l_i < z_i < u_i  => F_i(z) = 0
%             l_i = z          => F(z) >= 0
%             u_i = z          => F(z) <= 0
%
%        cpfunjac is the name of the m-file for 
%        evaluating the function F and its Jacobian J (without .m-extension).
%
%        The following m-file must be supplied (where default name is
%        'cpfunjac.m' unless stated otherwise in the variable cpfunjac).
%
%        'cpfunjac.m' contains function [f,J,domerr]=cpfunjac(z,jacflag) 
%                     that computes the function F 
%                     and if jacflag=1 the sparse Jacobian J at the point z.
%                     domerr returns the number of domain violations.
%
% If 'l','u' are omitted, z is assumed non-negative.
%
function z = pathmcp(z,l,u,cpfj)
%    ( M-file &  mex-interface written Oct 1995 by Michael Ferris)
%      updated Oct 1997 for Matlab 5.1 and Path 4.0
%

z = z(:); n = length(z);
Big = 1e20;

if nargin < 4
  cpfj = 'cpfunjac';
  if nargin < 3
    u = Big*ones(size(z));
    if nargin < 2
      l = zeros(size(z));
      if nargin < 1
	error('must have one input argument for pathmcp');
      end
    end
  end
end

l = l(:); u = u(:);
if(length(l) ~= n | length(u) ~= n)
  error('Input arguments are of incompatible sizes');
end 
l = max(l,-Big*ones(size(l)));
u = min(u,Big*ones(size(u)));

% this is a fix, nnz may be bigger than this
[f,J,domerr] = feval(cpfj,z+1e-5*ones(size(z))+1e-5*abs(z),1);
if ~issparse(J)
  error([cpfj ' must return a sparse Jacobian']);
end
if (domerr > 0)
  [f,J,domerr] = feval(cpfj,z,1);
end 
nnzJ = nzmax(J);
  
[status,ttime] = mcppath(n, nnzJ, z, l, u, cpfj);
if (status ~= 1) 
  status,
  error('Path fails to solve problem');
end
return;
