solve_one_boundary.m 20.3 KB
Newer Older
1
function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, ...
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
2
                                        y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo, forward_backward, is_dynamic, verbose, M, options, oo)
3
4
5
6
7
8
9
10
11
% Computes the deterministic simulation of a block of equation containing
% lead or lag variables 
%
% 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
12
%   steady_state        [vector]        steady state of the model
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
%   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
%   is_linear           [integer]       if is_linear=1 the block is linear
%                                       if is_linear=0 the block is not linear
%   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
%                                       Newton algorithm : 
%                                            - 1 sparse LU
%                                            - 2 GMRES
%                                            - 3 BicGStab
%                                            - 4 Optimal path length
%   forward_backward    [integer]       The block has to be solve forward
%                                       (1) or backward (0)
%   is_dynamic          [integer]       (1) The block belong to the dynamic
%                                           file and the oo_.deterministic_simulation field has to be uptated
%                                       (0) The block belong to the static
%                                           file and th oo_.deteerminstic
%                                           field remains unchanged
%   verbose            [integer]        (0) iterations are not printed
%                                       (1) iterations are printed
43
44
45
46
%   indirect_call      [integer]        (0) direct call to the fname
%                                       (1) indirect call via the
%                                       local_fname wrapper
% OUTPUTS                                    
47
%   y                  [matrix]         All endogenous variables of the model      
48
49
%   info               [integer]        >=0 no error
%                                       <0 error
50
51
52
53
54
55
56
57
%  
% ALGORITHM
%   Newton with LU or GMRES or BicGstab for dynamic block
%    
% SPECIAL REQUIREMENTS
%   none.
%  

58
% Copyright (C) 1996-2013 Dynare Team
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
%
% 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/>.


Blck_size=size(y_index_eq,2);
g2 = [];
g3 = [];
correcting_factor=0.01;
80
ilu_setup.droptol=1e-10;
81
82
83
84
85
86
87
88
89
90
91
max_resa=1e100;
reduced = 0;
if(forward_backward)
    incr = 1;
    start = y_kmin+1;
    finish = periods+y_kmin;
else
    incr = -1;
    start = periods+y_kmin;
    finish = y_kmin+1;
end
92
%lambda=1;
93
94
95
96
for it_=start:incr:finish
    cvg=0;
    iter=0;
    g1=spalloc( Blck_size, Blck_size, nze);
97
    while ~(cvg==1 || iter>maxit_),
98
        if(is_dynamic)
99
100
            [r, y, g1, g2, g3] = feval(fname, y, x, params, steady_state, ...
                                       it_, 0);
101
        else
102
            [r, y, g1] = feval(fname, y, x, params);
103
104
105
106
107
108
109
110
111
112
113
114
115
116
        end;
        if(~isreal(r))
            max_res=(-(max(max(abs(r))))^2)^0.5;
        else
            max_res=max(max(abs(r)));
        end;
        %['max_res=' num2str(max_res) ' Block_Num=' int2str(Block_Num) ' it_=' int2str(it_)]
        %disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)]);
        
        %     fjac = zeros(Blck_size, Blck_size);
        %     disp(['Blck_size=' int2str(Blck_size) ' it_=' int2str(it_)]);
        %     dh = max(abs(y(it_, y_index_eq)),options_.gstep*ones(1, Blck_size))*eps^(1/3);
        %     fvec = r;
        %       for j = 1:Blck_size
