2023-06-06 18:13:10 +02:00
function [ y , T , success , max_res , iter ] = solve_one_boundary ( fh , y , x , params , steady_state , T , ...
y_index_eq , nze , periods , is_linear , Block_Num , y_kmin , maxit_ , solve_tolf , cutoff , stack_solve_algo , is_forward , is_dynamic , verbose , M , options )
2023-01-10 18:01:27 +01:00
% Computes the deterministic simulation or the steady state for a block of equations containing
% only lags or only leads (but not both).
2009-12-16 18:17:34 +01:00
%
% INPUTS
2023-01-10 18:01:27 +01:00
% fh [handle] function handle to the static/dynamic file for the block
2009-12-16 18:17:34 +01:00
% 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
% cutoff [double] cutoff to correct the direction in Newton in case
% of singular jacobian matrix
2022-06-15 14:46:44 +02:00
% stack_solve_algo [integer] linear solver method used in the Newton algorithm
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
2023-06-06 18:13:10 +02:00
% is_dynamic [logical] Whether this is a deterministic simulation
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
2023-06-06 18:13:10 +02:00
% y [matrix] All endogenous variables of the model
% T [matrix] Temporary terms
% success [logical] Whether a solution was found
% max_res [double] ∞-norm of the residual
% iter [integer] Number of iterations
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
2023-01-10 18:01:27 +01:00
% Copyright © 1996-2023 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 ;
2022-03-31 17:39:43 +02:00
lambda = 1 ; % Length of Newton step
2009-12-16 18:17:34 +01:00
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
2022-03-23 16:49:05 +01:00
cvg = false ;
2009-12-16 18:17:34 +01:00
iter = 0 ;
g1 = spalloc ( Blck_size , Blck_size , nze ) ;
2022-03-23 16:49:05 +01:00
while ~ ( cvg || iter > maxit_ )
2015-02-20 14:44:45 +01:00
if is_dynamic
2023-01-10 18:01:27 +01:00
[ yy , T ( : , it_ ) , r , g1 ] = fh ( dynendo ( y , it_ , M ) , x ( it_ , : ) , params , steady_state , ...
M . block_structure . block ( Block_Num ) . g1_sparse_rowval , ...
M . block_structure . block ( Block_Num ) . g1_sparse_colval , ...
M . block_structure . block ( Block_Num ) . g1_sparse_colptr , T ( : , it_ ) ) ;
y ( : , it_ ) = yy ( M . endo_nbr + ( 1 : M . endo_nbr ) ) ;
2009-12-16 18:17:34 +01:00
else
2023-01-10 18:01:27 +01:00
[ y , T , r , g1 ] = fh ( y , x , params , M . block_structure_stat . block ( Block_Num ) . g1_sparse_rowval , ...
M . block_structure_stat . block ( Block_Num ) . g1_sparse_colval , ...
M . block_structure_stat . block ( Block_Num ) . g1_sparse_colptr , 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
2022-03-16 13:00:50 +01:00
disp ( [ char ( M . endo_names { y_index_eq } ) repmat ( ' ' , numel ( y_index_eq ) , 1 ) num2str ( [ y ( y_index_eq , it_ ) r g1 ] ) ] )
2009-12-16 18:17:34 +01:00
else
2022-03-16 13:00:50 +01:00
disp ( [ char ( M . endo_names { y_index_eq } ) repmat ( ' ' , numel ( y_index_eq ) , 1 ) num2str ( [ y ( y_index_eq ) r g1 ] ) ] )
2015-02-20 14:44:45 +01:00
end
end
if ~ isreal ( max_res ) || isnan ( max_res )
2022-03-23 16:49:05 +01:00
cvg = false ;
2015-02-20 14:44:45 +01:00
elseif is_linear && iter > 0
2022-03-23 16:49:05 +01:00
cvg = true ;
2009-12-16 18:17:34 +01:00
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
2022-03-23 16:50:16 +01:00
disp ( [ ' The Jacobian matrix is singular, det(Jacobian)=' num2str ( detJ , ' %f' ) ' .' ] )
2015-02-20 14:44:45 +01:00
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
2023-06-06 18:13:10 +02:00
success = false ;
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
2023-06-15 15:49:17 +02:00
fprintf ( ' Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "maxit".\n' , Block_Num , it_ , iter ) ;
2015-02-20 14:44:45 +01:00
else
2023-06-15 15:49:17 +02:00
fprintf ( ' Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "maxit" or set "cutoff=0" in model options.\n' , Block_Num , it_ , iter ) ;
2015-02-20 14:44:45 +01:00
end
end
2023-06-06 18:13:10 +02:00
success = false ;
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 , ...
2023-01-10 18:01:27 +01:00
@ lnsrch1_wrapper_one_boundary , nn , ...
nn , options . solve_tolx , y_index_eq , fh , Block_Num , y , x , params , steady_state , T ( : , it_ ) , it_ , M ) ;
2022-02-11 16:47:58 +01:00
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
2023-01-10 18:01:27 +01:00
[ ~ , T ( : , it_ ) ] = fh ( dynendo ( y , it_ , M ) , x ( it_ , : ) , params , steady_state , ...
M . block_structure . block ( Block_Num ) . g1_sparse_rowval , ...
M . block_structure . block ( Block_Num ) . g1_sparse_colval , ...
M . block_structure . block ( Block_Num ) . g1_sparse_colptr , T ( : , it_ ) ) ;
2022-06-15 14:46:44 +02:00
elseif ( is_dynamic && ( stack_solve_algo == 1 || stack_solve_algo == 0 || stack_solve_algo == 6 ) ) || ( ~ 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 ) ;
2022-03-16 13:02:40 +01:00
[ dx , flag1 ] = bicgstab ( g1 , - r , 1e-6 , Blck_size , L1 , U1 ) ;
2015-02-20 14:44:45 +01:00
if flag1 > 0 || reduced
if verbose
if ( flag1 == 1 )
2022-03-16 13:02:40 +01:00
disp ( [ ' Error in simul: No convergence inside BiCGStab after ' num2str ( iter , ' %6d' ) ' iterations, in block' num2str ( Block_Num , ' %3d' ) ] )
2015-02-20 14:44:45 +01:00
elseif ( flag1 == 2 )
disp ( [ ' Error in simul: Preconditioner is ill-conditioned, in block' num2str ( Block_Num , ' %3d' ) ] )
elseif ( flag1 == 3 )
2022-03-16 13:02:40 +01:00
disp ( [ ' Error in simul: BiCGStab stagnated (Two consecutive iterates were the same.), in block' num2str ( Block_Num , ' %3d' ) ] )
2015-02-20 14:44:45 +01:00
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
2022-03-23 17:33:05 +01:00
if is_dynamic
error ( [ ' options_.stack_solve_algo = ' num2str ( stack_solve_algo ) ' not implemented' ] )
else
error ( [ ' options_.solve_algo = ' num2str ( options . solve_algo ) ' not implemented' ] )
2015-02-20 14:44:45 +01:00
end
end
2009-12-16 18:17:34 +01:00
iter = iter + 1 ;
max_resa = max_res ;
end
end
2022-03-23 16:49:05 +01:00
if ~ cvg
2015-02-20 14:44:45 +01:00
if verbose
if cutoff == 0
2023-06-15 15:49:17 +02:00
fprintf ( ' Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "maxit".\n' , Block_Num , it_ , iter ) ;
2015-02-20 14:44:45 +01:00
else
2023-06-15 15:49:17 +02:00
fprintf ( ' Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "maxit" or set "cutoff=0" in model options.\n' , Block_Num , it_ , iter ) ;
2015-02-20 14:44:45 +01:00
end
end
2023-06-06 18:13:10 +02:00
success = false ;
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
2023-06-06 18:13:10 +02:00
success = true ;
2023-01-10 18:01:27 +01:00
function y3n = dynendo ( y, it_, M)
if it_ > 1 && it_ < size ( y , 2 )
y3n = reshape ( y ( : , it_ + ( - 1 : 1 ) ) , 3 * M . endo_nbr , 1 ) ;
elseif it_ > 1 % Purely backward model (in last period)
y3n = [ reshape ( y ( : , it_ + ( - 1 : 0 ) ) , 2 * M . endo_nbr , 1 ) ; NaN ( M_ . endo_nbr , 1 ) ] ;
elseif it_ < size ( y , 2 ) % Purely forward model (in first period)
y3n = [ NaN ( M_ . endo_nbr , 1 ) ; reshape ( y ( : , it_ + ( 0 : 1 ) ) , 2 * M . endo_nbr , 1 ) ] ;
else % Static model
y3n = [ NaN ( M_ . endo_nbr , 1 ) ; y ( : , it_ ) ; NaN ( M_ . endo_nbr , 1 ) ]
end
function r = lnsrch1_wrapper_one_boundary ( ya, y_index, fh, Block_Num, y, x, params, steady_state, T, it_, M)
y ( y_index , it_ ) = ya ;
[ ~ , ~ , r ] = fh ( dynendo ( y , it_ , M ) , x ( it_ , : ) , params , steady_state , ...
M . block_structure . block ( Block_Num ) . g1_sparse_rowval , ...
M . block_structure . block ( Block_Num ) . g1_sparse_colval , ...
M . block_structure . block ( Block_Num ) . g1_sparse_colptr , T ) ;