From 450d7b099a10c5da1b9df01667a9c889c8b3606b Mon Sep 17 00:00:00 2001
From: Michel Juillard <michel.juillard@mjui.fr>
Date: Sun, 15 May 2016 17:33:59 +0200
Subject: [PATCH] fixing error status for PATH hook

---
 matlab/dynare_solve.m |   6 +-
 matlab/pathlcp.m      | 172 ++++++++++++++++++++++++++++++++++++
 matlab/pathmcp.m      | 197 ++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 373 insertions(+), 2 deletions(-)
 create mode 100644 matlab/pathlcp.m
 create mode 100644 matlab/pathmcp.m

diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 6540534fb..09afd3df7 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -206,7 +206,8 @@ elseif options.solve_algo == 11
     end
     olcppath = options.lcppath;
     [junk,M] = func(x,varargin{:});
-    [x,mu] = pathlcp(fjac,olcppath.q,olcppath.lb,olcppath.ub,x,olcppath.A,olcppath.b,olcppath.t,olcppath.mu0);
+    [x,mu,status] = pathlcp(fjac,olcppath.q,olcppath.lb,olcppath.ub,x,olcppath.A,olcppath.b,olcppath.t,olcppath.mu0);
+    info = ~status;
 elseif options.solve_algo == 12
     % PATH mixed complementary problem
     % PATH linear mixed complementary problem
@@ -219,7 +220,8 @@ elseif options.solve_algo == 12
     global mcp_data
     mcp_data.func = func;
     mcp_data.args = varargin;
-    [x,mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
+    [x,fval,jac,mu,status] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
+    info = ~status;
 else
     error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10:12]')
 end
