2020-06-02 17:45:28 +02:00
function [ y , T , oo_ , info ] = solve_one_boundary ( fname , y , x , params , steady_state , T , ...
y_index_eq , nze , periods , is_linear , Block_Num , y_kmin , maxit_ , solve_tolf , lambda , cutoff , stack_solve_algo , is_forward , is_dynamic , verbose , M , options , oo_ )
2009-12-16 18:17:34 +01:00
% Computes the deterministic simulation of a block of equation containing
2017-05-16 15:10:20 +02:00
% lead or lag variables
2009-12-16 18:17:34 +01:00
%
% INPUTS
% fname [string] name of the file containing the block
% to simulate
% y [matrix] All the endogenous variables of the model
% x [matrix] All the exogenous variables of the model
% params [vector] All the parameters of the model
2011-03-13 17:06:57 +01:00
% steady_state [vector] steady state of the model
2020-05-25 18:36:55 +02:00
% T [matrix] Temporary terms
2009-12-16 18:17:34 +01:00
% y_index_eq [vector of int] The index of the endogenous variables of
% the block
% nze [integer] number of non-zero elements in the
% jacobian matrix
% periods [integer] number of simulation periods
2020-05-25 16:58:11 +02:00
% is_linear [logical] whether the block is linear
2009-12-16 18:17:34 +01:00
% Block_Num [integer] block number
% y_kmin [integer] maximum number of lag in the model
% maxit_ [integer] maximum number of iteration in Newton
% solve_tolf [double] convergence criteria
% lambda [double] initial value of step size in
% Newton
% cutoff [double] cutoff to correct the direction in Newton in case
% of singular jacobian matrix
% stack_solve_algo [integer] linear solver method used in the
2017-05-16 15:10:20 +02:00
% Newton algorithm :
2009-12-16 18:17:34 +01:00
% - 1 sparse LU
% - 2 GMRES
% - 3 BicGStab
% - 4 Optimal path length
2020-05-25 16:58:11 +02:00
% is_forward [logical] Whether the block has to be solved forward
% If false, the block is solved backward
% is_dynamic [logical] If true, the block belongs to the dynamic file
2009-12-16 18:17:34 +01:00
% file and the oo_.deterministic_simulation field has to be uptated
2020-05-25 16:58:11 +02:00
% If false, the block belongs to the static
% file and the oo_.detereministic_simulation
2009-12-16 18:17:34 +01:00
% field remains unchanged
2020-05-25 16:58:11 +02:00
% verbose [logical] Whether iterations are to be printed
%
2017-05-16 15:10:20 +02:00
% OUTPUTS
% y [matrix] All endogenous variables of the model
2020-05-26 16:23:46 +02:00
% T [matrix] Temporary terms
2010-10-18 17:28:21 +02:00
% info [integer] >=0 no error
% <0 error
2017-05-16 15:10:20 +02:00
%
2009-12-16 18:17:34 +01:00
% ALGORITHM
% Newton with LU or GMRES or BicGstab for dynamic block
2017-05-16 15:10:20 +02:00
%
2009-12-16 18:17:34 +01:00
% SPECIAL REQUIREMENTS
% none.
2017-05-16 15:10:20 +02:00
%
2009-12-16 18:17:34 +01:00
2022-02-11 16:47:58 +01:00
% Copyright (C) 1996-2022 Dynare Team
2009-12-16 18:17:34 +01: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/>.
2009-12-16 18:17:34 +01:00
Blck_size = size ( y_index_eq , 2 ) ;
correcting_factor = 0.01 ;
For stack_solve_algo={2,3} and solve_algo={7,8}, rather use a “crout” ILU decomposition for the preconditioner
In 3025a14ed9074c684a96b428be939b60a3c2d15a, the call to the deprecated “luinc”
has been replaced by a call to “ilu”.
However, the type of “ilu” decomposition had not been specified. The default,
“nofill”, does not make use of the “droptol” option which was used with
“luinc”. Also, under Octave, it does not work when there is a zero on the
diagonal of the input matrix.
Rather use the “crout” type, which addresses these two issues.
2022-02-14 12:34:50 +01:00
ilu_setup . type = ' crout' ;
2013-03-25 12:03:09 +01:00
ilu_setup . droptol = 1e-10 ;
2009-12-16 18:17:34 +01:00
max_resa = 1e100 ;
reduced = 0 ;
2020-05-25 16:58:11 +02:00
if is_forward
2009-12-16 18:17:34 +01:00
incr = 1 ;
start = y_kmin + 1 ;
finish = periods + y_kmin ;
else
incr = - 1 ;
start = periods + y_kmin ;
finish = y_kmin + 1 ;
end
2015-02-20 14:44:45 +01:00
2009-12-16 18:17:34 +01:00
for it_ = start : incr : finish
cvg = 0 ;
iter = 0 ;
g1 = spalloc ( Blck_size , Blck_size , nze ) ;
2015-02-20 14:44:45 +01:00
while ~ ( cvg == 1 || iter > maxit_ )
if is_dynamic
2020-10-06 18:31:13 +02:00
[ r , yy , T ( : , it_ ) , g1 ] = feval ( fname , Block_Num , dynvars_from_endo_simul ( y , it_ , M ) , x , params , steady_state , T ( : , it_ ) , it_ , false ) ;
y ( : , it_ ) = yy ( M . lead_lag_incidence ( M . maximum_endo_lag + 1 , : ) ) ;
2009-12-16 18:17:34 +01:00
else
2020-10-06 18:31:13 +02:00
[ r , y , T , g1 ] = feval ( fname , Block_Num , y , x , params , T ) ;
2015-02-20 14:44:45 +01:00
end
if ~ isreal ( r )
2009-12-16 18:17:34 +01:00
max_res = ( - ( max ( max ( abs ( r ) ) ) ) ^2 ) ^0.5 ;
else
max_res = max ( max ( abs ( r ) ) ) ;
2015-02-20 14:44:45 +01:00
end
2020-05-25 16:58:11 +02:00
if verbose
2015-02-20 14:44:45 +01:00
disp ( [ ' iteration : ' int2str ( iter + 1 ) ' => ' num2str ( max_res ) ' time = ' int2str ( it_ ) ] )
if is_dynamic
2020-06-17 17:28:59 +02:00
disp ( [ M . endo_names { y_index_eq } num2str ( [ y ( y_index_eq , it_ ) r g1 ] ) ] )
2009-12-16 18:17:34 +01:00
else
2017-10-10 10:05:59 +02:00
disp ( [ M . endo_names { y_index_eq } num2str ( [ y ( y_index_eq ) r g1 ] ) ] )
2015-02-20 14:44:45 +01:00
end
end
if ~ isreal ( max_res ) || isnan ( max_res )
2009-12-16 18:17:34 +01:00
cvg = 0 ;
2015-02-20 14:44:45 +01:00
elseif is_linear && iter > 0
2009-12-16 18:17:34 +01:00
cvg = 1 ;
else
cvg = ( max_res < solve_tolf ) ;
2015-02-20 14:44:45 +01:00
end
if ~ cvg
if iter > 0
if ~ isreal ( max_res ) || isnan ( max_res ) || ( max_resa < max_res && iter > 1 )
if isnan ( max_res ) || ( max_resa < max_res && iter > 0 )
2009-12-16 18:17:34 +01:00
detJ = det ( g1a ) ;
if ( abs ( detJ ) < 1e-7 )
max_factor = max ( max ( abs ( g1a ) ) ) ;
ze_elem = sum ( diag ( g1a ) < cutoff ) ;
2015-02-20 14:44:45 +01:00
if verbose
disp ( [ num2str ( full ( ze_elem ) , ' %d' ) ' elements on the Jacobian diagonal are below the cutoff (' num2str ( cutoff , ' %f' ) ' )' ] )
end
if correcting_factor < max_factor
2009-12-16 18:17:34 +01:00
correcting_factor = correcting_factor * 4 ;
2015-02-20 14:44:45 +01:00
if verbose
disp ( [ ' The Jacobain matrix is singular, det(Jacobian)=' num2str ( detJ , ' %f' ) ' .' ] )
disp ( [ ' trying to correct the Jacobian matrix:' ] )
disp ( [ ' correcting_factor=' num2str ( correcting_factor , ' %f' ) ' max(Jacobian)=' num2str ( full ( max_factor ) , ' %f' ) ] )
end
2009-12-16 18:17:34 +01:00
dx = - r / ( g1 + correcting_factor * speye ( Blck_size ) ) ;
2020-06-17 17:28:59 +02:00
y ( y_index_eq , it_ ) = ya_save + lambda * dx ;
2015-02-20 14:44:45 +01:00
continue
2009-12-16 18:17:34 +01:00
else
2015-02-20 14:44:45 +01:00
if verbose
disp ( ' The singularity of the jacobian matrix could not be corrected' )
end
2009-12-16 18:17:34 +01:00
info = - Block_Num * 10 ;
2015-02-20 14:44:45 +01:00
return
end
end
elseif lambda > 1e-8
2009-12-16 18:17:34 +01:00
lambda = lambda / 2 ;
reduced = 1 ;
2015-02-20 14:44:45 +01:00
if verbose
disp ( [ ' reducing the path length: lambda=' num2str ( lambda , ' %f' ) ] )
end
if is_dynamic
2020-06-17 17:28:59 +02:00
y ( y_index_eq , it_ ) = ya_save - lambda * dx ;
2009-12-16 18:17:34 +01:00
else
y ( y_index_eq ) = ya_save - lambda * dx ;
2015-02-20 14:44:45 +01:00
end
continue
2009-12-16 18:17:34 +01:00
else
2015-02-20 14:44:45 +01:00
if verbose
if cutoff == 0
fprintf ( ' Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit".\n' , Block_Num , it_ , iter ) ;
else
fprintf ( ' Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n' , Block_Num , it_ , iter ) ;
end
end
if is_dynamic
2020-06-02 17:45:28 +02:00
oo_ . deterministic_simulation . status = false ;
2009-12-16 18:17:34 +01:00
oo_ . deterministic_simulation . error = max_res ;
oo_ . deterministic_simulation . iterations = iter ;
2020-06-02 17:45:28 +02:00
oo_ . deterministic_simulation . block ( Block_Num ) . status = false ; % Convergency failed.
2009-12-16 18:17:34 +01:00
oo_ . deterministic_simulation . block ( Block_Num ) . error = max_res ;
oo_ . deterministic_simulation . block ( Block_Num ) . iterations = iter ;
2015-02-20 14:44:45 +01:00
end
2009-12-16 18:17:34 +01:00
info = - Block_Num * 10 ;
2015-02-20 14:44:45 +01:00
return
end
2009-12-16 18:17:34 +01:00
else
2015-02-20 14:44:45 +01:00
if lambda < 1
2009-12-16 18:17:34 +01:00
lambda = max ( lambda * 2 , 1 ) ;
2015-02-20 14:44:45 +01:00
end
end
end
if is_dynamic
2020-06-17 17:28:59 +02:00
ya = y ( y_index_eq , it_ ) ;
2009-12-16 18:17:34 +01:00
else
ya = y ( y_index_eq ) ;
2015-02-20 14:44:45 +01:00
end
2009-12-16 18:17:34 +01:00
ya_save = ya ;
g1a = g1 ;
2020-05-25 17:29:24 +02:00
if is_dynamic && stack_solve_algo == 4
2009-12-16 18:17:34 +01:00
stpmx = 100 ;
2022-02-11 16:47:58 +01:00
stpmax = stpmx * max ( [ sqrt ( ya ' * ya ) ; size ( y_index_eq , 2 ) ] ) ;
2009-12-16 18:17:34 +01:00
nn = 1 : size ( y_index_eq , 2 ) ;
g = ( r ' * g1 ) ' ;
f = 0.5 * r ' * r ;
p = - g1 \ r ;
2020-06-16 18:14:09 +02:00
[ ya , f , r , check ] = lnsrch1 ( ya , f , g , p , stpmax , ...
2020-05-25 17:29:24 +02:00
' lnsrch1_wrapper_one_boundary' , nn , ...
2022-02-11 16:47:58 +01:00
nn , options . solve_tolx , y_index_eq , fname , Block_Num , y , x , params , steady_state , T ( : , it_ ) , it_ , M ) ;
dx = ya - y ( y_index_eq , it_ ) ;
y ( y_index_eq , it_ ) = ya ;
2020-05-26 16:23:46 +02:00
%% Recompute temporary terms, since they are not given as output of lnsrch1
2020-06-19 18:34:48 +02:00
[ ~ , ~ , T ( : , it_ ) ] = feval ( fname , Block_Num , dynvars_from_endo_simul ( y , it_ , M ) , x , params , steady_state , T ( : , it_ ) , it_ , false ) ;
2020-05-25 17:29:24 +02:00
elseif ( is_dynamic && ( stack_solve_algo == 1 || stack_solve_algo == 0 ) ) || ( ~ is_dynamic && options . solve_algo == 6 )
2020-05-25 16:58:11 +02:00
if verbose && ~ is_dynamic
2015-02-20 14:44:45 +01:00
disp ( ' steady: Sparse LU ' )
2011-02-04 17:17:48 +01:00
end
2009-12-16 18:17:34 +01:00
dx = g1 \ r ;
ya = ya - lambda * dx ;
2015-02-20 14:44:45 +01:00
if is_dynamic
2020-06-17 17:28:59 +02:00
y ( y_index_eq , it_ ) = ya ;
2009-12-16 18:17:34 +01:00
else
y ( y_index_eq ) = ya ;
2015-02-20 14:44:45 +01:00
end
elseif ( stack_solve_algo == 2 && is_dynamic ) || ( options . solve_algo == 7 && ~ is_dynamic )
2009-12-16 18:17:34 +01:00
flag1 = 1 ;
2020-05-25 16:58:11 +02:00
if verbose && ~ is_dynamic
2015-02-20 14:44:45 +01:00
disp ( ' steady: GMRES ' )
2011-02-04 17:17:48 +01:00
end
2015-02-20 14:44:45 +01:00
while flag1 > 0
2013-03-25 12:03:09 +01:00
[ L1 , U1 ] = ilu ( g1 , ilu_setup ) ;
2009-12-21 11:29:21 +01:00
[ dx , flag1 ] = gmres ( g1 , - r , Blck_size , 1e-6 , Blck_size , L1 , U1 ) ;
2015-02-20 14:44:45 +01:00
if flag1 > 0 || reduced
if verbose
if flag1 == 1
disp ( [ ' Error in simul: No convergence inside GMRES after ' num2str ( iter , ' %6d' ) ' iterations, in block' num2str ( Block_Num , ' %3d' ) ] )
elseif ( flag1 == 2 )
disp ( [ ' Error in simul: Preconditioner is ill-conditioned, in block' num2str ( Block_Num , ' %3d' ) ] )
elseif ( flag1 == 3 )
disp ( [ ' Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str ( Block_Num , ' %3d' ) ] )
end
end
2013-03-25 12:03:09 +01:00
ilu_setup . droptol = ilu_setup . droptol / 10 ;
2009-12-16 18:17:34 +01:00
reduced = 0 ;
else
ya = ya + lambda * dx ;
2015-02-20 14:44:45 +01:00
if is_dynamic
2020-06-17 17:28:59 +02:00
y ( y_index_eq , it_ ) = ya ;
2009-12-16 18:17:34 +01:00
else
y ( y_index_eq ) = ya ' ;
2015-02-20 14:44:45 +01:00
end
end
end
elseif ( stack_solve_algo == 3 && is_dynamic ) || ( options . solve_algo == 8 && ~ is_dynamic )
2009-12-16 18:17:34 +01:00
flag1 = 1 ;
2020-05-25 16:58:11 +02:00
if verbose && ~ is_dynamic
2015-02-20 14:44:45 +01:00
disp ( ' steady: BiCGStab' )
2011-02-04 17:17:48 +01:00
end
2015-02-20 14:44:45 +01:00
while flag1 > 0
2013-03-25 12:03:09 +01:00
[ L1 , U1 ] = ilu ( g1 , ilu_setup ) ;
2011-02-04 17:17:48 +01:00
phat = ya - U1 \ ( L1 \ r ) ;
2015-02-20 14:44:45 +01:00
if is_dynamic
2020-06-17 17:28:59 +02:00
y ( y_index_eq , it_ ) = phat ;
2011-01-14 19:22:29 +01:00
else
y ( y_index_eq ) = phat ;
2015-02-20 14:44:45 +01:00
end
if is_dynamic
2020-06-19 18:34:48 +02:00
[ r , ~ , T ( : , it_ ) , g1 ] = feval ( fname , Block_Num , dynvars_from_endo_simul ( y , it_ , M ) , x , params , ...
2020-06-17 16:48:40 +02:00
steady_state , T ( : , it_ ) , it_ , false ) ;
2011-01-14 19:22:29 +01:00
else
2020-06-19 18:34:48 +02:00
[ r , ~ , T , g1 ] = feval ( fname , Block_Num , y , x , params , T ) ;
2015-02-20 14:44:45 +01:00
end
if max ( abs ( r ) ) > = options . solve_tolf
2011-01-14 14:45:48 +01:00
[ dx , flag1 ] = bicgstab ( g1 , - r , 1e-7 , Blck_size , L1 , U1 ) ;
2011-01-14 19:22:29 +01:00
else
2011-02-04 17:17:48 +01:00
flag1 = 0 ;
dx = phat - ya ;
2015-02-20 14:44:45 +01:00
end
if flag1 > 0 || reduced
if verbose
if ( flag1 == 1 )
disp ( [ ' Error in simul: No convergence inside BICGSTAB after ' num2str ( iter , ' %6d' ) ' iterations, in block' num2str ( Block_Num , ' %3d' ) ] )
elseif ( flag1 == 2 )
disp ( [ ' Error in simul: Preconditioner is ill-conditioned, in block' num2str ( Block_Num , ' %3d' ) ] )
elseif ( flag1 == 3 )
disp ( [ ' Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str ( Block_Num , ' %3d' ) ] )
end
end
2013-03-25 12:03:09 +01:00
ilu_setup . droptol = ilu_setup . droptol / 10 ;
2009-12-16 18:17:34 +01:00
reduced = 0 ;
else
ya = ya + lambda * dx ;
2015-02-20 14:44:45 +01:00
if is_dynamic
2020-06-17 17:28:59 +02:00
y ( y_index_eq , it_ ) = ya ;
2009-12-16 18:17:34 +01:00
else
y ( y_index_eq ) = ya ' ;
2015-02-20 14:44:45 +01:00
end
end
end
2009-12-16 18:17:34 +01:00
else
2015-02-20 14:44:45 +01:00
if verbose
disp ( ' unknown option : ' )
if is_dynamic
disp ( [ ' options_.stack_solve_algo = ' num2str ( stack_solve_algo ) ' not implemented' ] )
else
disp ( [ ' options_.solve_algo = ' num2str ( options . solve_algo ) ' not implemented' ] )
end
end
2009-12-16 18:17:34 +01:00
info = - Block_Num * 10 ;
2015-02-20 14:44:45 +01:00
return
end
2009-12-16 18:17:34 +01:00
iter = iter + 1 ;
max_resa = max_res ;
end
end
if cvg == 0
2015-02-20 14:44:45 +01:00
if verbose
if cutoff == 0
fprintf ( ' Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit\".\n' , Block_Num , it_ , iter ) ;
else
fprintf ( ' Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n' , Block_Num , it_ , iter ) ;
end
end
2020-05-25 16:58:11 +02:00
if is_dynamic
2020-06-02 17:45:28 +02:00
oo_ . deterministic_simulation . status = false ;
2009-12-16 18:17:34 +01:00
oo_ . deterministic_simulation . error = max_res ;
oo_ . deterministic_simulation . iterations = iter ;
2020-06-02 17:45:28 +02:00
oo_ . deterministic_simulation . block ( Block_Num ) . status = false ; % Convergency failed.
2009-12-16 18:17:34 +01:00
oo_ . deterministic_simulation . block ( Block_Num ) . error = max_res ;
oo_ . deterministic_simulation . block ( Block_Num ) . iterations = iter ;
2017-05-16 12:42:01 +02:00
end
2009-12-16 18:17:34 +01:00
info = - Block_Num * 10 ;
2015-02-20 14:44:45 +01:00
return
2009-12-16 18:17:34 +01:00
end
end
2015-02-20 14:44:45 +01:00
if is_dynamic
2010-10-22 16:49:47 +02:00
info = 1 ;
2020-06-02 17:45:28 +02:00
oo_ . deterministic_simulation . status = true ;
2009-12-16 18:17:34 +01:00
oo_ . deterministic_simulation . error = max_res ;
oo_ . deterministic_simulation . iterations = iter ;
2020-06-02 17:45:28 +02:00
oo_ . deterministic_simulation . block ( Block_Num ) . status = true ;
2009-12-16 18:17:34 +01:00
oo_ . deterministic_simulation . block ( Block_Num ) . error = max_res ;
oo_ . deterministic_simulation . block ( Block_Num ) . iterations = iter ;
2010-10-22 16:49:47 +02:00
else
2011-02-04 17:17:48 +01:00
info = 0 ;
2015-02-20 14:44:45 +01:00
end