% PATH - mixed complementarity problems by the PATH method.
% Package written in C-language.  Sparsity is preserved.
%
% 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
%
% Calling syntax: z = pathsol(z,l,u,cpfunjac)
%
%        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 = pathsol(z,l,u,cpfj)
%    ( M-file &  mex-interface written Oct 1995 by Michael Ferris)
%

[n,m]=size(z);
if min(m,n) ~= 1
  error('Input point is not a vector');
end
if m > n 
  z = z'; n = m;
end
Big = 1e20;

if nargin < 4
  cpfj = 'cpfunjac';
  if nargin < 3
    u = Big*ones(size(z));
    if nargin < 2
      l = zeros(size(z));
    end
  end
end

[m1,n1]=size(l);
if(n1 ~= 1)
  l = l';
  if(n1 ~= n | m1 ~=1)
    error('Input arguments are of incompatible sizes');
  end 
elseif (m1 ~= n)
  error('Input arguments are of incompatible sizes');
end
l = max(l,-Big*ones(size(l)));
  
[m1,n1]=size(u);
if(n1 ~= 1)
  u = u';
  if(n1 ~= n | m1 ~=1)
    error('Input arguments are of incompatible sizes');
  end 
elseif (m1 ~= n)
  error('Input arguments are of incompatible sizes');
end 
u = min(u,Big*ones(size(u)));

J = sparse(1,1);
eval(['[f,J,domerr] = ' cpfj '(z,1);']);
if ~issparse(J)
  error([cpfj ' must return a sparse Jacobian']);
end
nnzJ = nzmax(J)+n;
  
clear pathmex
status = pathmex(n, nnzJ, z, l, u, cpfj);
if (status ~= 1) 
  error('Path fails to solve problem');
end
return;