117
        %         ydh = y ;
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
        %           ydh(it_, y_index_eq(j)) = ydh(it_, y_index_eq(j)) + dh(j)  ;
        %           [t, y1, g11, g21, g31]=feval(fname, ydh, x, params, it_, 0);
        %           fjac(:,j) = (t - fvec)./dh(j) ;
        %       end;
        %     diff = g1 -fjac;
        %     diff
        %     disp('g1');
        %     disp([num2str(g1,'%4.5f')]);
        %     disp('fjac');
        %     disp([num2str(fjac,'%4.5f')]);
        %     [c_max, i_c_max] = max(abs(diff));
        %     [l_c_max, i_r_max] = max(c_max);
        %     disp(['maximum element row=' int2str(i_c_max(i_r_max)) ' and column=' int2str(i_r_max) ' value = ' num2str(l_c_max)]);
        %     equation = i_c_max(i_r_max);
        %     variable = i_r_max;
        %     variable
        %     mod(variable, Blck_size)
        %     disp(['equation ' int2str(equation) ' and variable ' int2str(y_index_eq(variable)) ' ' M_.endo_names(y_index_eq(variable), :)]);
        %     disp(['g1(' int2str(equation) ', ' int2str(variable) ')=' num2str(g1(equation, variable),'%3.10f') ' fjac(' int2str(equation) ', ' int2str(variable) ')=' num2str(fjac(equation, variable), '%3.10f') ' y(' int2str(it_) ', ' int2str(variable) ')=' num2str(y(it_, variable))]);
        %     %return;
        %     %g1 = fjac;

        
        
        
        
        if(verbose==1)
            disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)]);
            if(is_dynamic)
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
147
                disp([M.endo_names(y_index_eq,:) num2str([y(it_,y_index_eq)' r g1])]);
148
            else
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
149
                disp([M.endo_names(y_index_eq,:) num2str([y(y_index_eq) r g1])]);
150
151
            end;
        end;
152
        if(~isreal(max_res) || isnan(max_res))
153
            cvg = 0;
154
        elseif(is_linear && iter>0)
155
156
157
158
159
160
            cvg = 1;
        else
            cvg=(max_res<solve_tolf);
        end;
        if(~cvg)
            if(iter>0)
161
162
                if(~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1))
                    if(isnan(max_res) || (max_resa<max_res && iter>0))
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
                        detJ=det(g1a);
                        if(abs(detJ)<1e-7)
                            max_factor=max(max(abs(g1a)));
                            ze_elem=sum(diag(g1a)<cutoff);
                            disp([num2str(full(ze_elem),'%d') ' elements on the Jacobian diagonal are below the cutoff (' num2str(cutoff,'%f') ')']);
                            if(correcting_factor<max_factor)
                                correcting_factor=correcting_factor*4;
                                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')]);
                                dx = - r/(g1+correcting_factor*speye(Blck_size));
                                %dx = -b'/(g1+correcting_factor*speye(Blck_size))-ya_save;
                                y(it_,y_index_eq)=ya_save+lambda*dx;
                                continue;
                            else
                                disp('The singularity of the jacobian matrix could not be corrected');
                                info = -Block_Num*10;
                                return;
                            end;
                        end;
                    elseif(lambda>1e-8)
                        lambda=lambda/2;
                        reduced = 1;
                        disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);
                        if(is_dynamic)
                            y(it_,y_index_eq)=ya_save-lambda*dx;
                        else
                            y(y_index_eq)=ya_save-lambda*dx;
                        end;
                        continue;
                    else
                        if(cutoff == 0)
195
                            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);
196
                        else
197
                            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);
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
                        end;
                        if(is_dynamic)
                            oo_.deterministic_simulation.status = 0;
                            oo_.deterministic_simulation.error = max_res;
                            oo_.deterministic_simulation.iterations = iter;
                            oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed.
                            oo_.deterministic_simulation.block(Block_Num).error = max_res;
                            oo_.deterministic_simulation.block(Block_Num).iterations = iter;
                        end;
                        info = -Block_Num*10;
                        return;
                    end;
                else
                    if(lambda<1)
                        lambda=max(lambda*2, 1);
                    end;
                end;
            end;
            if(is_dynamic)
                ya = y(it_,y_index_eq)';
            else
                ya = y(y_index_eq);
            end;
            ya_save=ya;
            g1a=g1;
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
223
            if(~is_dynamic && options.solve_algo == 0)
224
225
226
                if (verbose == 1)
                    disp('steady: fsolve');
                end
227
                if ~exist('OCTAVE_VERSION')
228
                    if ~user_has_matlab_license('optimization_toolbox')
