2020-02-20 13:07:54 +01:00
|
|
|
|
function [x, errorflag, fvec, fjac] = dynare_solve(f, x, options, varargin)
|
|
|
|
|
|
|
|
|
|
% Solves a nonlinear system of equations, f(x) = 0 with n unknowns
|
|
|
|
|
% and n equations.
|
2008-01-03 13:12:07 +01:00
|
|
|
|
%
|
|
|
|
|
% INPUTS
|
2020-02-20 13:07:54 +01:00
|
|
|
|
% - f [char, fhandle] function to be solved
|
|
|
|
|
% - x [double] n×1 vector, initial guess.
|
|
|
|
|
% - options [struct] Dynare options, aka options_.
|
|
|
|
|
% - varargin list of additional arguments to be passed to func.
|
2014-09-10 16:30:17 +02:00
|
|
|
|
%
|
2008-01-03 13:12:07 +01:00
|
|
|
|
% OUTPUTS
|
2020-02-20 13:07:54 +01:00
|
|
|
|
% - x [double] n×1 vector, solution.
|
|
|
|
|
% - errorflag [logical] scalar, true iff the model can not be solved.
|
|
|
|
|
% - fvec [double] n×1 vector, function value at x (f(x), used for debugging when errorflag is true).
|
|
|
|
|
% - fjac [double] n×n matrix, Jacobian value at x (J(x), used for debugging when errorflag is true).
|
2020-09-25 23:45:21 +02:00
|
|
|
|
|
2021-04-08 13:28:07 +02:00
|
|
|
|
% Copyright © 2001-2021 Dynare Team
|
2008-08-01 14:40:33 +02:00
|
|
|
|
%
|
|
|
|
|
% 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
|
2021-06-09 17:33:48 +02:00
|
|
|
|
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
2008-01-03 13:12:07 +01:00
|
|
|
|
|
2020-02-20 13:07:54 +01:00
|
|
|
|
jacobian_flag = options.jacobian_flag; % true iff Jacobian is returned by f routine (as a second output argument).
|
2009-12-16 18:17:34 +01:00
|
|
|
|
|
2015-02-18 23:52:17 +01:00
|
|
|
|
% Set tolerance parameter depending the the caller function.
|
|
|
|
|
stack = dbstack;
|
2016-11-10 17:58:47 +01:00
|
|
|
|
if isoctave
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[~, name, ext]=fileparts(stack(2).file);
|
2016-11-10 17:58:47 +01:00
|
|
|
|
caller_file_name=[name,ext];
|
|
|
|
|
else
|
2021-06-04 17:58:48 +02:00
|
|
|
|
if size(stack,1)>1
|
|
|
|
|
caller_file_name=stack(2).file;
|
|
|
|
|
else
|
|
|
|
|
caller_file_name=stack(1).file;
|
|
|
|
|
end
|
2016-11-10 17:58:47 +01:00
|
|
|
|
end
|
2021-04-08 20:26:23 +02:00
|
|
|
|
if strcmp(caller_file_name, 'solve_stacked_problem.m') || strcmp(caller_file_name, 'sim1_purely_backward.m')
|
2015-05-25 10:06:13 +02:00
|
|
|
|
tolf = options.dynatol.f;
|
2020-07-06 13:11:23 +02:00
|
|
|
|
tolx = options.dynatol.x;
|
2015-02-18 23:52:17 +01:00
|
|
|
|
else
|
2015-05-25 10:06:13 +02:00
|
|
|
|
tolf = options.solve_tolf;
|
2020-07-06 13:11:23 +02:00
|
|
|
|
tolx = options.solve_tolx;
|
2015-02-18 23:52:17 +01:00
|
|
|
|
end
|
|
|
|
|
|
2016-11-10 17:58:47 +01:00
|
|
|
|
if strcmp(caller_file_name,'dyn_ramsey_static.m')
|
2016-08-02 17:06:51 +02:00
|
|
|
|
maxit = options.ramsey.maxit;
|
|
|
|
|
else
|
|
|
|
|
maxit = options.steady.maxit;
|
|
|
|
|
end
|
|
|
|
|
|
2020-02-20 13:07:54 +01:00
|
|
|
|
errorflag = false;
|
2012-11-01 09:26:05 +01:00
|
|
|
|
nn = size(x,1);
|
2012-10-31 11:08:22 +01:00
|
|
|
|
|
2017-05-10 23:27:41 +02:00
|
|
|
|
% Get status of the initial guess (default values?)
|
|
|
|
|
if any(x)
|
|
|
|
|
% The current initial guess is not the default for all the variables.
|
2020-02-20 13:07:54 +01:00
|
|
|
|
idx = find(x); % Indices of the variables with default initial guess values.
|
2017-05-10 23:27:41 +02:00
|
|
|
|
in0 = length(idx);
|
|
|
|
|
else
|
|
|
|
|
% The current initial guess is the default for all the variables.
|
|
|
|
|
idx = transpose(1:nn);
|
|
|
|
|
in0 = nn;
|
|
|
|
|
end
|
|
|
|
|
|
2020-02-20 13:07:54 +01:00
|
|
|
|
% Get first element of varargin if solve_algo ∈ {12,14} and rename varargin.
|
|
|
|
|
if ismember(options.solve_algo, [12, 14])
|
|
|
|
|
isloggedlhs = varargin{1};
|
2020-02-27 22:08:21 +01:00
|
|
|
|
isauxdiffloggedrhs = varargin{2};
|
|
|
|
|
endo_names = varargin{3};
|
|
|
|
|
lhs = varargin{4};
|
|
|
|
|
arguments = varargin(5:end);
|
2020-02-20 13:07:54 +01:00
|
|
|
|
else
|
|
|
|
|
arguments = varargin;
|
|
|
|
|
end
|
|
|
|
|
|
2012-10-31 11:08:22 +01:00
|
|
|
|
% checking initial values
|
|
|
|
|
if jacobian_flag
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[fvec, fjac] = feval(f, x, arguments{:});
|
2017-05-10 23:27:41 +02:00
|
|
|
|
wrong_initial_guess_flag = false;
|
2017-09-29 10:19:13 +02:00
|
|
|
|
if ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) ...
|
|
|
|
|
|| any(~isreal(fvec)) || any(~isreal(fjac(:)))
|
2019-12-12 18:50:46 +01:00
|
|
|
|
if max(abs(fvec)) < tolf %return if initial value solves problem
|
2019-12-20 16:28:06 +01:00
|
|
|
|
info = 0;
|
2019-12-12 18:50:46 +01:00
|
|
|
|
return;
|
2019-12-20 16:28:06 +01:00
|
|
|
|
end
|
2017-09-29 10:19:13 +02:00
|
|
|
|
disp('Randomize initial guess...')
|
2017-05-10 23:27:41 +02:00
|
|
|
|
% Let's try random numbers for the variables initialized with the default value.
|
|
|
|
|
wrong_initial_guess_flag = true;
|
|
|
|
|
% First try with positive numbers.
|
|
|
|
|
tentative_number = 0;
|
|
|
|
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
|
|
|
|
tentative_number = tentative_number+1;
|
|
|
|
|
x(idx) = rand(in0, 1)*10;
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[fvec, fjac] = feval(f, x, arguments{:});
|
2017-05-10 23:27:41 +02:00
|
|
|
|
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
|
|
|
|
|
end
|
|
|
|
|
% If all previous attempts failed, try with real numbers.
|
|
|
|
|
tentative_number = 0;
|
|
|
|
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
|
|
|
|
tentative_number = tentative_number+1;
|
|
|
|
|
x(idx) = randn(in0, 1)*10;
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[fvec, fjac] = feval(f, x, arguments{:});
|
2017-05-10 23:27:41 +02:00
|
|
|
|
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
|
|
|
|
|
end
|
|
|
|
|
% Last tentative, ff all previous attempts failed, try with negative numbers.
|
|
|
|
|
tentative_number = 0;
|
|
|
|
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
|
|
|
|
tentative_number = tentative_number+1;
|
|
|
|
|
x(idx) = -rand(in0, 1)*10;
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[fvec, fjac] = feval(f, x, arguments{:});
|
2017-05-10 23:27:41 +02:00
|
|
|
|
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
|
|
|
|
|
end
|
2013-06-20 18:53:35 +02:00
|
|
|
|
end
|
2012-10-31 11:08:22 +01:00
|
|
|
|
else
|
2020-02-20 13:07:54 +01:00
|
|
|
|
fvec = feval(f, x, arguments{:});
|
|
|
|
|
fjac = zeros(nn, nn);
|
2017-05-10 23:27:41 +02:00
|
|
|
|
wrong_initial_guess_flag = false;
|
|
|
|
|
if ~all(isfinite(fvec))
|
|
|
|
|
% Let's try random numbers for the variables initialized with the default value.
|
|
|
|
|
wrong_initial_guess_flag = true;
|
|
|
|
|
% First try with positive numbers.
|
|
|
|
|
tentative_number = 0;
|
|
|
|
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
|
|
|
|
tentative_number = tentative_number+1;
|
|
|
|
|
x(idx) = rand(in0, 1)*10;
|
2020-02-20 13:07:54 +01:00
|
|
|
|
fvec = feval(f, x, arguments{:});
|
2017-05-10 23:27:41 +02:00
|
|
|
|
wrong_initial_guess_flag = ~all(isfinite(fvec));
|
|
|
|
|
end
|
|
|
|
|
% If all previous attempts failed, try with real numbers.
|
|
|
|
|
tentative_number = 0;
|
|
|
|
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
|
|
|
|
tentative_number = tentative_number+1;
|
|
|
|
|
x(idx) = randn(in0, 1)*10;
|
2020-02-20 13:07:54 +01:00
|
|
|
|
fvec = feval(f, x, arguments{:});
|
2017-05-10 23:27:41 +02:00
|
|
|
|
wrong_initial_guess_flag = ~all(isfinite(fvec));
|
|
|
|
|
end
|
|
|
|
|
% Last tentative, ff all previous attempts failed, try with negative numbers.
|
|
|
|
|
tentative_number = 0;
|
|
|
|
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
|
|
|
|
tentative_number = tentative_number+1;
|
|
|
|
|
x(idx) = -rand(in0, 1)*10;
|
2020-02-20 13:07:54 +01:00
|
|
|
|
fvec = feval(f, x, arguments{:});
|
2017-05-10 23:27:41 +02:00
|
|
|
|
wrong_initial_guess_flag = ~all(isfinite(fvec));
|
|
|
|
|
end
|
|
|
|
|
end
|
2012-10-31 11:08:22 +01:00
|
|
|
|
end
|
|
|
|
|
|
2017-05-10 23:27:41 +02:00
|
|
|
|
% Exit with error if no initial guess has been found.
|
|
|
|
|
if wrong_initial_guess_flag
|
2020-02-20 13:07:54 +01:00
|
|
|
|
errorflag = true;
|
2014-05-26 16:16:22 +02:00
|
|
|
|
x = NaN(size(fvec));
|
2017-05-10 23:27:41 +02:00
|
|
|
|
return
|
2012-10-31 11:08:22 +01:00
|
|
|
|
end
|
|
|
|
|
|
2015-03-07 10:10:35 +01:00
|
|
|
|
% this test doesn't check complementarity conditions and is not used for
|
|
|
|
|
% mixed complementarity problems
|
2016-08-19 11:20:14 +02:00
|
|
|
|
if (~ismember(options.solve_algo,[10,11])) && (max(abs(fvec)) < tolf)
|
2012-10-31 11:08:22 +01:00
|
|
|
|
return ;
|
|
|
|
|
end
|
|
|
|
|
|
2015-05-25 10:06:13 +02:00
|
|
|
|
if options.solve_algo == 0
|
2013-11-04 10:54:45 +01:00
|
|
|
|
if ~isoctave
|
2012-01-09 11:06:26 +01:00
|
|
|
|
if ~user_has_matlab_license('optimization_toolbox')
|
2011-12-19 12:36:35 +01:00
|
|
|
|
error('You can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
|
|
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2015-05-29 15:19:20 +02:00
|
|
|
|
options4fsolve=optimset('fsolve');
|
|
|
|
|
options4fsolve.MaxFunEvals = 50000;
|
2016-08-02 17:06:51 +02:00
|
|
|
|
options4fsolve.MaxIter = maxit;
|
2015-05-29 15:19:20 +02:00
|
|
|
|
options4fsolve.TolFun = tolf;
|
2016-10-09 22:15:48 +02:00
|
|
|
|
if options.debug==1
|
|
|
|
|
options4fsolve.Display = 'final';
|
|
|
|
|
else
|
|
|
|
|
options4fsolve.Display = 'off';
|
|
|
|
|
end
|
2008-09-22 15:13:04 +02:00
|
|
|
|
if jacobian_flag
|
2015-07-07 11:11:28 +02:00
|
|
|
|
options4fsolve.Jacobian = 'on';
|
2008-09-22 15:13:04 +02:00
|
|
|
|
else
|
2015-07-07 11:11:28 +02:00
|
|
|
|
options4fsolve.Jacobian = 'off';
|
2008-09-22 15:13:04 +02:00
|
|
|
|
end
|
2013-11-04 10:54:45 +01:00
|
|
|
|
if ~isoctave
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[x, ~, exitval] = fsolve(f, x, options4fsolve, arguments{:});
|
2010-10-28 12:07:33 +02:00
|
|
|
|
else
|
|
|
|
|
% Under Octave, use a wrapper, since fsolve() does not have a 4th arg
|
2020-02-20 13:07:54 +01:00
|
|
|
|
if ischar(f)
|
|
|
|
|
f2 = str2func(f);
|
2016-08-12 11:59:57 +02:00
|
|
|
|
else
|
2020-02-20 13:07:54 +01:00
|
|
|
|
f2 = f;
|
2016-08-12 11:59:57 +02:00
|
|
|
|
end
|
2020-02-20 13:07:54 +01:00
|
|
|
|
f = @(x) f2(x, arguments{:});
|
2011-09-19 12:36:37 +02:00
|
|
|
|
% The Octave version of fsolve does not converge when it starts from the solution
|
2020-02-20 13:07:54 +01:00
|
|
|
|
fvec = feval(f, x);
|
2012-10-31 11:08:22 +01:00
|
|
|
|
if max(abs(fvec)) >= tolf
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[x, ~,exitval] = fsolve(f, x, options4fsolve);
|
2011-01-14 12:27:13 +01:00
|
|
|
|
else
|
|
|
|
|
exitval = 3;
|
2017-05-16 12:42:01 +02:00
|
|
|
|
end
|
2010-10-28 12:07:33 +02:00
|
|
|
|
end
|
2014-09-10 16:30:17 +02:00
|
|
|
|
|
2015-11-16 11:11:58 +01:00
|
|
|
|
if exitval == 1
|
2020-02-20 13:07:54 +01:00
|
|
|
|
errorflag = false;
|
2015-11-16 11:11:58 +01:00
|
|
|
|
elseif exitval > 1
|
2020-02-20 13:07:54 +01:00
|
|
|
|
if ischar(f)
|
|
|
|
|
f2 = str2func(f);
|
2016-08-12 11:59:57 +02:00
|
|
|
|
else
|
2020-02-20 13:07:54 +01:00
|
|
|
|
f2 = f;
|
2016-08-12 11:59:57 +02:00
|
|
|
|
end
|
2020-02-20 13:07:54 +01:00
|
|
|
|
f = @(x) f2(x, arguments{:});
|
|
|
|
|
fvec = feval(f, x);
|
2016-05-05 10:46:39 +02:00
|
|
|
|
if max(abs(fvec)) >= tolf
|
2020-02-20 13:07:54 +01:00
|
|
|
|
errorflag = true;
|
2015-11-16 11:11:58 +01:00
|
|
|
|
else
|
2020-02-20 13:07:54 +01:00
|
|
|
|
errorflag = false;
|
2015-11-16 11:11:58 +01:00
|
|
|
|
end
|
2008-09-22 15:13:04 +02:00
|
|
|
|
else
|
2020-02-20 13:07:54 +01:00
|
|
|
|
errorflag = true;
|
2008-09-22 15:13:04 +02:00
|
|
|
|
end
|
2020-02-20 13:07:54 +01:00
|
|
|
|
elseif options.solve_algo==1
|
|
|
|
|
[x, errorflag] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, ...
|
2021-04-08 13:28:07 +02:00
|
|
|
|
tolf, tolx, ...
|
2021-07-23 22:16:43 +02:00
|
|
|
|
maxit, [], options.debug, arguments{:});
|
2020-02-20 13:07:54 +01:00
|
|
|
|
elseif options.solve_algo==9
|
2021-07-22 16:46:10 +02:00
|
|
|
|
[x, errorflag] = trust_region(f,x, 1:nn, 1:nn, jacobian_flag, options.gstep, ...
|
|
|
|
|
tolf, tolx, maxit, ...
|
|
|
|
|
options.trust_region_initial_step_bound_factor, ...
|
|
|
|
|
options.debug, arguments{:});
|
2020-07-16 18:20:07 +02:00
|
|
|
|
elseif ismember(options.solve_algo, [2, 12, 4])
|
2020-02-20 13:07:54 +01:00
|
|
|
|
if ismember(options.solve_algo, [2, 12])
|
2014-02-04 17:55:55 +01:00
|
|
|
|
solver = @solve1;
|
|
|
|
|
else
|
|
|
|
|
solver = @trust_region;
|
|
|
|
|
end
|
2020-07-16 18:20:07 +02:00
|
|
|
|
specializedunivariateblocks = options.solve_algo == 12;
|
2006-01-18 17:50:33 +01:00
|
|
|
|
if ~jacobian_flag
|
2009-12-16 18:17:34 +01:00
|
|
|
|
fjac = zeros(nn,nn) ;
|
2020-02-20 13:07:54 +01:00
|
|
|
|
dh = max(abs(x), options.gstep(1)*ones(nn,1))*eps^(1/3);
|
2009-12-16 18:17:34 +01:00
|
|
|
|
for j = 1:nn
|
|
|
|
|
xdh = x ;
|
|
|
|
|
xdh(j) = xdh(j)+dh(j) ;
|
2020-02-20 13:07:54 +01:00
|
|
|
|
fjac(:,j) = (feval(f, xdh, arguments{:})-fvec)./dh(j) ;
|
2009-12-16 18:17:34 +01:00
|
|
|
|
end
|
2006-01-18 17:50:33 +01:00
|
|
|
|
end
|
2020-02-27 22:08:21 +01:00
|
|
|
|
[j1,j2,r,s] = dmperm(fjac);
|
2020-03-11 12:38:59 +01:00
|
|
|
|
JAC = abs(fjac(j1,j2))>0;
|
2015-05-25 10:06:13 +02:00
|
|
|
|
if options.debug
|
2020-07-16 18:20:07 +02:00
|
|
|
|
disp(['DYNARE_SOLVE (solve_algo=2|4|12): number of blocks = ' num2str(length(r)-1)]);
|
2008-09-16 19:24:11 +02:00
|
|
|
|
end
|
2020-02-27 22:08:21 +01:00
|
|
|
|
l = 0;
|
|
|
|
|
fre = false;
|
2005-02-18 20:54:39 +01:00
|
|
|
|
for i=length(r)-1:-1:1
|
2020-02-20 13:07:54 +01:00
|
|
|
|
blocklength = r(i+1)-r(i);
|
2015-05-25 10:06:13 +02:00
|
|
|
|
if options.debug
|
2020-07-16 18:20:07 +02:00
|
|
|
|
dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u of size %u.', i, blocklength);
|
2009-12-16 18:17:34 +01:00
|
|
|
|
end
|
2020-02-20 13:07:54 +01:00
|
|
|
|
j = r(i):r(i+1)-1;
|
|
|
|
|
if specializedunivariateblocks
|
|
|
|
|
if options.debug
|
2020-07-16 18:20:07 +02:00
|
|
|
|
dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u by evaluating RHS.', i);
|
2020-02-20 13:07:54 +01:00
|
|
|
|
end
|
|
|
|
|
if isequal(blocklength, 1)
|
2020-02-27 22:08:21 +01:00
|
|
|
|
if i<length(r)-1
|
|
|
|
|
if fre || any(JAC(r(i), s(i)+(1:l)))
|
|
|
|
|
% Reevaluation of the residuals is required because the current RHS depends on
|
|
|
|
|
% variables that potentially have been updated previously.
|
|
|
|
|
z = feval(f, x, arguments{:});
|
|
|
|
|
l = 0;
|
|
|
|
|
fre = false;
|
|
|
|
|
end
|
2020-02-20 13:07:54 +01:00
|
|
|
|
else
|
2020-02-27 22:08:21 +01:00
|
|
|
|
% First iteration requires the evaluation of the residuals.
|
|
|
|
|
z = feval(f, x, arguments{:});
|
|
|
|
|
end
|
|
|
|
|
l = l+1;
|
|
|
|
|
if isequal(lhs{j1(j)}, endo_names{j2(j)}) || isequal(lhs{j1(j)}, sprintf('log(%s)', endo_names{j2(j)}))
|
|
|
|
|
if isloggedlhs(j1(j))
|
|
|
|
|
x(j2(j)) = exp(log(x(j2(j)))-z(j1(j)));
|
|
|
|
|
else
|
|
|
|
|
x(j2(j)) = x(j2(j))-z(j1(j));
|
|
|
|
|
end
|
|
|
|
|
else
|
|
|
|
|
if options.debug
|
|
|
|
|
dprintf('LHS variable is not determined by RHS expression (%u).', j1(j))
|
|
|
|
|
dprintf('%s -> %s', lhs{j1(j)}, endo_names{j2(j)})
|
|
|
|
|
end
|
|
|
|
|
if ~isempty(regexp(lhs{j1(j)}, '\<AUX_DIFF_(\d*)\>', 'once'))
|
|
|
|
|
if isauxdiffloggedrhs(j1(j))
|
|
|
|
|
x(j2(j)) = exp(log(x(j2(j)))+z(j1(j)));
|
|
|
|
|
else
|
|
|
|
|
x(j2(j)) = x(j2(j))+z(j1(j));
|
|
|
|
|
end
|
|
|
|
|
else
|
|
|
|
|
error('Algorithm solve_algo=%u cannot be used with this nonlinear problem.', options.solve_algo)
|
|
|
|
|
end
|
2020-02-20 13:07:54 +01:00
|
|
|
|
end
|
|
|
|
|
continue
|
|
|
|
|
end
|
|
|
|
|
else
|
|
|
|
|
if options.debug
|
2020-07-16 18:20:07 +02:00
|
|
|
|
dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u with trust_region routine.', i);
|
2020-02-20 13:07:54 +01:00
|
|
|
|
end
|
|
|
|
|
end
|
|
|
|
|
[x, errorflag] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
|
2020-02-27 22:08:21 +01:00
|
|
|
|
options.gstep, ...
|
2021-07-22 16:46:10 +02:00
|
|
|
|
tolf, options.solve_tolx, maxit, ...
|
|
|
|
|
options.trust_region_initial_step_bound_factor, ...
|
|
|
|
|
options.debug, arguments{:});
|
2020-02-27 22:08:21 +01:00
|
|
|
|
fre = true;
|
2020-02-20 13:07:54 +01:00
|
|
|
|
if errorflag
|
2009-12-16 18:17:34 +01:00
|
|
|
|
return
|
|
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
|
end
|
2020-02-20 13:07:54 +01:00
|
|
|
|
fvec = feval(f, x, arguments{:});
|
|
|
|
|
if max(abs(fvec))>tolf
|
2020-03-11 12:47:02 +01:00
|
|
|
|
disp('Call solver on the full nonlinear problem.')
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[x, errorflag] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
|
2021-07-22 16:46:10 +02:00
|
|
|
|
options.gstep, tolf, options.solve_tolx, maxit, ...
|
|
|
|
|
options.trust_region_initial_step_bound_factor, ...
|
|
|
|
|
options.debug, arguments{:});
|
2007-01-04 15:42:27 +01:00
|
|
|
|
end
|
2020-02-20 13:07:54 +01:00
|
|
|
|
elseif options.solve_algo==3
|
2007-11-27 15:16:13 +01:00
|
|
|
|
if jacobian_flag
|
2021-06-02 16:21:16 +02:00
|
|
|
|
[x, errorflag] = csolve(f, x, f, tolf, maxit, arguments{:});
|
2007-11-27 15:16:13 +01:00
|
|
|
|
else
|
2021-06-02 16:21:16 +02:00
|
|
|
|
[x, errorflag] = csolve(f, x, [], tolf, maxit, arguments{:});
|
2008-09-22 15:13:04 +02:00
|
|
|
|
end
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[fvec, fjac] = feval(f, x, arguments{:});
|
|
|
|
|
elseif options.solve_algo==10
|
2016-05-05 20:40:14 +02:00
|
|
|
|
% LMMCP
|
2015-05-25 10:06:13 +02:00
|
|
|
|
olmmcp = options.lmmcp;
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[x, ~, exitflag] = lmmcp(f, x, olmmcp.lb, olmmcp.ub, olmmcp, arguments{:});
|
|
|
|
|
if exitflag==1
|
|
|
|
|
errorflag = false;
|
2014-05-12 17:15:49 +02:00
|
|
|
|
else
|
2020-02-20 13:07:54 +01:00
|
|
|
|
errorflag = true;
|
2014-05-12 17:15:49 +02:00
|
|
|
|
end
|
2016-05-05 20:40:14 +02:00
|
|
|
|
elseif options.solve_algo == 11
|
|
|
|
|
% PATH mixed complementary problem
|
|
|
|
|
% PATH linear mixed complementary problem
|
|
|
|
|
if ~exist('mcppath')
|
|
|
|
|
error(['PATH can''t be provided with Dynare. You need to install it ' ...
|
|
|
|
|
'yourself and add its location to Matlab/Octave path before ' ...
|
|
|
|
|
'running Dynare'])
|
|
|
|
|
end
|
|
|
|
|
omcppath = options.mcppath;
|
|
|
|
|
global mcp_data
|
2020-02-20 13:07:54 +01:00
|
|
|
|
mcp_data.func = f;
|
|
|
|
|
mcp_data.args = arguments;
|
2016-07-04 10:29:45 +02:00
|
|
|
|
try
|
2020-02-20 13:07:54 +01:00
|
|
|
|
[x, fval, jac, mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
|
2016-07-04 10:29:45 +02:00
|
|
|
|
catch
|
2020-02-20 13:07:54 +01:00
|
|
|
|
errorflag = true;
|
2016-07-04 10:29:45 +02:00
|
|
|
|
end
|
2020-07-16 18:20:07 +02:00
|
|
|
|
elseif ismember(options.solve_algo, [13, 14])
|
|
|
|
|
if ~jacobian_flag
|
|
|
|
|
error('DYNARE_SOLVE: option solve_algo=13|14 needs computed Jacobian')
|
|
|
|
|
end
|
|
|
|
|
auxstruct = struct();
|
|
|
|
|
if options.solve_algo == 14
|
|
|
|
|
auxstruct.lhs = lhs;
|
|
|
|
|
auxstruct.endo_names = endo_names;
|
|
|
|
|
auxstruct.isloggedlhs = isloggedlhs;
|
|
|
|
|
auxstruct.isauxdiffloggedrhs = isauxdiffloggedrhs;
|
|
|
|
|
end
|
2021-07-22 16:13:35 +02:00
|
|
|
|
[x, errorflag] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, auxstruct, arguments{:});
|
2020-07-16 18:20:07 +02:00
|
|
|
|
[fvec, fjac] = feval(f, x, arguments{:});
|
2009-12-16 18:17:34 +01:00
|
|
|
|
else
|
2020-07-16 18:20:07 +02:00
|
|
|
|
error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11,12,13,14]')
|
2009-12-16 18:17:34 +01:00
|
|
|
|
end
|