diff --git a/matlab/pathlcp.m b/matlab/pathlcp.m
new file mode 100644
index 000000000..30ddc367e
--- /dev/null
+++ b/matlab/pathlcp.m
@@ -0,0 +1,172 @@
+function [z,mu,status] = pathlcp(M,q,l,u,z,A,b,t,mu)
+% pathlcp(M,q,l,u,z,A,b,t,mu)
+%
+% Solve the standard linear complementarity problem using PATH:
+%     z >= 0, Mz + q >= 0, z'*(Mz + q) = 0
+%
+% Required input:
+%     M(n,n)  - matrix
+%     q(n)    - vector
+% 
+% Output:
+%     z(n)    - solution
+%     mu(m)   - multipliers (if polyhedral constraints are present)
+%
+% Optional input:
+%     l(n)    - lower bounds                       default: zero
+%     u(n)    - upper bounds                       default: infinity
+%     z(n)    - starting point                     default: zero
+%     A(m,n)  - polyhedral constraint matrix       default: empty
+%     b(m)    - polyhedral right-hand side         default: empty
+%     t(m)    - type of polyhedral constraint      default: 1
+%                  < 0: less than or equal
+%                    0: equation
+%                  > 0: greater than or equal
+%     mu(m)   - starting value for multipliers     default: zero
+%
+% The optional lower and upper bounds are used to define a linear mixed 
+% complementarity problem (box constrained variational inequality).
+%       l <= z <= u
+%       where l_i < z_i < u_i  => (Mz + q)_i = 0
+%             l_i = z          => (Mz + q)_i >= 0
+%             u_i = z          => (Mz + q)_i <= 0
+% 
+% The optional constraints are used to define a polyhedrally constrained
+% variational inequality.  These are transformed internally to a standard
+% mixed complementarity problem.  The polyhedral constraints are of the
+% form
+%       Ax ? b
+% where ? can be <=, =, or >= depending on the type specified for each
+% constraint.
+
+Big = 1e20;
+
+if (nargin < 2)
+  error('two input arguments required for lcp(M, q)');
+end
+
+if (~issparse(M))
+  M = sparse(M);	% Make sure M is sparse
+end
+q = full(q(:)); 	% Make sure q is a column vector
+
+[mm,mn] = size(M);	% Get the size of the inputs
+n = length(q);
+
+if (mm ~= mn | mm ~= n) 
+  error('dimensions of M and q must match');
+end
+
+if (n == 0)
+  error('empty model');
+end
+
+if (nargin < 3 | isempty(l))
+  l = zeros(n,1);
+end
+
+if (nargin < 4 | isempty(u))
+  u = Big*ones(n,1);
+end
+
+if (nargin < 5 | isempty(z))
+  z = zeros(n,1);
+end
+
+z = full(z(:)); l = full(l(:)); u = full(u(:));
+if (length(z) ~= n | length(l) ~= n | length(u) ~= n)
+  error('Input arguments are of incompatible sizes');
+end
+
+l = max(l,-Big*ones(n,1));
+u = min(u,Big*ones(n,1));
+z = min(max(z,l),u);
+
+m = 0;
+if (nargin > 5)
+  if (nargin < 7)
+    error('Polyhedral constraints require A and b');
+  end
+
+  if (~issparse(A))
+    A = sparse(A);
+  end
+  b = full(b(:));
+
+  m = length(b);
+
+  if (m > 0)
+
+    [am, an] = size(A);
+
+    if (am ~= m | an ~= n)
+      error('Polyhedral constraints of incompatible sizes');
+    end
+
+    if (nargin < 8 | isempty(t))
+      t = ones(m,1);
+    end
+
+    if (nargin < 9 | isempty(mu))
+      mu = zeros(m,1);
+    end
+
+    t = full(t(:)); mu = full(mu(:));
+    if (length(t) ~= m | length(mu) ~= m)
+      error('Polyhedral input arguments are of incompatible sizes');
+    end
+
+    l_p = -Big*ones(m,1);
+    u_p =  Big*ones(m,1);
+
+    idx = find(t > 0);
+    if (length(idx) > 0)
+      l_p(idx) = zeros(length(idx),1);
+    end
+
+    idx = find(t < 0);
+    if (length(idx) > 0)
+      u_p(idx) = zeros(length(idx),1);
+    end
+
+    mu = min(max(mu,l_p),u_p);
+
+    M = [M -A'; A sparse(m,m)];
+    q = [q; -b];
+
+    z = [z; mu];
+    l = [l; l_p];
+    u = [u; u_p];
+  else
+    if (nargin >= 9 & ~isempty(mu))
+      error('No polyhedral constraints -- multipliers set.');
+    end
+
+    if (nargin >= 8 & ~isempty(t))
+      error('No polyhedral constraints -- equation types set.');
+    end
+  end
+end
+
+idx = find(l > u);
+if length(idx) > 0
+  error('Bounds infeasible.');
+end
+
+nnzJ = nnz(M);
+
+[status, ttime] = lcppath(n+m, nnzJ, z, l, u, M, q);
+
+%if (status ~= 1) 
+%  status,
+%  error('Path fails to solve problem');
+%end
+
+mu = [];
+if (m > 0)
+  mu = z(n+1:n+m);
+  z = z(1:n);
+end
+
+return;
+
diff --git a/matlab/pathmcp.m b/matlab/pathmcp.m
new file mode 100644
index 000000000..1db4dccb1
--- /dev/null
+++ b/matlab/pathmcp.m
@@ -0,0 +1,197 @@
+function [z,f,J,mu,status] = pathmcp(z,l,u,cpfj,A,b,t,mu)
+% pathmcp(z,l,u,cpfj,A,b,t,mu)
+%
+% Solve a polyhedrally constrained variational inequality using PATH
+%
+% Calling syntax: [z,f,J] = pathmcp(z,l,u,cpfunjac,A,b,t,mu)
+%
+% Input:
+%  z - starting point
+%  l - lower bounds on z
+%  u - upper bounds on z
+%
+%  cpfunjac - 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
+%             'mcp_funjac.m' unless stated otherwise in the variable cpfunjac).
+%
+%             'mcp_funjac.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.
+%
+%  A - constraint matrix
+%  b - right hand side of the constraints
+%  t - types of the constraints
+%      <0 : less than or equal
+%      =0 : equal to
+%      >0 : greater than or equal
+%
+%      We have Ax ? b, ? is the type of constraint
+%
+% Output:
+%  z      - solution
+%  mu     - multipliers on the constraints
+%  f      - function evaluation at the solution
+%  J      - jacobian evaluation at the solution
+
+
+Big = 1e20;
+
+if (nargin < 1)
+  error('one input arguments required for mcp(z)');
+end
+
+z = full(z(:)); 
+n = length(z);
+
+if (n == 0)
+  error('empty model');
+end
+
+if (nargin < 2 | isempty(l))
+  l = zeros(n,1);
+end
+
+if (nargin < 3 | isempty(u))
+  u = Big*ones(n,1);
+end
+
+l = full(l(:)); u = full(u(:));
+if (length(l) ~= n | length(u) ~= n)
+  error('Input arguments are of incompatible sizes');
+end
+
+l = max(l,-Big*ones(n,1));
+u = min(u,Big*ones(n,1));
+z = min(max(z,l),u);
+
+if (nargin < 4 | isempty(cpfj))
+  cpfj = 'mcp_funjac';
+end
+
+m   = 0;
+mu  = [];
+l_p = [];
+u_p = [];
+
+if (nargin > 4)
+  if (nargin < 6)
+    error('Polyhedral constraints require A and b');
+  end
+
+  if (~issparse(A))
+    A = sparse(A);
+  end
+  b = full(b(:));
+
+  m = length(b);
+
+  if (m > 0)
+
+    [am, an] = size(A);
+
+    if (am ~= m | an ~= n)
+      error('Polyhedral constraints of incompatible sizes');
+    end
+
+    if (nargin < 7 | isempty(t))
+      t = ones(m,1);
+    end
+
+    if (nargin < 8 | isempty(mu))
+      mu = zeros(m,1);
+    end
+
+    t = full(t(:)); mu = full(mu(:));
+    if (length(t) ~= m | length(mu) ~= m)
+      error('Polyhedral input arguments are of incompatible sizes');
+    end
+
+    l_p = -Big*ones(m,1);
+    u_p =  Big*ones(m,1);
+
+    idx = find(t > 0);
+    if (length(idx) > 0)
+      l_p(idx) = zeros(length(idx),1);
+    end
+
+    idx = find(t < 0);
+    if (length(idx) > 0)
+      u_p(idx) = zeros(length(idx),1);
+    end
+
+    mu = min(max(mu,l_p),u_p);
+  else
+    if (nargin >= 8 & ~isempty(mu))
+      error('No polyhedral constraints -- multipliers set.');
+    end
+
+    if (nargin >= 7 & ~isempty(t))
+      error('No polyhedral constraints -- equation types set.');
+    end
+  end
+else
+  A = [];
+end
+
+% 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 (domerr > 0)
+  [f,J,domerr] = feval(cpfj,z,1);
+end 
+
+if (domerr > 0)
+  error([cpfj ' not defined at starting point']);
+end
+
+if ~issparse(J)
+  error([cpfj ' must return a sparse Jacobian']);
+end
+
+nnzJ = nzmax(J);
+
+row = n + m;
+ele = nnzJ + 2*nzmax(A);
+
+init = [z; mu];
+low  = [l; l_p];
+upp  = [u; u_p];
+
+if m > 0
+  global mcp_vifunc;
+  global mcp_viconn;
+  global mcp_viconm; 
+  global mcp_viconA; 
+  global mcp_viconb;
+
+  mcp_vifunc = cpfj;
+  mcp_viconn = n;
+  mcp_viconm = m;
+  mcp_viconA = A;
+  mcp_viconb = b;
+
+  [status, ttime, f, J] = mcppath(row, ele, init, low, upp, 'mcp_vifunjac');
+else
+  [status, ttime, f, J] = mcppath(row, ele, init, low, upp, cpfj);
+end
+
+%if (status ~= 1) 
+%  status,
+%  error('Path fails to solve problem');
+%end
+
+mu = [];
+z = init;
+
+if m > 0
+  mu = init(n+1:n+m);
+  z = init(1:n);
+
+  J = J(1:n,1:n);
+  f = f(1:n) + A'*mu;
+end
+
+return;