229
230
                        error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
                    end
231
232
233
234
235
236
237
                end
                options=optimset('fsolve');
                options.MaxFunEvals = 50000;
                options.MaxIter = 2000;
                options.TolFun=1e-8;
                options.Display = 'iter';
                options.Jacobian = 'on';
238
239
240
241
242
243
244
245
                if ~exist('OCTAVE_VERSION')
                    [yn,fval,exitval,output] = fsolve(@local_fname, y(y_index_eq), ...
                                                      options, x, params, steady_state, y, y_index_eq, fname, 0);
                else
                    % Under Octave, use a wrapper, since fsolve() does not have a 4th arg
                    func = @(z) local_fname(z, x, params, steady_state, y, y_index_eq, fname, 0);
                    % The Octave version of fsolve does not converge when it starts from the solution
                    fvec = feval(func,y(y_index_eq));
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
246
                    if max(abs(fvec)) >= options.solve_tolf
247
248
249
250
251
252
253
                        [yn,fval,exitval,output] = fsolve(func,y(y_index_eq),options);
                    else
                        yn = y(y_index_eq);
                        exitval = 3;
                    end;
                end
                    
254
255
256
257
258
259
                y(y_index_eq) = yn;
                if exitval > 0
                    info = 0;
                else
                    info = -Block_Num*10;
                end
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
260
            elseif((~is_dynamic && options.solve_algo==2) || (is_dynamic && stack_solve_algo==4))
261
                if (verbose == 1 && ~is_dynamic)
262
263
                    disp('steady: LU + lnsrch1');
                end
264
265
266
                lambda=1;
                stpmx = 100 ;
                if (is_dynamic)
sebastien's avatar
sebastien committed
267
                    stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
268
                else
sebastien's avatar
sebastien committed
269
                    stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
