solve_one_boundary.m

function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, ...
                                        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)
% 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.
%  

% Copyright (C) 1996-2013 Dynare Team
%
% 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;
ilu_setup.droptol=1e-10;
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
%lambda=1;
for it_=start:incr:finish
    cvg=0;
    iter=0;
    g1=spalloc( Blck_size, Blck_size, nze);
    while ~(cvg==1 || iter>maxit_),
        if(is_dynamic)
            [r, y, g1, g2, g3] = feval(fname, y, x, params, steady_state, ...
                                       it_, 0);
        else
            [r, y, g1] = feval(fname, y, x, params);
        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
        %         ydh = y ;
        %           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)
                disp([M.endo_names(y_index_eq,:) num2str([y(it_,y_index_eq)' r g1])]);
            else
                disp([M.endo_names(y_index_eq,:) num2str([y(y_index_eq) r g1])]);
            end;
        end;
        if(~isreal(max_res) || isnan(max_res))
            cvg = 0;
        elseif(is_linear && iter>0)
            cvg = 1;
        else
            cvg=(max_res<solve_tolf);
        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))
                        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)
                            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);
                        else
                            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);
                        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;
            if(~is_dynamic && options.solve_algo == 0)
                if (verbose == 1)
                    disp('steady: fsolve');
                end
                if ~exist('OCTAVE_VERSION')
                    if ~user_has_matlab_license('optimization_toolbox')
                        error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
                    end
                end
                options=optimset('fsolve');
                options.MaxFunEvals = 50000;
                options.MaxIter = 2000;
                options.TolFun=1e-8;
                options.Display = 'iter';
                options.Jacobian = 'on';
                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));
                    if max(abs(fvec)) >= options.solve_tolf
                        [yn,fval,exitval,output] = fsolve(func,y(y_index_eq),options);
                    else
                        yn = y(y_index_eq);
                        exitval = 3;
                    end;
                end
                    
                y(y_index_eq) = yn;
                if exitval > 0
                    info = 0;
                else
                    info = -Block_Num*10;
                end
            elseif((~is_dynamic && options.solve_algo==2) || (is_dynamic && stack_solve_algo==4))
                if (verbose == 1 && ~is_dynamic)
                    disp('steady: LU + lnsrch1');
                end
                lambda=1;
                stpmx = 100 ;
                if (is_dynamic)
                    stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
                else
                    stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
                end;
                nn=1:size(y_index_eq,2);
                g = (r'*g1)';
                f = 0.5*r'*r;
                p = -g1\r ;
                if (is_dynamic)
                    [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_);
                    dx = ya' - y(it_, :);
                else
                    [ya,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, ...
                                           params, steady_state,0);
                    dx = ya - y(y_index_eq);
                end;
                
                if(is_dynamic)
                    y(it_,:) = ya';
                else
                    y = ya';
                end;
            elseif(~is_dynamic && options.solve_algo==3)
                if (verbose == 1)
                    disp('steady: csolve');
                end
                [yn,info] = csolve(@local_fname, y(y_index_eq),@ ...
                                   local_fname,1e-6,500, x, params, steady_state, y, y_index_eq, fname, 1);
                dx = ya - yn;
                y(y_index_eq) = yn;
            elseif((stack_solve_algo==1 && is_dynamic) || (stack_solve_algo==0 && is_dynamic) || (~is_dynamic && (options.solve_algo==1 || options.solve_algo==6))),
                if (verbose == 1 && ~is_dynamic)
                    disp('steady: Sparse LU ');
                end
                dx =  g1\r;
                ya = ya - lambda*dx;
                if(is_dynamic)
                    y(it_,y_index_eq) = ya';
                else
                    y(y_index_eq) = ya;
                end;
            elseif((stack_solve_algo==2 && is_dynamic) || (options.solve_algo==7 && ~is_dynamic)),
                flag1=1;
                if exist('OCTAVE_VERSION')
                    error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=7 since GMRES is not implemented in Octave')
                end
                if (verbose == 1 && ~is_dynamic)
                    disp('steady: GMRES ');
                end
                while(flag1>0)
                    [L1, U1]=ilu(g1,ilu_setup);
                    [dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
                    if (flag1>0 || reduced)
                        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;
                        ilu_setup.droptol = ilu_setup.droptol/10;
                        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;
            elseif((stack_solve_algo==3 && is_dynamic) || (options.solve_algo==8 && ~is_dynamic)),
                flag1=1;
                if (verbose == 1 && ~is_dynamic)
                    disp('steady: BiCGStab');
                end
                while(flag1>0)
                    [L1, U1]=ilu(g1,ilu_setup);
                    phat = ya - U1 \ (L1 \ r);
                    if(is_dynamic)
                        y(it_,y_index_eq) = phat;
                    else
                        y(y_index_eq) = phat;
                    end;
                    if(is_dynamic)
                        [r, y, g1, g2, g3] = feval(fname, y, x, params, ...
                                                   steady_state, it_, 0);
                    else
                        [r, y, g1] = feval(fname, y, x, params);
                    end;
                    if max(abs(r)) >= options.solve_tolf
                        [dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
                    else
                        flag1 = 0;
                        dx = phat - ya;
                    end;
                    if (flag1>0 || reduced)
                        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;
                        ilu_setup.droptol = ilu_setup.droptol/10;
                        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
                    disp(['options_.solve_algo = ' num2str(options.solve_algo) ' not implemented']);
                end;
                info = -Block_Num*10;
                return;
            end;
            iter=iter+1;
            max_resa = max_res;
        end
    end
    if cvg==0
        if(cutoff == 0)
            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);
        else
            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);
        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)
    info = 1;
    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;
else
    info = 0;
end;
return;

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