diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 3033f1143..1b68c688b 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -164,6 +164,15 @@ elseif options_.solve_algo == 3
else
[x,info] = csolve(func,x,[],1e-6,500,varargin{:});
end
+elseif options_.solve_algo == 10
+ olmmcp = options_.lmmcp;
+ olmmcp.Display = 'iter';
+ [x,fval,exitflag] = lmmcp(func,x,olmmcp.lb,olmmcp.ub,olmmcp,varargin{:});
+ if exitflag == 1
+ info = 0;
+ else
+ info = 1;
+ end
else
error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9]')
end
diff --git a/matlab/lmmcp/catstruct.m b/matlab/lmmcp/catstruct.m
new file mode 100644
index 000000000..8248366e4
--- /dev/null
+++ b/matlab/lmmcp/catstruct.m
@@ -0,0 +1,148 @@
+function A = catstruct(varargin)
+% CATSTRUCT Concatenate or merge structures with different fieldnames
+% X = CATSTRUCT(S1,S2,S3,...) merges the structures S1, S2, S3 ...
+% into one new structure X. X contains all fields present in the various
+% structures. An example:
+%
+% A.name = 'Me' ;
+% B.income = 99999 ;
+% X = catstruct(A,B)
+% % -> X.name = 'Me' ;
+% % X.income = 99999 ;
+%
+% If a fieldname is not unique among structures (i.e., a fieldname is
+% present in more than one structure), only the value from the last
+% structure with this field is used. In this case, the fields are
+% alphabetically sorted. A warning is issued as well. An axample:
+%
+% S1.name = 'Me' ;
+% S2.age = 20 ; S3.age = 30 ; S4.age = 40 ;
+% S5.honest = false ;
+% Y = catstruct(S1,S2,S3,S4,S5) % use value from S4
+%
+% The inputs can be array of structures. All structures should have the
+% same size. An example:
+%
+% C(1).bb = 1 ; C(2).bb = 2 ;
+% D(1).aa = 3 ; D(2).aa = 4 ;
+% CD = catstruct(C,D) % CD is a 1x2 structure array with fields bb and aa
+%
+% The last input can be the string 'sorted'. In this case,
+% CATSTRUCT(S1,S2, ..., 'sorted') will sort the fieldnames alphabetically.
+% To sort the fieldnames of a structure A, you could use
+% CATSTRUCT(A,'sorted') but I recommend ORDERFIELDS for doing that.
+%
+% When there is nothing to concatenate, the result will be an empty
+% struct (0x0 struct array with no fields).
+%
+% NOTE: To concatenate similar arrays of structs, you can use simple
+% concatenation:
+% A = dir('*.mat') ; B = dir('*.m') ; C = [A ; B] ;
+%
+% See also CAT, STRUCT, FIELDNAMES, STRUCT2CELL, ORDERFIELDS
+
+% for Matlab R13 and up
+% version 3.0 (mar 2013)
+% (c) Jos van der Geest
+% email: jos@jasen.nl
+
+% (C) 2013 Christophe Gouel
+
+% History
+% Created in 2005
+% Revisions
+% 2.0 (sep 2007) removed bug when dealing with fields containing cell
+% arrays (Thanks to Rene Willemink)
+% 2.1 (sep 2008) added warning and error identifiers
+% 2.2 (oct 2008) fixed error when dealing with empty structs (Thanks to
+% Lars Barring)
+% 3.0 (mar 2013) fixed problem when the inputs were array of structures
+% (thanks to Tor Inge Birkenes for pointing this out).
+% Rephrased the help section as well.
+
+error(nargchk(1,Inf,nargin)) ;
+N = nargin ;
+
+if ~isstruct(varargin{end}),
+ if isequal(varargin{end},'sorted'),
+ sorted = 1 ;
+ N = N-1 ;
+ error(nargchk(1,Inf,N)) ;
+ else
+ error('catstruct:InvalidArgument','Last argument should be a structure, or the string "sorted".') ;
+ end
+else
+ sorted = 0 ;
+end
+
+sz0 = [] ; % used to check that all inputs have the same size
+
+% used to check for a few trivial cases
+NonEmptyInputs = false(N,1) ;
+NonEmptyInputsN = 0 ;
+
+% used to collect the fieldnames and the inputs
+FN = cell(N,1) ;
+VAL = cell(N,1) ;
+
+% parse the inputs
+for ii=1:N,
+ X = varargin{ii} ;
+ if ~isstruct(X),
+ error('catstruct:InvalidArgument',['Argument #' num2str(ii) ' is not a structure.']) ;
+ end
+
+ if ~isempty(X),
+ % empty structs are ignored
+ if ii > 1 && ~isempty(sz0)
+ if ~isequal(size(X), sz0)
+ error('catstruct:UnequalSizes','All structures should have the same size.') ;
+ end
+ else
+ sz0 = size(X) ;
+ end
+ NonEmptyInputsN = NonEmptyInputsN + 1 ;
+ NonEmptyInputs(ii) = true ;
+ FN{ii} = fieldnames(X) ;
+ VAL{ii} = struct2cell(X) ;
+ end
+end
+
+if NonEmptyInputsN == 0
+ % all structures were empty
+ A = struct([]) ;
+elseif NonEmptyInputsN == 1,
+ % there was only one non-empty structure
+ A = varargin{NonEmptyInputs} ;
+ if sorted,
+ A = orderfields(A) ;
+ end
+else
+ % there is actually something to concatenate
+ FN = cat(1,FN{:}) ;
+ VAL = cat(1,VAL{:}) ;
+ FN = squeeze(FN) ;
+ VAL = squeeze(VAL) ;
+ MatlabVersion = version;
+ if str2double(MatlabVersion(end-5:end-2))<2013 % Equivalent to, but faster than if verLessThan('matlab','8.1')
+ [UFN,ind] = unique(FN) ;
+ else
+ [UFN,ind] = unique(FN,'legacy') ;
+ end
+
+ if numel(UFN) ~= numel(FN),
+ warning('catstruct:DuplicatesFound','Fieldnames are not unique between structures.') ;
+ sorted = 1 ;
+ end
+
+ if sorted,
+ VAL = VAL(ind,:) ;
+ FN = FN(ind,:) ;
+ end
+
+ A = cell2struct(VAL, FN);
+ A = reshape(A, sz0) ; % reshape into original format
+end
+
+
+
diff --git a/matlab/lmmcp/dyn_lmmcp.m b/matlab/lmmcp/dyn_lmmcp.m
new file mode 100644
index 000000000..ce3db1c7b
--- /dev/null
+++ b/matlab/lmmcp/dyn_lmmcp.m
@@ -0,0 +1,77 @@
+function [endo_simul,info] = dyn_lmmcp(M,options,oo)
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+[lb,ub,eq_index] = get_complementarity_conditions(M);
+
+lead_lag_incidence = M.lead_lag_incidence;
+
+ny = M.endo_nbr;
+
+max_lag = M.maximum_endo_lag;
+
+nyp = nnz(lead_lag_incidence(1,:)) ;
+iyp = find(lead_lag_incidence(1,:)>0) ;
+ny0 = nnz(lead_lag_incidence(2,:)) ;
+iy0 = find(lead_lag_incidence(2,:)>0) ;
+nyf = nnz(lead_lag_incidence(3,:)) ;
+iyf = find(lead_lag_incidence(3,:)>0) ;
+
+nd = nyp+ny0+nyf;
+nrc = nyf+1 ;
+isp = [1:nyp] ;
+is = [nyp+1:ny+nyp] ;
+isf = iyf+nyp ;
+isf1 = [nyp+ny+1:nyf+nyp+ny+1] ;
+stop = 0 ;
+iz = [1:ny+nyp+nyf];
+
+periods = options.periods;
+steady_state = oo.steady_state;
+params = M.params;
+endo_simul = oo.endo_simul;
+exo_simul = oo.exo_simul;
+i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)');
+i_cols_A1 = find(lead_lag_incidence(2:3,:)');
+i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
+i_cols_j = 1:nd;
+i_upd = ny+(1:periods*ny);
+
+x = endo_simul(:);
+
+model_dynamic = str2func([M.fname,'_dynamic']);
+z = x(find(lead_lag_incidence'));
+[res,A] = model_dynamic(z, exo_simul, params, steady_state,2);
+nnzA = nnz(A);
+
+LB = repmat(lb,periods,1);
+UB = repmat(ub,periods,1);
+
+Y0 = endo_simul(:,1);
+YT = endo_simul(:,end);
+x = endo_simul(:,2:end-1);
+x = x(:);
+
+func_handle = @(x) dyn_lmmcp_func(x,model_dynamic, Y0, YT, exo_simul, ...
+ params, steady_state, periods, ny, ...
+ lead_lag_incidence, i_cols_A1, i_cols_1, ...
+ i_cols_T, i_cols_j,nnzA,eq_index);
+
+[x, info] = lmmcp(func_handle,x,LB,UB);
+
+endo_simul = [Y0 reshape(x,ny,periods) YT];
\ No newline at end of file
diff --git a/matlab/lmmcp/dyn_lmmcp_func.m b/matlab/lmmcp/dyn_lmmcp_func.m
new file mode 100644
index 000000000..936f98306
--- /dev/null
+++ b/matlab/lmmcp/dyn_lmmcp_func.m
@@ -0,0 +1,56 @@
+function [F,A] = dyn_lmmcp_func(x, model_dynamic, Y0, YT, exo_simul, params, ...
+ steady_state, periods, ny, lead_lag_incidence, ...
+ i_cols_A1, i_cols_1, i_cols_T, i_cols_j, ...
+ nnzA,eq_index)
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+ Y = [Y0; x; YT];
+
+ F = zeros(periods*ny,1);
+ if nargout == 2
+ A = sparse([],[],[],periods*ny,periods*ny,periods*nnzA);
+ end
+
+ i_rows = 1:ny;
+ i_cols = find(lead_lag_incidence');
+ i_cols_A = i_cols;
+
+ for it = 2:(periods+1)
+
+ [res,jacobian] = model_dynamic(Y(i_cols),exo_simul, params, ...
+ steady_state,it);
+ F(i_rows) = res(eq_index);
+
+ if nargout == 2
+ if it == 2
+ A(i_rows,i_cols_A1) = jacobian(eq_index,i_cols_1);
+ elseif it == periods+1
+ A(i_rows,i_cols_A(i_cols_T)) = jacobian(eq_index,i_cols_T);
+ else
+ A(i_rows,i_cols_A) = jacobian(eq_index,i_cols_j);
+ end
+ end
+
+ i_rows = i_rows + ny;
+ i_cols = i_cols + ny;
+ if nargout == 2 && it > 2
+ i_cols_A = i_cols_A + ny;
+ end
+ end
+
diff --git a/matlab/lmmcp/get_complementarity_conditions.m b/matlab/lmmcp/get_complementarity_conditions.m
new file mode 100644
index 000000000..a58f8e337
--- /dev/null
+++ b/matlab/lmmcp/get_complementarity_conditions.m
@@ -0,0 +1,53 @@
+function [lb,ub,eq_index] = get_complementarity_condition(M)
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+
+etags = M.equations_tags;
+ub = inf(M.endo_nbr,1);
+lb = -ub;
+eq_index = (1:M.endo_nbr)';
+for i=1:size(etags,1)
+ if strcmp(etags{i,2},'mcp')
+ [b,m] = strsplit(etags{i,3},{' ','<','>'});
+ if length(m) == 1
+ if any(m{1} == '<')
+ k = find(strcmp(b{1},cellstr(M.endo_names)));
+ if isempty(k)
+ error(sprintf(['Complementarity condition %s: variable %s is ' ...
+ 'not recognized',etags{i,3},b{1}]))
+ end
+ ub(k) = str2num(b{2});
+ eq_index(etags{i,1}) = k;
+ eq_index(k) = etags{i,1};
+ elseif any(m{1} == '>')
+ k = find(strcmp(b{1},cellstr(M.endo_names)));
+ if isempty(k)
+ error(sprintf(['Complementarity condition %s: variable %s is ' ...
+ 'not recognized',etags{i},b{1}]))
+ end
+ lb(k) = str2num(b{2});
+ eq_index(etags{i,1}) = k;
+ eq_index(k) = etags{i,1};
+ end
+ else
+ error(sprintf('Complementarity condition %s can''t be parsed',etags{i,3}))
+ end
+ end
+end
+
diff --git a/matlab/lmmcp/lmmcp.m b/matlab/lmmcp/lmmcp.m
new file mode 100644
index 000000000..1122b83ed
--- /dev/null
+++ b/matlab/lmmcp/lmmcp.m
@@ -0,0 +1,615 @@
+function [x,FVAL,EXITFLAG,OUTPUT,JACOB] = lmmcp(FUN,x,lb,ub,options,varargin)
+% LMMCP solves a mixed complementarity problem.
+%
+% LMMCP uses a semismooth least squares formulation. The method applies a
+% Levenberg-Marquardt/Gauss-Newton algorithm to a least-squares formulation.
+%
+% X = LMMCP(FUN,X0) tries to solve the system of nonlinear equations F(X)=0 and
+% starts at the vector X0. FUN accepts a vector X and return a vector F of equation
+% values F evaluated at X and, as second output if required, a matrix J, the
+% Jacobian evaluated at X.
+%
+% X = LMMCP(FUN,X0,LB,UB) solves the mixed complementarity problem of the form:
+% LB =X => F(X)>0,
+% LB<=X<=UB => F(X)=0,
+% X =UB => F(X)<0.
+%
+% X = LMMCP(FUN,X0,LB,UB,OPTIONS) solves the MCP problem using the options
+% defined in the structure OPTIONS. Main fields are
+% Display : control the display of iterations, 'none' (default),
+% 'iter-detailed' or 'final-detailed'
+% Switch from phase I to phase II
+% preprocess : activate preprocessor for phase I (default = 1)
+% presteps : number of iterations in phase I (default = 20)
+% Termination parameters
+% epsilon2 : termination value of the merit function (default = 1E-16)
+% MaxIter : maximum number of iterations (default = 500)
+% tmin : safeguard stepsize (default = 1E-12)
+% Stepsize parameters
+% m : number of previous function values to use in the nonmonotone
+% line search rule (default = 10)
+% kwatch : maximum number of steps (default = 20 and should not be
+% smaller than m)
+% watchdog : activate the watchdog strategy (default = 1)
+% Ther are other minor parameters. Please see the code for their default values
+% and interpretation.
+%
+% [X,FVAL] = LMMCP(FUN,X0,...) returns the value of the equations FUN at X.
+%
+% [X,FVAL,EXITFLAG] = LMMCP(FUN,X0,...) returns EXITFLAG that describes the exit
+% conditions. Possible values are
+% 1 : LMMCP converged to a root
+% 0 : Too many iterations
+% -1 :
+%
+% [X,FVAL,EXITFLAG,OUTPUT] = LMMCP(FUN,X0,...) returns the structure OUTPUT that
+% contains the number of iterations (OUTPUT.iterations), the value of the merit
+% function (OUTPUT.Psix), and the norm of the derivative of the merit function
+% (OUTPUT.normDPsix).
+%
+% [X,FVAL,EXITFLAG,OUTPUT,JACOB] = LMMCP(FUN,X0,...) returns JACOB the Jacobian
+% of FUN evaluated at X.
+%
+% More details of the main program may be found in the following paper:
+%
+% Christian Kanzow and Stefania Petra: On a semismooth least squares formulation of
+% complementarity problems with gap reduction. Optimization Methods and Software
+% 19, 2004, pp. 507-525.
+%
+% In addition, the current implementation uses a preprocessor which is the
+% projected Levenberg-Marquardt step from the following preprint:
+%
+% Christian Kanzow and Stefania Petra: Projected filter trust region methods for a
+% semismooth least squares formulation of mixed complementarity
+% problems. Optimization Methods and Software
+% 22, 2007, pp. 713-735.
+%
+% A user's guide is also available:
+%
+% Christian Kanzow and Stefania Petra (2005).
+% LMMCP --- A Levenberg-Marquardt-type MATLAB Solver for Mixed Complementarity Problems.
+% University of Wuerzburg.
+% http://www.mathematik.uni-wuerzburg.de/~kanzow/software/UserGuide.pdf
+%
+% This is a modification by Christophe Gouel of the original files, which can be
+% downloaded from:
+% http://www.mathematik.uni-wuerzburg.de/~kanzow/software/LMMCP.zip
+%
+% copyright: Christian Kanzow and Stefania Petra
+% Institute of Applied Mathematics and Statistics
+% University of Wuerzburg
+% Am Hubland
+% 97074 Wuerzburg
+% GERMANY
+%
+% e-mail: kanzow@mathematik.uni-wuerzburg.de
+% petra@mathematik.uni-wuerzburg.de
+
+% ------------Initialization----------------
+defaultopt = struct(...
+ 'beta', 0.55,...
+ 'Big', 1e10,...
+ 'delta', 5,...
+ 'deltamin', 1,...
+ 'Display', 'none',...
+ 'epsilon1', 1e-6,...
+ 'epsilon2', 1e-16,...
+ 'eta', 0.95,...
+ 'kwatch', 20,...
+ 'lambda1', 0.1,...
+ 'm', 10,...
+ 'MaxIter', 500,...
+ 'null', 1e-10,...
+ 'preprocess', 1,...
+ 'presteps', 20,...
+ 'sigma', 1e-4,...
+ 'sigma1', 0.5,...
+ 'sigma2', 2,...
+ 'tmin', 1e-12,...
+ 'watchdog', 1);
+
+if nargin < 4
+ ub = inf(size(x));
+ if nargin < 3
+ lb = -inf(size(x));
+ end
+end
+
+if nargin < 5 || isempty(options) || ~isstruct(options)
+ options = defaultopt;
+else
+ warning('off','catstruct:DuplicatesFound')
+ options = catstruct(defaultopt,options);
+end
+
+
+switch options.Display
+ case {'off','none'}
+ verbosity = 0;
+ case {'iter','iter-detailed'}
+ verbosity = 2;
+ case {'final','final-detailed'}
+ verbosity = 1;
+ otherwise
+ verbosity = 0;
+end
+
+% parameter settings
+eps1 = options.epsilon1;
+eps2 = options.epsilon2;
+null = options.null;
+Big = options.Big;
+
+% maximal number of iterations
+kmax = options.MaxIter;
+
+% choice of lambda
+lambda1 = options.lambda1;
+lambda2 = 1-lambda1;
+
+% steplength parameters
+beta = options.beta;
+sigma = options.sigma;
+tmin = options.tmin;
+
+% parameters watchdog and nonmonotone line search; redefined later
+m = options.m;
+kwatch = options.kwatch;
+watchdog = options.watchdog; % 1=watchdog strategy active, otherwise not
+
+% parameters for preprocessor
+preprocess = options.preprocess; % 1=preprocessor used, otherwise not
+presteps = options.presteps; % maximum number of preprocessing steps
+
+% trust-region parameters for preprocessor
+delta = options.delta;
+deltamin = options.deltamin;
+sigma1 = options.sigma1;
+sigma2 = options.sigma2;
+eta = options.eta;
+
+% initializations
+k = 0;
+k_main = 0;
+
+% compute a feasible starting point by projection
+x = max(lb,min(x,ub));
+
+n = length(x);
+OUTPUT.Dim = n;
+
+% definition of index sets I_l, I_u and I_lu
+Indexset = zeros(n,1);
+I_l = lb>-Big & ub>Big;
+I_u = lb<-Big & ub-Big & ub 1
+ fprintf(' k Psi(x) || DPsi(x) || stepsize\n');
+ disp('====================================================================')
+ disp('********************* Output at starting point *********************')
+ fprintf('%4.0f %24.5e %24.5e\n',k,Psix,normDPsix);
+end
+
+%
+% Preprocessor using local method
+%
+
+if preprocess==1
+
+ if verbosity > 1
+ disp('************************** Preprocessor ****************************')
+ end
+
+ normpLM=1;
+ while (k < presteps) && (Psix > eps2) && (normpLM>null)
+ k=k+1;
+
+ % choice of Levenberg-Marquardt parameter, note that we do not use
+ % the condition estimator for large-scale problems, although this
+ % may cause numerical problems in some examples
+
+ i=0;
+ mu=0;
+ if n<100000
+ i=1;
+ mu=1e-16;
+ if condest(DPhix'*DPhix)>1e25
+ mu=1e-6/(k+1);
+ end
+ end
+ if i==1
+ pLM= [ DPhix ; sqrt(mu)*speye(n)]\[-Phix;sparse(n,1)];
+ else
+ pLM=-DPhix\Phix;
+ end
+ normpLM=norm(pLM);
+
+ % compute the projected Levenberg-Marquard step onto box Xk
+ lbnew=max(min(lb-x,0),-delta);
+ ubnew=min(max(ub-x,0),delta);
+ d=max(lbnew,min(pLM,ubnew));
+ xnew=x+d;
+
+ % function evaluations etc.
+ [Fxnew,DFxnew] = feval(FUN,xnew,varargin{:});
+ Phixnew=Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
+ Psixnew=0.5*(Phixnew'*Phixnew);
+ normPhixnew=norm(Phixnew);
+
+ % update of delta
+ if normPhixnew<=eta*normPhix
+ delta=max(deltamin,sigma2*delta);
+ elseif normPhixnew>5*eta*normPhix
+ delta=max(deltamin,sigma1*delta);
+ end
+
+ % update
+ x=xnew;
+ Fx=Fxnew;
+ DFx=DFxnew;
+ Phix=Phixnew;
+ Psix=Psixnew;
+ normPhix=normPhixnew;
+ DPhix=DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset);
+ DPsix=DPhix'*Phix;
+ normDPsix=norm(DPsix,inf);
+
+ % output at each iteration
+ t=1;
+ if verbosity > 1
+ fprintf('%4.0f %24.5e %24.5e %11.7g\n',k,Psix,normDPsix,t);
+ end
+ end
+end
+
+% terminate program or redefine current iterate as original initial point
+if preprocess==1 && Psix 0
+ fprintf('Psix = %1.4e\nnormDPsix = %1.4e\n',Psix,normDPsix);
+ disp('Approximate solution found.')
+ end
+ EXITFLAG = 1;
+ FVAL = Fx;
+ OUTPUT.iterations = k;
+ OUTPUT.Psix = Psix;
+ OUTPUT.normDPsix = normDPsix;
+ JACOB = DFx;
+ return
+elseif preprocess==1 && Psix>=eps2
+ x=x0;
+ Phix=Phix0;
+ Psix=Psix0;
+ DPhix=DPhix0;
+ DPsix=DPsix0;
+ if verbosity > 1
+ disp('******************** Restart with initial point ********************')
+ fprintf('%4.0f %24.5e %24.5e\n',k_main,Psix0,normDPsix0);
+ end
+end
+
+%
+% Main algorithm
+%
+
+if verbosity > 1
+ disp('************************** Main program ****************************')
+end
+
+while (k < kmax) && (Psix > eps2)
+
+ % choice of Levenberg-Marquardt parameter, note that we do not use
+ % the condition estimator for large-scale problems, although this
+ % may cause numerical problems in some examples
+
+ i=0;
+ if n<100
+ i=1;
+ mu=1e-16;
+ if condest(DPhix'*DPhix)>1e25
+ mu=1e-1/(k+1);
+ end
+ end
+
+ % compute a Levenberg-Marquard direction
+
+ if i==1
+ d= [ DPhix ; sqrt(mu)*eye(n)]\[-Phix;zeros(n,1)];
+ else
+ d=-DPhix\Phix;
+ end
+
+ % computation of steplength t using the nonmonotone Armijo-rule
+ % starting with the 6-th iteration
+
+ % computation of steplength t using the monotone Armijo-rule if
+ % d is a 'good' descent direction or k<=5
+
+ t = 1;
+ xnew = x+d;
+ Fxnew = feval(FUN,xnew,varargin{:});
+ Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
+ Psixnew = 0.5*(Phixnew'*Phixnew);
+ const = sigma*DPsix'*d;
+
+ while (Psixnew > MaxPsi + const*t) && (t > tmin)
+ t = t*beta;
+ xnew = x+t*d;
+ Fxnew = feval(FUN,xnew,varargin{:});
+ Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
+ Psixnew = 0.5*(Phixnew'*Phixnew);
+ end
+
+ % updatings
+ x=xnew;
+ Fx=Fxnew;
+ Phix=Phixnew;
+ Psix=Psixnew;
+ [~,DFx]=feval(FUN,x,varargin{:});
+ DPhix=DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset);
+ DPsix=DPhix'*Phix;
+ normDPsix=norm(DPsix);
+ k=k+1;
+ k_main=k_main+1;
+
+ if k_main<=5
+ aux(mod(k_main,m)+1)=Psix;
+ MaxPsi=Psix;
+ else
+ aux(mod(k_main,m)+1) = Psix;
+ MaxPsi = max(aux);
+ end;
+
+ % updatings for the watchdog strategy
+ if watchdog ==1
+ if Psixkwatch
+ x=xbest;
+ Phix=Phibest;
+ Psix=Psibest;
+ DPhix=DPhibest;
+ DPsix=DPsibest;
+ normDPsix=normDPsibest;
+ MaxPsi=Psix;
+ end;
+ end;
+
+ if verbosity > 1
+ % output at each iteration
+ fprintf('%4.0f %24.5e %24.5e %11.7g\n',k,Psix,normDPsix,t);
+ end
+end;
+
+% final output
+if Psix<=eps2
+ EXITFLAG = 1;
+ if verbosity > 0, disp('Approximate solution found.'); end
+elseif k>=kmax
+ EXITFLAG = 0;
+ if verbosity > 0, disp('Maximum iteration number reached.'); end
+elseif normDPsix<=eps1
+ EXITFLAG = -1; % Provisoire
+ if verbosity > 0, disp('Approximate stationary point found.'); end
+else
+ EXITFLAG = -1; % Provisoire
+ if verbosity > 0, disp('No solution found.'); end
+end
+
+FVAL = Fx;
+OUTPUT.iterations = k;
+OUTPUT.Psix = Psix;
+OUTPUT.normDPsix = normDPsix;
+JACOB = DFx;
+
+% Subfunctions
+
+function y = Phi(x,Fx,lb,ub,lambda1,lambda2,n,Indexset)
+
+y = zeros(2*n,1);
+phi_u = sqrt((ub-x).^2+Fx.^2)-ub+x+Fx;
+LZ = false(n,1); % logical zero
+
+I0 = Indexset==0;
+y(I0) = -lambda1*Fx(I0);
+y([LZ; I0]) = -lambda2*Fx(I0);
+
+I1 = Indexset==1;
+y(I1) = lambda1*(-x(I1)+lb(I1)-Fx(I1)+sqrt((x(I1)-lb(I1)).^2+Fx(I1).^2));
+y([LZ; I1]) = lambda2*max(0,x(I1)-lb(I1)).*max(0,Fx(I1));
+
+I2 = Indexset==2;
+y(I2) = -lambda1*phi_u(I2);
+y([LZ; I2]) = lambda2*max(0,ub(I2)-x(I2)).*max(0,-Fx(I2));
+
+I3 = Indexset==3;
+y(I3) = lambda1*(sqrt((x(I3)-lb(I3)).^2+phi_u(I3).^2)-x(I3)+lb(I3)-phi_u(I3));
+y([LZ; I3]) = lambda2*(max(0,x(I3)-lb(I3)).*max(0,Fx(I3))+max(0,ub(I3)-x(I3)).*max(0,-Fx(I3)));
+
+
+function H = DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset)
+% DPHI evaluates an element of the C-subdifferential of operator Phi
+
+null = 1e-8;
+beta_l = zeros(n,1);
+beta_u = zeros(n,1);
+alpha_l = zeros(n,1);
+alpha_u = zeros(n,1);
+
+
+z = zeros(n,1);
+H2 = sparse(n,n);
+
+I = abs(x-lb)<=null & abs(Fx)<=null;
+beta_l(I) = 1;
+z(I) = 1;
+
+I = abs(ub-x)<=null & abs(Fx)<=null;
+beta_u(I) = 1;
+z(I) = 1;
+
+I = x-lb>=-null & Fx>=-null;
+alpha_l(I) = 1;
+
+I = ub-x>=-null & Fx<=null;
+alpha_u(I) = 1;
+
+Da = zeros(n,1);
+Db = zeros(n,1);
+
+I = 1:n;
+
+I0 = Indexset==0;
+Da(I0) = 0;
+Db(I0) = -1;
+H2(I0,:) = -DFx(I0,:);
+
+I1 = Indexset==1;
+denom1 = zeros(n,1);
+denom2 = zeros(n,1);
+if any(I1)
+ denom1(I1) = max(null,sqrt((x(I1)-lb(I1)).^2+Fx(I1).^2));
+ denom2(I1) = max(null,sqrt(z(I1).^2+(DFx(I1,:)*z).^2));
+end
+
+I1b = Indexset==1 & beta_l==0;
+Da(I1b) = (x(I1b)-lb(I1b))./denom1(I1b)-1;
+Db(I1b) = Fx(I1b)./denom1(I1b)-1;
+I1b = Indexset==1 & beta_l~=0;
+if any(I1b)
+ Da(I1b) = z(I1b)./denom2(I1b)-1;
+ Db(I1b) = (DFx(I1b,:)*z)./denom2(I1b)-1;
+end
+
+I1a = I(Indexset==1 & alpha_l==1);
+if any(I1a)
+ H2(I1a,:) = repmat(x(I1a)-lb(I1a),1,n).*DFx(I1a,:)+sparse(1:length(I1a),I1a, ...
+ Fx(I1a),length(I1a),n,length(I1a));
+end
+
+I2 = Indexset==2;
+denom1 = zeros(n,1);
+denom2 = zeros(n,1);
+if any(I2)
+ denom1(I2) = max(null,sqrt((ub(I2)-x(I2)).^2+Fx(I2).^2));
+ denom2(I2) = max(null,sqrt(z(I2).^2+(DFx(I2,:)*z).^2));
+end
+
+I2b = Indexset==2 & beta_u==0;
+Da(I2b) = (ub(I2b)-x(I2b))./denom1(I2b)-1;
+Db(I2b) = -Fx(I2b)./denom1(I2b)-1;
+I2b = Indexset==2 & beta_u~=0;
+if any(I2b)
+ Da(I2b) = -z(I2b)./denom2(I2b)-1;
+ Db(I2b) = -(DFx(I2b,:)*z)./denom2(I2b)-1;
+end
+
+I2a = I(Indexset==2 & alpha_u==1);
+if any(I2a)
+ H2(I2a,:) = repmat(x(I2a)-ub(I2a),1,n).*DFx(I2a,:)+sparse(1:length(I2a),I2a, ...
+ Fx(I2a),length(I2a),n,length(I2a));
+end
+
+I3 = Indexset==3;
+phi = zeros(n,1);
+ai = zeros(n,1);
+bi = zeros(n,1);
+ci = zeros(n,1);
+di = zeros(n,1);
+denom1 = zeros(n,1);
+denom2 = zeros(n,1);
+denom3 = zeros(n,1);
+denom4 = zeros(n,1);
+if any(I3)
+ phi(I3) = -ub(I3)+x(I3)+Fx(I3)+sqrt((ub(I3)-x(I3)).^2+Fx(I3).^2);
+ denom1(I3) = max(null,sqrt((x(I3)-lb(I3)).^2+phi(I3).^2));
+ denom2(I3) = max(null,sqrt(z(I3).^2+(DFx(I3,:)*z).^2));
+ denom3(I3) = max(null,sqrt((ub(I3)-x(I3)).^2+Fx(I3).^2));
+ denom4(I3) = max(null,sqrt(z(I3).^2));
+end
+
+I3bu = Indexset==3 & beta_u==0;
+ci(I3bu) = (x(I3bu)-ub(I3bu))./denom3(I3bu)+1;
+di(I3bu) = Fx(I3bu)./denom3(I3bu)+1;
+I3bu = Indexset==3 & beta_u~=0;
+if any(I3bu)
+ ci(I3bu) = 1+z(I3bu)./denom2(I3bu);
+ di(I3bu) = 1+(DFx(I3bu,:)*z)./denom2(I3bu);
+end
+
+I3bl = Indexset==3 & beta_l==0;
+ai(I3bl) = (x(I3bl)-lb(I3bl))./denom1(I3bl)-1;
+bi(I3bl) = phi(I3bl)./denom1(I3bl)-1;
+I3bl = Indexset==3 & beta_l~=0;
+if any(I3bl)
+ ai(I3bl) = z(I3bl)./denom4(I3bl)-1;
+ bi(I3bl) = (ci(I3bl).*z(I3bl)+(di(I3bl,ones(1,n)).*DFx(I3bl,:))*z)./denom4(I3bl)-1;
+end
+
+Da(I3) = ai(I3)+bi(I3).*ci(I3);
+Db(I3) = bi(I3).*di(I3);
+
+I3a = I(Indexset==3 & alpha_l==1 & alpha_u==1);
+if any(I3a)
+ H2(I3a,:) = repmat(-lb(I3a)-ub(I3a)+2*x(I3a),1,n).*DFx(I3a,:)+...
+ 2*sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
+end
+I3a = I(Indexset==3 & alpha_l==1 & alpha_u~=1);
+if any(I3a)
+ H2(I3a,:) = repmat(x(I3a)-lb(I3a),1,n).*DFx(I3a,:)+...
+ sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
+end
+I3a = I(Indexset==3 & alpha_l~=1 & alpha_u==1);
+if any(I3a)
+ H2(I3a,:) = repmat(x(I3a)-ub(I3a),1,n).*DFx(I3a,:)+...
+ sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
+end
+
+%H1 = sparse(1:n,1:n,Da,n,n,n)+Db(:,ones(n,1)).*DFx;
+H1 = bsxfun(@times,Db,DFx);
+H1 = spdiags(diag(H1)+Da,0,H1);
+H = [lambda1*H1; lambda2*H2];