solve_one_boundary

 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
   steady_state        [vector]        steady state of the model
   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
   indirect_call      [integer]        (0) direct call to the fname
                                       (1) indirect call via the
                                       local_fname wrapper
 OUTPUTS                                    
   y                  [matrix]         All endogenous variables of the model      
   info               [integer]        >=0 no error
                                       <0 error
  
 ALGORITHM
   Newton with LU or GMRES or BicGstab for dynamic block
    
 SPECIAL REQUIREMENTS
   none.

0001 function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, ...
0002                                         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)
0003 % Computes the deterministic simulation of a block of equation containing
0004 % lead or lag variables
0005 %
0006 % INPUTS
0007 %   fname               [string]        name of the file containing the block
0008 %                                       to simulate
0009 %   y                   [matrix]        All the endogenous variables of the model
0010 %   x                   [matrix]        All the exogenous variables of the model
0011 %   params              [vector]        All the parameters of the model
0012 %   steady_state        [vector]        steady state of the model
0013 %   y_index_eq          [vector of int] The index of the endogenous variables of
0014 %                                       the block
0015 %   nze                 [integer]       number of non-zero elements in the
0016 %                                       jacobian matrix
0017 %   periods             [integer]       number of simulation periods
0018 %   is_linear           [integer]       if is_linear=1 the block is linear
0019 %                                       if is_linear=0 the block is not linear
0020 %   Block_Num           [integer]       block number
0021 %   y_kmin              [integer]       maximum number of lag in the model
0022 %   maxit_              [integer]       maximum number of iteration in Newton
0023 %   solve_tolf          [double]        convergence criteria
0024 %   lambda              [double]        initial value of step size in
0025 %   Newton
0026 %   cutoff              [double]        cutoff to correct the direction in Newton in case
0027 %                                       of singular jacobian matrix
0028 %   stack_solve_algo    [integer]       linear solver method used in the
0029 %                                       Newton algorithm :
0030 %                                            - 1 sparse LU
0031 %                                            - 2 GMRES
0032 %                                            - 3 BicGStab
0033 %                                            - 4 Optimal path length
0034 %   forward_backward    [integer]       The block has to be solve forward
0035 %                                       (1) or backward (0)
0036 %   is_dynamic          [integer]       (1) The block belong to the dynamic
0037 %                                           file and the oo_.deterministic_simulation field has to be uptated
0038 %                                       (0) The block belong to the static
0039 %                                           file and th oo_.deteerminstic
0040 %                                           field remains unchanged
0041 %   verbose            [integer]        (0) iterations are not printed
0042 %                                       (1) iterations are printed
0043 %   indirect_call      [integer]        (0) direct call to the fname
0044 %                                       (1) indirect call via the
0045 %                                       local_fname wrapper
0046 % OUTPUTS
0047 %   y                  [matrix]         All endogenous variables of the model
0048 %   info               [integer]        >=0 no error
0049 %                                       <0 error
0050 %
0051 % ALGORITHM
0052 %   Newton with LU or GMRES or BicGstab for dynamic block
0053 %
0054 % SPECIAL REQUIREMENTS
0055 %   none.
0056 %
0057 
0058 % Copyright (C) 1996-2012 Dynare Team
0059 %
0060 % This file is part of Dynare.
0061 %
0062 % Dynare is free software: you can redistribute it and/or modify
0063 % it under the terms of the GNU General Public License as published by
0064 % the Free Software Foundation, either version 3 of the License, or
0065 % (at your option) any later version.
0066 %
0067 % Dynare is distributed in the hope that it will be useful,
0068 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0069 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0070 % GNU General Public License for more details.
0071 %
0072 % You should have received a copy of the GNU General Public License
0073 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0074 
0075 
0076 Blck_size=size(y_index_eq,2);
0077 g2 = [];
0078 g3 = [];
0079 correcting_factor=0.01;
0080 luinc_tol=1e-10;
0081 max_resa=1e100;
0082 reduced = 0;
0083 if(forward_backward)
0084     incr = 1;
0085     start = y_kmin+1;
0086     finish = periods+y_kmin;
0087 else
0088     incr = -1;
0089     start = periods+y_kmin;
0090     finish = y_kmin+1;
0091 end
0092 %lambda=1;
0093 for it_=start:incr:finish
0094     cvg=0;
0095     iter=0;
0096     g1=spalloc( Blck_size, Blck_size, nze);
0097     while ~(cvg==1 || iter>maxit_),
0098         if(is_dynamic)
0099             [r, y, g1, g2, g3] = feval(fname, y, x, params, steady_state, ...
0100                                        it_, 0);
0101         else
0102             [r, y, g1] = feval(fname, y, x, params);
0103         end;
0104         if(~isreal(r))
0105             max_res=(-(max(max(abs(r))))^2)^0.5;
0106         else
0107             max_res=max(max(abs(r)));
0108         end;
0109         %['max_res=' num2str(max_res) ' Block_Num=' int2str(Block_Num) ' it_=' int2str(it_)]
0110         %disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)]);
0111         
0112         %     fjac = zeros(Blck_size, Blck_size);
0113         %     disp(['Blck_size=' int2str(Blck_size) ' it_=' int2str(it_)]);
0114         %     dh = max(abs(y(it_, y_index_eq)),options_.gstep*ones(1, Blck_size))*eps^(1/3);
0115         %     fvec = r;
0116         %       for j = 1:Blck_size
0117         %         ydh = y ;
0118         %           ydh(it_, y_index_eq(j)) = ydh(it_, y_index_eq(j)) + dh(j)  ;
0119         %           [t, y1, g11, g21, g31]=feval(fname, ydh, x, params, it_, 0);
0120         %           fjac(:,j) = (t - fvec)./dh(j) ;
0121         %       end;
0122         %     diff = g1 -fjac;
0123         %     diff
0124         %     disp('g1');
0125         %     disp([num2str(g1,'%4.5f')]);
0126         %     disp('fjac');
0127         %     disp([num2str(fjac,'%4.5f')]);
0128         %     [c_max, i_c_max] = max(abs(diff));
0129         %     [l_c_max, i_r_max] = max(c_max);
0130         %     disp(['maximum element row=' int2str(i_c_max(i_r_max)) ' and column=' int2str(i_r_max) ' value = ' num2str(l_c_max)]);
0131         %     equation = i_c_max(i_r_max);
0132         %     variable = i_r_max;
0133         %     variable
0134         %     mod(variable, Blck_size)
0135         %     disp(['equation ' int2str(equation) ' and variable ' int2str(y_index_eq(variable)) ' ' M_.endo_names(y_index_eq(variable), :)]);
0136         %     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))]);
0137         %     %return;
0138         %     %g1 = fjac;
0139 
0140         
0141         
0142         
0143         
0144         if(verbose==1)
0145             disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)]);
0146             if(is_dynamic)
0147                 disp([M.endo_names(y_index_eq,:) num2str([y(it_,y_index_eq)' r g1])]);
0148             else
0149                 disp([M.endo_names(y_index_eq,:) num2str([y(y_index_eq) r g1])]);
0150             end;
0151         end;
0152         if(~isreal(max_res) || isnan(max_res))
0153             cvg = 0;
0154         elseif(is_linear && iter>0)
0155             cvg = 1;
0156         else
0157             cvg=(max_res<solve_tolf);
0158         end;
0159         if(~cvg)
0160             if(iter>0)
0161                 if(~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1))
0162                     if(isnan(max_res) || (max_resa<max_res && iter>0))
0163                         detJ=det(g1a);
0164                         if(abs(detJ)<1e-7)
0165                             max_factor=max(max(abs(g1a)));
0166                             ze_elem=sum(diag(g1a)<cutoff);
0167                             disp([num2str(full(ze_elem),'%d') ' elements on the Jacobian diagonal are below the cutoff (' num2str(cutoff,'%f') ')']);
0168                             if(correcting_factor<max_factor)
0169                                 correcting_factor=correcting_factor*4;
0170                                 disp(['The Jacobain matrix is singular, det(Jacobian)=' num2str(detJ,'%f') '.']);
0171                                 disp(['    trying to correct the Jacobian matrix:']);
0172                                 disp(['    correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')]);
0173                                 dx = - r/(g1+correcting_factor*speye(Blck_size));
0174                                 %dx = -b'/(g1+correcting_factor*speye(Blck_size))-ya_save;
0175                                 y(it_,y_index_eq)=ya_save+lambda*dx;
0176                                 continue;
0177                             else
0178                                 disp('The singularity of the jacobian matrix could not be corrected');
0179                                 info = -Block_Num*10;
0180                                 return;
0181                             end;
0182                         end;
0183                     elseif(lambda>1e-8)
0184                         lambda=lambda/2;
0185                         reduced = 1;
0186                         disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);
0187                         if(is_dynamic)
0188                             y(it_,y_index_eq)=ya_save-lambda*dx;
0189                         else
0190                             y(y_index_eq)=ya_save-lambda*dx;
0191                         end;
0192                         continue;
0193                     else
0194                         if(cutoff == 0)
0195                             fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.maxit_".\n',Block_Num, it_, iter);
0196                         else
0197                             fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.maxit_" or set "cutoff=0" in model options.\n',Block_Num, it_, iter);
0198                         end;
0199                         if(is_dynamic)
0200                             oo_.deterministic_simulation.status = 0;
0201                             oo_.deterministic_simulation.error = max_res;
0202                             oo_.deterministic_simulation.iterations = iter;
0203                             oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed.
0204                             oo_.deterministic_simulation.block(Block_Num).error = max_res;
0205                             oo_.deterministic_simulation.block(Block_Num).iterations = iter;
0206                         end;
0207                         info = -Block_Num*10;
0208                         return;
0209                     end;
0210                 else
0211                     if(lambda<1)
0212                         lambda=max(lambda*2, 1);
0213                     end;
0214                 end;
0215             end;
0216             if(is_dynamic)
0217                 ya = y(it_,y_index_eq)';
0218             else
0219                 ya = y(y_index_eq);
0220             end;
0221             ya_save=ya;
0222             g1a=g1;
0223             if(~is_dynamic && options.solve_algo == 0)
0224                 if (verbose == 1)
0225                     disp('steady: fsolve');
0226                 end
0227                 if ~exist('OCTAVE_VERSION')
0228                     if ~user_has_matlab_license('optimization_toolbox')
0229                         error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
0230                     end
0231                 end
0232                 options=optimset('fsolve');
0233                 options.MaxFunEvals = 50000;
0234                 options.MaxIter = 2000;
0235                 options.TolFun=1e-8;
0236                 options.Display = 'iter';
0237                 options.Jacobian = 'on';
0238                 if ~exist('OCTAVE_VERSION')
0239                     [yn,fval,exitval,output] = fsolve(@local_fname, y(y_index_eq), ...
0240                                                       options, x, params, steady_state, y, y_index_eq, fname, 0);
0241                 else
0242                     % Under Octave, use a wrapper, since fsolve() does not have a 4th arg
0243                     func = @(z) local_fname(z, x, params, steady_state, y, y_index_eq, fname, 0);
0244                     % The Octave version of fsolve does not converge when it starts from the solution
0245                     fvec = feval(func,y(y_index_eq));
0246                     if max(abs(fvec)) >= options.solve_tolf
0247                         [yn,fval,exitval,output] = fsolve(func,y(y_index_eq),options);
0248                     else
0249                         yn = y(y_index_eq);
0250                         exitval = 3;
0251                     end;
0252                 end
0253                     
0254                 y(y_index_eq) = yn;
0255                 if exitval > 0
0256                     info = 0;
0257                 else
0258                     info = -Block_Num*10;
0259                 end
0260             elseif((~is_dynamic && options.solve_algo==2) || (is_dynamic && stack_solve_algo==4))
0261                 if (verbose == 1 && ~is_dynamic)
0262                     disp('steady: LU + lnsrch1');
0263                 end
0264                 lambda=1;
0265                 stpmx = 100 ;
0266                 if (is_dynamic)
0267                     stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
0268                 else
0269                     stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
0270                 end;
0271                 nn=1:size(y_index_eq,2);
0272                 g = (r'*g1)';
0273                 f = 0.5*r'*r;
0274                 p = -g1\r ;
0275                 if (is_dynamic)
0276                     [ya,f,r,check]=lnsrch1(y(it_,:)',f,g,p,stpmax, ...
0277                                            'lnsrch1_wrapper_one_boundary',nn, ...
0278                                            y_index_eq, y_index_eq, fname, y, x, params, steady_state, it_);
0279                     dx = ya' - y(it_, :);
0280                 else
0281                     [ya,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, ...
0282                                            params, steady_state,0);
0283                     dx = ya - y(y_index_eq);
0284                 end;
0285                 
0286                 if(is_dynamic)
0287                     y(it_,:) = ya';
0288                 else
0289                     y = ya';
0290                 end;
0291             elseif(~is_dynamic && options.solve_algo==3)
0292                 if (verbose == 1)
0293                     disp('steady: csolve');
0294                 end
0295                 [yn,info] = csolve(@local_fname, y(y_index_eq),@ ...
0296                                    local_fname,1e-6,500, x, params, steady_state, y, y_index_eq, fname, 1);
0297                 dx = ya - yn;
0298                 y(y_index_eq) = yn;
0299             elseif((stack_solve_algo==1 && is_dynamic) || (stack_solve_algo==0 && is_dynamic) || (~is_dynamic && (options.solve_algo==1 || options.solve_algo==6))),
0300                 if (verbose == 1 && ~is_dynamic)
0301                     disp('steady: Sparse LU ');
0302                 end
0303                 dx =  g1\r;
0304                 ya = ya - lambda*dx;
0305                 if(is_dynamic)
0306                     y(it_,y_index_eq) = ya';
0307                 else
0308                     y(y_index_eq) = ya;
0309                 end;
0310             elseif((stack_solve_algo==2 && is_dynamic) || (options.solve_algo==7 && ~is_dynamic)),
0311                 flag1=1;
0312                 if exist('OCTAVE_VERSION')
0313                     error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=7 since GMRES is not implemented in Octave')
0314                 end
0315                 if (verbose == 1 && ~is_dynamic)
0316                     disp('steady: GMRES ');
0317                 end
0318                 while(flag1>0)
0319                     [L1, U1]=luinc(g1,luinc_tol);
0320                     [dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
0321                     if (flag1>0 || reduced)
0322                         if(flag1==1)
0323                             disp(['Error in simul: No convergence inside GMRES after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]);
0324                         elseif(flag1==2)
0325                             disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')]);
0326                         elseif(flag1==3)
0327                             disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')]);
0328                         end;
0329                         luinc_tol = luinc_tol/10;
0330                         reduced = 0;
0331                     else
0332                         ya = ya + lambda*dx;
0333                         if(is_dynamic)
0334                             y(it_,y_index_eq) = ya';
0335                         else
0336                             y(y_index_eq) = ya';
0337                         end;
0338                     end;
0339                 end;
0340             elseif((stack_solve_algo==3 && is_dynamic) || (options.solve_algo==8 && ~is_dynamic)),
0341                 flag1=1;
0342                 if (verbose == 1 && ~is_dynamic)
0343                     disp('steady: BiCGStab');
0344                 end
0345                 while(flag1>0)
0346                     [L1, U1]=luinc(g1,luinc_tol);
0347                     phat = ya - U1 \ (L1 \ r);
0348                     if(is_dynamic)
0349                         y(it_,y_index_eq) = phat;
0350                     else
0351                         y(y_index_eq) = phat;
0352                     end;
0353                     if(is_dynamic)
0354                         [r, y, g1, g2, g3] = feval(fname, y, x, params, ...
0355                                                    steady_state, it_, 0);
0356                     else
0357                         [r, y, g1] = feval(fname, y, x, params);
0358                     end;
0359                     if max(abs(r)) >= options.solve_tolf
0360                         [dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
0361                     else
0362                         flag1 = 0;
0363                         dx = phat - ya;
0364                     end;
0365                     if (flag1>0 || reduced)
0366                         if(flag1==1)
0367                             disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]);
0368                         elseif(flag1==2)
0369                             disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')]);
0370                         elseif(flag1==3)
0371                             disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')]);
0372                         end;
0373                         luinc_tol = luinc_tol/10;
0374                         reduced = 0;
0375                     else
0376                         ya = ya + lambda*dx;
0377                         if(is_dynamic)
0378                             y(it_,y_index_eq) = ya';
0379                         else
0380                             y(y_index_eq) = ya';
0381                         end;
0382                     end;
0383                 end;
0384             else
0385                 disp('unknown option : ');
0386                 if(is_dynamic)
0387                     disp(['options_.stack_solve_algo = ' num2str(stack_solve_algo) ' not implemented']);
0388                 else
0389                     disp(['options_.solve_algo = ' num2str(options.solve_algo) ' not implemented']);
0390                 end;
0391                 info = -Block_Num*10;
0392                 return;
0393             end;
0394             iter=iter+1;
0395             max_resa = max_res;
0396         end
0397     end
0398     if cvg==0
0399         if(cutoff == 0)
0400             fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.maxit_\".\n',Block_Num, it_,iter);
0401         else
0402             fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.maxit_" or set "cutoff=0" in model options.\n',Block_Num, it_,iter);
0403         end;
0404         if(is_dynamic)
0405             oo_.deterministic_simulation.status = 0;
0406             oo_.deterministic_simulation.error = max_res;
0407             oo_.deterministic_simulation.iterations = iter;
0408             oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed.
0409             oo_.deterministic_simulation.block(Block_Num).error = max_res;
0410             oo_.deterministic_simulation.block(Block_Num).iterations = iter;
0411         end;
0412         info = -Block_Num*10;
0413         return;
0414     end
0415 end
0416 if(is_dynamic)
0417     info = 1;
0418     oo_.deterministic_simulation.status = 1;
0419     oo_.deterministic_simulation.error = max_res;
0420     oo_.deterministic_simulation.iterations = iter;
0421     oo_.deterministic_simulation.block(Block_Num).status = 1;
0422     oo_.deterministic_simulation.block(Block_Num).error = max_res;
0423     oo_.deterministic_simulation.block(Block_Num).iterations = iter;
0424 else
0425     info = 0;
0426 end;
0427 return;
0428 
0429 function [err, G]=local_fname(yl, x, params, steady_state, y, y_index_eq, fname, is_csolve)
0430 y(y_index_eq) = yl;
0431 [err, y, G] = feval(fname, y, x, params, steady_state, 0);
0432 if(is_csolve)
0433     G = full(G);
0434 end;