270
271
272
273
274
275
                end;
                nn=1:size(y_index_eq,2);
                g = (r'*g1)';
                f = 0.5*r'*r;
                p = -g1\r ;
                if (is_dynamic)
276
277
278
                    [ya,f,r,check]=lnsrch1(y(it_,:)',f,g,p,stpmax, ...
                                           'lnsrch1_wrapper_one_boundary',nn, ...
                                           y_index_eq, y_index_eq, fname, y, x, params, steady_state, it_);
279
                    dx = ya' - y(it_, :);
280
                else
281
282
                    [ya,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, ...
                                           params, steady_state,0);
283
                    dx = ya - y(y_index_eq);
284
                end;
285
                
sebastien's avatar
sebastien committed
286
287
288
289
290
                if(is_dynamic)
                    y(it_,:) = ya';
                else
                    y = ya';
                end;
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
291
            elseif(~is_dynamic && options.solve_algo==3)
292
293
294
                if (verbose == 1)
                    disp('steady: csolve');
                end
295
296
                [yn,info] = csolve(@local_fname, y(y_index_eq),@ ...
                                   local_fname,1e-6,500, x, params, steady_state, y, y_index_eq, fname, 1);
297
298
                dx = ya - yn;
                y(y_index_eq) = yn;
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
299
            elseif((stack_solve_algo==1 && is_dynamic) || (stack_solve_algo==0 && is_dynamic) || (~is_dynamic && (options.solve_algo==1 || options.solve_algo==6))),
300
                if (verbose == 1 && ~is_dynamic)
301
302
                    disp('steady: Sparse LU ');
                end
303
304
305
306
307
308
309
                dx =  g1\r;
                ya = ya - lambda*dx;
                if(is_dynamic)
                    y(it_,y_index_eq) = ya';
                else
                    y(y_index_eq) = ya;
                end;
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
310
            elseif((stack_solve_algo==2 && is_dynamic) || (options.solve_algo==7 && ~is_dynamic)),
311
                flag1=1;
312
                if exist('OCTAVE_VERSION')
313
314
                    error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=7 since GMRES is not implemented in Octave')
                end
315
                if (verbose == 1 && ~is_dynamic)
316
317
                    disp('steady: GMRES ');
                end
318
                while(flag1>0)
319
                    [L1, U1]=ilu(g1,ilu_setup);
sebastien's avatar
sebastien committed
320
                    [dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
321
                    if (flag1>0 || reduced)
322
323
324
325
326
327
328
                        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;
329
                        ilu_setup.droptol = ilu_setup.droptol/10;
330
331
332
333
334
335
336
337
338
339
                        reduced = 0;
                    else
                        ya = ya + lambda*dx;
                        if(is_dynamic)
                            y(it_,y_index_eq) = ya';
                        else
                            y(y_index_eq) = ya';
                        end;
                    end;
                end;
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
340
            elseif((stack_solve_algo==3 && is_dynamic) || (options.solve_algo==8 && ~is_dynamic)),
341
                flag1=1;
342
                if (verbose == 1 && ~is_dynamic)
343
344
                    disp('steady: BiCGStab');
                end
345
                while(flag1>0)
346
                    [L1, U1]=ilu(g1,ilu_setup);
347
                    phat = ya - U1 \ (L1 \ r);
348
349
350
351
352
                    if(is_dynamic)
                        y(it_,y_index_eq) = phat;
                    else
                        y(y_index_eq) = phat;
                    end;
353
                    if(is_dynamic)
354
355
                        [r, y, g1, g2, g3] = feval(fname, y, x, params, ...
                                                   steady_state, it_, 0);
356
                    else
357
                        [r, y, g1] = feval(fname, y, x, params);
358
                    end;
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
359
                    if max(abs(r)) >= options.solve_tolf
360
                        [dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
361
                    else
362
363
364
                        flag1 = 0;
                        dx = phat - ya;
                    end;
365
                    if (flag1>0 || reduced)
366
367
368
369
370
371
372
                        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;
373
                        ilu_setup.droptol = ilu_setup.droptol/10;
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
                        reduced = 0;
                    else
                        ya = ya + lambda*dx;
                        if(is_dynamic)
                            y(it_,y_index_eq) = ya';
                        else
                            y(y_index_eq) = ya';
                        end;
                    end;
                end;
            else
                disp('unknown option : ');
                if(is_dynamic)
                    disp(['options_.stack_solve_algo = ' num2str(stack_solve_algo) ' not implemented']);
                else
Ferhat Mihoubi's avatar
Ferhat Mihoubi committed
389
                    disp(['options_.solve_algo = ' num2str(options.solve_algo) ' not implemented']);
390
391
392
393
394
395
396
397
398
399
                end;
                info = -Block_Num*10;
                return;
            end;
            iter=iter+1;
            max_resa = max_res;
        end
    end
    if cvg==0
        if(cutoff == 0)
400
            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);
401
        else
402
            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);
403
404
405
406
407
408
409
410
411
412
413
414
415
416
        end;
        if(is_dynamic)
            oo_.deterministic_simulation.status = 0;
            oo_.deterministic_simulation.error = max_res;
            oo_.deterministic_simulation.iterations = iter;
            oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed.
            oo_.deterministic_simulation.block(Block_Num).error = max_res;
            oo_.deterministic_simulation.block(Block_Num).iterations = iter;
        end;
        info = -Block_Num*10;
        return;
    end
end
if(is_dynamic)
417
    info = 1;
418
419
420
421
422
423
    oo_.deterministic_simulation.status = 1;
    oo_.deterministic_simulation.error = max_res;
    oo_.deterministic_simulation.iterations = iter;
    oo_.deterministic_simulation.block(Block_Num).status = 1;
    oo_.deterministic_simulation.block(Block_Num).error = max_res;
    oo_.deterministic_simulation.block(Block_Num).iterations = iter;
424
else
425
    info = 0;
426
427
428
end;
return;

429
function [err, G]=local_fname(yl, x, params, steady_state, y, y_index_eq, fname, is_csolve)
430
y(y_index_eq) = yl;
431
[err, y, G] = feval(fname, y, x, params, steady_state, 0);
432
433
434
if(is_csolve)
    G = full(G);
end;