1function [y, oo]= solve_two_boundaries(fname, y, x, params, steady_state, y_index, nze, periods, y_kmin_l, y_kmax_l, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo,options,M, oo) 2% Computes the deterministic simulation of a block of equation containing 3% both lead and lag variables using relaxation methods 4% 5% INPUTS 6% fname [string] name of the file containing the block 7% to simulate 8% y [matrix] All the endogenous variables of the model 9% x [matrix] All the exogenous variables of the model 10% params [vector] All the parameters of the model 11% steady_state [vector] steady state of the model 12% y_index [vector of int] The index of the endogenous variables of 13% the block 14% nze [integer] number of non-zero elements in the 15% jacobian matrix 16% periods [integer] number of simulation periods 17% y_kmin_l [integer] maximum number of lag in the block 18% y_kmax_l [integer] maximum number of lead in the block 19% is_linear [integer] if is_linear=1 the block is linear 20% if is_linear=0 the block is not linear 21% Block_Num [integer] block number 22% y_kmin [integer] maximum number of lag in the model 23% maxit_ [integer] maximum number of iteration in Newton 24% solve_tolf [double] convergence criteria 25% lambda [double] initial value of step size in 26% Newton 27% cutoff [double] cutoff to correct the direction in Newton in case 28% of singular jacobian matrix 29% stack_solve_algo [integer] linear solver method used in the 30% Newton algorithm : 31% - 1 sprse LU 32% - 2 GMRES 33% - 3 BicGStab 34% - 4 Optimal path length 35% M [structure] Model description 36% oo [structure] Results 37% 38% OUTPUTS 39% y [matrix] All endogenous variables of the model 40% oo [structure] Results 41% 42% ALGORITHM 43% Newton with LU or GMRES or BicGstab 44% 45% SPECIAL REQUIREMENTS 46% none. 47% 48 49% Copyright (C) 1996-2018 Dynare Team 50% 51% This file is part of Dynare. 52% 53% Dynare is free software: you can redistribute it and/or modify 54% it under the terms of the GNU General Public License as published by 55% the Free Software Foundation, either version 3 of the License, or 56% (at your option) any later version. 57% 58% Dynare is distributed in the hope that it will be useful, 59% but WITHOUT ANY WARRANTY; without even the implied warranty of 60% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 61% GNU General Public License for more details. 62% 63% You should have received a copy of the GNU General Public License 64% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 65 66verbose = options.verbosity; 67 68cvg=0; 69iter=0; 70Per_u_=0; 71g2 = []; 72g3 = []; 73Blck_size=size(y_index,2); 74correcting_factor=0.01; 75ilu_setup.droptol=1e-10; 76ilu_setup.type = 'ilutp'; 77%ilu_setup.milu = 'col'; 78ilu_setup.milu = 'off'; 79ilu_setup.thresh = 1; 80ilu_setup.udiag = 0; 81max_resa=1e100; 82Jacobian_Size=Blck_size*(y_kmin+y_kmax_l +periods); 83g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods); 84reduced = 0; 85while ~(cvg==1 || iter>maxit_) 86 [r, y, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, periods, 0, y_kmin, Blck_size,options.periods); 87 preconditioner = 2; 88 g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size); 89 term1 = g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'; 90 term2 = g1(:, (periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)'; 91 b = b - term1 - term2; 92 [max_res, max_indx]=max(max(abs(r'))); 93 if ~isreal(r) 94 max_res = (-max_res^2)^0.5; 95 end 96 if ~isreal(max_res) || isnan(max_res) 97 cvg = 0; 98 elseif(is_linear && iter>0) 99 cvg = 1; 100 else 101 cvg=(max_res<solve_tolf); 102 end 103 if ~cvg 104 if iter>0 105 if ~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1) 106 if verbose && ~isreal(max_res) 107 disp(['Variable ' M.endo_names{max_indx} ' (' int2str(max_indx) ') returns an undefined value']); 108 end 109 if isnan(max_res) 110 detJ=det(g1aa); 111 if abs(detJ)<1e-7 112 max_factor=max(max(abs(g1aa))); 113 ze_elem=sum(diag(g1aa)<cutoff); 114 if verbose 115 disp([num2str(full(ze_elem),'%d') ' elements on the Jacobian diagonal are below the cutoff (' num2str(cutoff,'%f') ')']); 116 end 117 if correcting_factor<max_factor 118 correcting_factor=correcting_factor*4; 119 if verbose 120 disp(['The Jacobain matrix is singular, det(Jacobian)=' num2str(detJ,'%f') '.']); 121 disp([' trying to correct the Jacobian matrix:']); 122 disp([' correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')]); 123 end 124 dx = (g1aa+correcting_factor*speye(periods*Blck_size))\ba- ya; 125 y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)'; 126 continue 127 else 128 disp('The singularity of the jacobian matrix could not be corrected'); 129 return 130 end 131 end 132 elseif lambda>1e-8 133 lambda=lambda/2; 134 reduced = 1; 135 if verbose 136 disp(['reducing the path length: lambda=' num2str(lambda,'%f')]); 137 end 138 y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)'; 139 continue 140 else 141 if verbose 142 if cutoff==0 143 fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.simul.maxit".\n',Block_Num, iter); 144 else 145 fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, iter); 146 end 147 end 148 oo.deterministic_simulation.status = 0; 149 oo.deterministic_simulation.error = max_res; 150 oo.deterministic_simulation.iterations = iter; 151 oo.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed. 152 oo.deterministic_simulation.block(Block_Num).error = max_res; 153 oo.deterministic_simulation.block(Block_Num).iterations = iter; 154 return 155 end 156 else 157 if lambda<1 158 lambda=max(lambda*2, 1); 159 end 160 end 161 end 162 ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)'; 163 ya_save=ya; 164 g1aa=g1a; 165 ba=b; 166 max_resa=max_res; 167 if stack_solve_algo==0 168 dx = g1a\b- ya; 169 ya = ya + lambda*dx; 170 y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; 171 elseif stack_solve_algo==1 172 for t=1:periods 173 first_elem = (t-1)*Blck_size+1; 174 last_elem = t*Blck_size; 175 next_elem = (t+1)*Blck_size; 176 Elem = first_elem:last_elem; 177 Elem_1 = last_elem+1:next_elem; 178 B1_inv = inv(g1a(Elem, Elem)); 179 if (t < periods) 180 S1 = B1_inv * g1a(Elem, Elem_1); 181 end 182 g1a(Elem, Elem_1) = S1; 183 b(Elem) = B1_inv * b(Elem); 184 g1a(Elem, Elem) = ones(Blck_size, Blck_size); 185 if t<periods 186 g1a(Elem_1, Elem_1) = g1a(Elem_1, Elem_1) - g1a(Elem_1, Elem) * S1; 187 b(Elem_1) = b(Elem_1) - g1a(Elem_1, Elem) * b(Elem); 188 g1a(Elem_1, Elem) = zeros(Blck_size, Blck_size); 189 end 190 end 191 za = b(Elem); 192 zaa = za; 193 y_Elem = (periods - 1) * Blck_size + 1:(periods) * Blck_size; 194 dx = ya; 195 dx(y_Elem) = za - ya(y_Elem); 196 ya(y_Elem) = ya(y_Elem) + lambda*dx(y_Elem); 197 for t=periods-1:-1:1 198 first_elem = (t-1)*Blck_size+1; 199 last_elem = t*Blck_size; 200 next_elem = (t+1)*Blck_size; 201 Elem_1 = last_elem+1:next_elem; 202 Elem = first_elem:last_elem; 203 za = b(Elem) - g1a(Elem, Elem_1) * zaa; 204 zaa = za; 205 y_Elem = Blck_size * (t-1)+1:Blck_size * (t); 206 dx(y_Elem) = za - ya(y_Elem); 207 ya(y_Elem) = ya(y_Elem) + lambda*dx(y_Elem); 208 y(y_kmin + t, y_index) = ya(y_Elem); 209 end 210 elseif stack_solve_algo==2 211 flag1=1; 212 while flag1>0 213 if preconditioner==2 214 [L1, U1]=ilu(g1a,ilu_setup); 215 elseif preconditioner==3 216 Size = Blck_size; 217 gss1 = g1a(Size + 1: 2*Size,Size + 1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size); 218 [L1, U1]=lu(gss1); 219 L(1:Size,1:Size) = L1; 220 U(1:Size,1:Size) = U1; 221 gss2 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size); 222 [L2, U2]=lu(gss2); 223 L(Size+1:(periods-1)*Size,Size+1:(periods-1)*Size) = kron(eye(periods-2), L2); 224 U(Size+1:(periods-1)*Size,Size+1:(periods-1)*Size) = kron(eye(periods-2), U2); 225 gss2 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size); 226 [L3, U3]=lu(gss2); 227 L((periods-1)*Size+1:periods*Size,(periods-1)*Size+1:periods*Size) = L3; 228 U((periods-1)*Size+1:periods*Size,(periods-1)*Size+1:periods*Size) = U3; 229 L1 = L; 230 U1 = U; 231 elseif preconditioner==4 232 Size = Blck_size; 233 gss1 = g1a(1: 3*Size, 1: 3*Size); 234 [L, U] = lu(gss1); 235 L1 = kron(eye(ceil(periods/3)),L); 236 U1 = kron(eye(ceil(periods/3)),U); 237 L1 = L1(1:periods * Size, 1:periods * Size); 238 U1 = U1(1:periods * Size, 1:periods * Size); 239 end 240 [za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1); 241 if (flag1>0 || reduced) 242 if verbose 243 if flag1==1 244 disp(['Error in simul: No convergence inside GMRES after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]); 245 elseif flag1==2 246 disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]); 247 elseif flag1==3 248 disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]); 249 end 250 end 251 ilu_setup.droptol = ilu_setup.droptol/10; 252 reduced = 0; 253 else 254 dx = za - ya; 255 ya = ya + lambda*dx; 256 y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; 257 end 258 end 259 elseif stack_solve_algo==3 260 flag1=1; 261 while flag1>0 262 if preconditioner==2 263 [L1, U1]=ilu(g1a,ilu_setup); 264 [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1); 265 elseif preconditioner==3 266 Size = Blck_size; 267 gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size); 268 [L1, U1]=lu(gss0); 269 P1 = eye(size(gss0)); 270 Q1 = eye(size(gss0)); 271 L = kron(eye(periods),L1); 272 U = kron(eye(periods),U1); 273 P = kron(eye(periods),P1); 274 Q = kron(eye(periods),Q1); 275 [za,flag1] = bicgstab1(g1a,b,1e-7,Blck_size*periods,L,U, P, Q); 276 else 277 Size = Blck_size; 278 gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size); 279 [L1, U1]=lu(gss0); 280 L1 = kron(eye(periods),L1); 281 U1 = kron(eye(periods),U1); 282 [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1); 283 end 284 if flag1>0 || reduced 285 if verbose 286 if flag1==1 287 disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]); 288 elseif flag1==2 289 disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]); 290 elseif flag1==3 291 disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]); 292 end 293 end 294 ilu_setup.droptol = ilu_setup.droptol/10; 295 reduced = 0; 296 else 297 dx = za - ya; 298 ya = ya + lambda*dx; 299 y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; 300 end 301 end 302 elseif stack_solve_algo==4 303 ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*Blck_size, 1); 304 stpmx = 100 ; 305 stpmax = stpmx*max([sqrt(ya'*ya);size(y_index,2)]); 306 nn=1:size(ra,1); 307 g = (ra'*g1a)'; 308 f = 0.5*ra'*ra; 309 p = -g1a\ra; 310 [yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,'lnsrch1_wrapper_two_boundaries',nn,nn, options.solve_tolx, fname, y, y_index,x, params, steady_state, periods, y_kmin, Blck_size,options.periods); 311 dx = ya - yn; 312 y(1+y_kmin:periods+y_kmin,y_index)=reshape(yn',length(y_index),periods)'; 313 end 314 end 315 iter=iter+1; 316 if verbose 317 disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e')]); 318 end 319end 320 321if (iter>maxit_) 322 if verbose 323 printline(41) 324 %disp(['No convergence after ' num2str(iter,'%4d') ' iterations in Block ' num2str(Block_Num,'%d')]) 325 end 326 oo.deterministic_simulation.status = 0; 327 oo.deterministic_simulation.error = max_res; 328 oo.deterministic_simulation.iterations = iter; 329 oo.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed. 330 oo.deterministic_simulation.block(Block_Num).error = max_res; 331 oo.deterministic_simulation.block(Block_Num).iterations = iter; 332 return 333end 334 335oo.deterministic_simulation.status = 1; 336oo.deterministic_simulation.error = max_res; 337oo.deterministic_simulation.iterations = iter; 338oo.deterministic_simulation.block(Block_Num).status = 1;% Convergency obtained. 339oo.deterministic_simulation.block(Block_Num).error = max_res; 340oo.deterministic_simulation.block(Block_Num).iterations = iter;