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;