diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 2619c4cc256e5f0e9597bb3e6942af0e9f301ca1..2dfc5735becd8b8a041877369181c3215c18aa6f 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -61,7 +61,7 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
end
elseif options_.solve_algo == 1
nn = size(x,1);
- [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
+ [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,1,varargin{:});
elseif options_.solve_algo == 2 || options_.solve_algo == 4
nn = size(x,1) ;
tolf = options_.solve_tolf ;
@@ -101,27 +101,22 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
if options_.debug
disp(['DYNARE_SOLVE (solve_algo=2|4): number of blocks = ' num2str(length(r))]);
end
+
+ % Activate bad conditioning flag for solve_algo = 2, but not for solve_algo = 4
+ bad_cond_flag = (options_.solve_algo == 2);
for i=length(r)-1:-1:1
if options_.debug
disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]);
end
- if options_.solve_algo == 2
- [x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag,varargin{:});
- else % solve_algo=4
- [x,info]=solve2(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag,varargin{:});
- end
+ [x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, bad_cond_flag, varargin{:});
if info
return
end
end
fvec = feval(func,x,varargin{:});
if max(abs(fvec)) > tolf
- if options_.solve_algo == 2
- [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
- else % solve_algo=4
- [x,info]=solve2(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
- end
+ [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag, bad_cond_flag, varargin{:});
end
elseif options_.solve_algo == 3
if jacobian_flag
diff --git a/matlab/solve1.m b/matlab/solve1.m
index 9f6f84cbebd89ddd442f906dc4e305f1a28af5dd..e451eccbf846fe16e63d03114bad06d609e26a61 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -1,5 +1,5 @@
-function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin)
-% function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin)
+function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
+% function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
% Solves systems of non linear equations of several variables
%
% INPUTS
@@ -9,7 +9,9 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin)
% j2: unknown variables index
% jacobian_flag=1: jacobian given by the 'func' function
% jacobian_flag=0: jacobian obtained numerically
-% varargin: list of arguments following jacobian_flag
+% bad_cond_flag=1: when Jacobian is badly conditionned, use an
+% alternative formula to Newton step
+% varargin: list of arguments following bad_cond_flag
%
% OUTPUTS
% x: results
@@ -18,7 +20,7 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin)
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2001-2008 Dynare Team
+% Copyright (C) 2001-2009 Dynare Team
%
% This file is part of Dynare.
%
@@ -113,8 +115,7 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin)
fvec = q'*fvec;
p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
end
-% elseif cond(fjac) > 10*sqrt(eps)
- elseif cond(fjac) > 1/sqrt(eps)
+ elseif bad_cond_flag && cond(fjac) > 1/sqrt(eps)
fjac2=fjac'*fjac;
p=-(fjac2+sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec);
else
diff --git a/matlab/solve2.m b/matlab/solve2.m
deleted file mode 100644
index db0012fcd6e45863ea96476145a61222010dab01..0000000000000000000000000000000000000000
--- a/matlab/solve2.m
+++ /dev/null
@@ -1,171 +0,0 @@
-function [x,check] = solve2(func,x,j1,j2,jacobian_flag,varargin)
-% function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin)
-% Solves systems of non linear equations of several variables
-%
-% This solver is used for solve_algo = ...
-%
-% This is a modified version of solve1:
-% In solve1, before proceeding to the Newton step, we first check
-% that the condition number is reasonable, and we use an alternate
-% step formula if it is not the case.
-% Here, we first do the Newton step, then we check if the left
-% matrix division returned a warning (in case of badly scale or
-% nearly singular jacobian) in which case we use the alternate
-% step formula.
-%
-% INPUTS
-% func: name of the function to be solved
-% x: guess values
-% j1: equations index for which the model is solved
-% j2: unknown variables index
-% jacobian_flag=1: jacobian given by the 'func' function
-% jacobian_flag=0: jacobian obtained numerically
-% varargin: list of arguments following jacobian_flag
-%
-% OUTPUTS
-% x: results
-% check=1: the model can not be solved
-%
-% SPECIAL REQUIREMENTS
-% none
-
-% Copyright (C) 2001-2008 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 <http://www.gnu.org/licenses/>.
-
- global M_ options_ fjac
-
- nn = length(j1);
- fjac = zeros(nn,nn) ;
- g = zeros(nn,1) ;
-
- tolf = options_.solve_tolf ;
- tolx = options_.solve_tolx;
- tolmin = tolx ;
-
- stpmx = 100 ;
- maxit = options_.solve_maxit ;
-
- check = 0 ;
-
- fvec = feval(func,x,varargin{:});
- fvec = fvec(j1);
-
- i = find(~isfinite(fvec));
-
- if ~isempty(i)
- disp(['STEADY: numerical initial values incompatible with the following' ...
- ' equations'])
- disp(j1(i)')
- end
-
- f = 0.5*fvec'*fvec ;
-
- if max(abs(fvec)) < tolf
- return ;
- end
-
- stpmax = stpmx*max([sqrt(x'*x);nn]) ;
- first_time = 1;
- for its = 1:maxit
- if jacobian_flag
- [fvec,fjac] = feval(func,x,varargin{:});
- fvec = fvec(j1);
- fjac = fjac(j1,j2);
- else
- dh = max(abs(x(j2)),options_.gstep*ones(nn,1))*eps^(1/3);
-
- for j = 1:nn
- xdh = x ;
- xdh(j2(j)) = xdh(j2(j))+dh(j) ;
- t = feval(func,xdh,varargin{:});
- fjac(:,j) = (t(j1) - fvec)./dh(j) ;
- g(j) = fvec'*fjac(:,j) ;
- end
- end
-
- g = (fvec'*fjac)';
- if options_.debug
- disp(['cond(fjac) ' num2str(cond(fjac))])
- end
- M_.unit_root = 0;
- if M_.unit_root
- if first_time
- first_time = 0;
- [q,r,e]=qr(fjac);
- n = sum(abs(diag(r)) < 1e-12);
- fvec = q'*fvec;
- p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
- disp(' ')
- disp('STEADY with unit roots:')
- disp(' ')
- if n > 0
- disp([' The following variable(s) kept their value given in INITVAL' ...
- ' or ENDVAL'])
- disp(char(e(:,end-n+1:end)'*M_.endo_names))
- else
- disp(' STEADY can''t find any unit root!')
- end
- else
- [q,r]=qr(fjac*e);
- fvec = q'*fvec;
- p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
- end
- else
- lastwarn('');
- p = -fjac\fvec;
- if ~isempty(lastwarn)
- fjac2=fjac'*fjac;
- p=-(fjac2+sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec);
- end
- end
- xold = x ;
- fold = f ;
-
- [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin{:});
-
- if options_.debug
- disp([its f])
- disp([xold x])
- end
-
- if check > 0
- den = max([f;0.5*nn]) ;
- if max(abs(g).*max([abs(x(j2)') ones(1,nn)])')/den < tolmin
- return
- else
- disp (' ')
- disp (['SOLVE: Iteration ' num2str(its)])
- disp (['Spurious convergence.'])
- disp (x)
- return
- end
-
- if max(abs(x(j2)-xold(j2))./max([abs(x(j2)') ones(1,nn)])') < tolx
- disp (' ')
- disp (['SOLVE: Iteration ' num2str(its)])
- disp (['Convergence on dX.'])
- disp (x)
- return
- end
- elseif max(abs(fvec)) < tolf
- return
- end
- end
-
- check = 1;
- disp(' ')
- disp('SOLVE: maxit has been reached')