lnsrch1

 function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin)
 Computes the optimal step by minimizing the residual sum of squares

 INPUTS
   xold:     actual point
   fold:     residual sum of squares at the point xold
   g:        gradient
   p:        Newton direction
   stpmax:   maximum step
   func:     name of the function
   j1:       equations index to be solved
   j2:       unknowns index
   varargin: list of arguments following j2

 OUTPUTS
   x:        chosen point
   f:        residual sum of squares value for a given x
   fvec:     residuals vector
   check=1:  problem of the looping which continues indefinitely

 
 SPECIAL REQUIREMENTS
   none

0001 function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin)
0002 % function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin)
0003 % Computes the optimal step by minimizing the residual sum of squares
0004 %
0005 % INPUTS
0006 %   xold:     actual point
0007 %   fold:     residual sum of squares at the point xold
0008 %   g:        gradient
0009 %   p:        Newton direction
0010 %   stpmax:   maximum step
0011 %   func:     name of the function
0012 %   j1:       equations index to be solved
0013 %   j2:       unknowns index
0014 %   varargin: list of arguments following j2
0015 %
0016 % OUTPUTS
0017 %   x:        chosen point
0018 %   f:        residual sum of squares value for a given x
0019 %   fvec:     residuals vector
0020 %   check=1:  problem of the looping which continues indefinitely
0021 %
0022 %
0023 % SPECIAL REQUIREMENTS
0024 %   none
0025 
0026 % Copyright (C) 2001-2010 Dynare Team
0027 %
0028 % This file is part of Dynare.
0029 %
0030 % Dynare is free software: you can redistribute it and/or modify
0031 % it under the terms of the GNU General Public License as published by
0032 % the Free Software Foundation, either version 3 of the License, or
0033 % (at your option) any later version.
0034 %
0035 % Dynare is distributed in the hope that it will be useful,
0036 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0037 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0038 % GNU General Public License for more details.
0039 %
0040 % You should have received a copy of the GNU General Public License
0041 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0042 
0043 global options_
0044 
0045 alf = 1e-4 ;
0046 tolx = options_.solve_tolx;
0047 alam = 1;
0048 
0049 x = xold;
0050 nn = length(j2);
0051 summ = sqrt(sum(p.*p)) ;
0052 if ~isfinite(summ)
0053     error(['Some element of Newton direction isn''t finite. Jacobian maybe' ...
0054            ' singular or there is a problem with initial values'])
0055 end
0056 
0057 if summ > stpmax
0058     p=p.*stpmax/summ ;
0059 end
0060 
0061 slope = g'*p ;
0062 
0063 test = max(abs(p)'./max([abs(xold(j2))';ones(1,nn)])) ;
0064 alamin = tolx/test ;
0065 
0066 if alamin > 0.1
0067     alamin = 0.1;
0068 end
0069 
0070 while 1
0071     if alam < alamin
0072         check = 1 ;
0073         return
0074     end
0075     
0076     x(j2) = xold(j2) + (alam*p) ;
0077     fvec = feval(func,x,varargin{:}) ;
0078     fvec = fvec(j1);
0079     f = 0.5*fvec'*fvec ;
0080 
0081     if any(isnan(fvec))
0082         alam = alam/2 ;
0083         alam2 = alam ;
0084         f2 = f ;
0085         fold2 = fold ;
0086     else
0087 
0088         if f <= fold+alf*alam*slope
0089             check = 0;
0090             break ;
0091         else
0092             if alam == 1
0093                 tmplam = -slope/(2*(f-fold-slope)) ;
0094             else
0095                 rhs1 = f-fold-alam*slope ;
0096                 rhs2 = f2-fold2-alam2*slope ;
0097                 a = (rhs1/(alam^2)-rhs2/(alam2^2))/(alam-alam2) ;
0098                 b = (-alam2*rhs1/(alam^2)+alam*rhs2/(alam2^2))/(alam-alam2) ;
0099                 if a == 0
0100                     tmplam = -slope/(2*b) ;
0101                 else
0102                     disc = (b^2)-3*a*slope ;
0103 
0104                     if disc < 0
0105                         error ('Roundoff problem in nlsearch') ;
0106                     else
0107                         tmplam = (-b+sqrt(disc))/(3*a) ;
0108                     end
0109 
0110                 end
0111 
0112                 if tmplam > 0.5*alam
0113                     tmplam = 0.5*alam;
0114                 end
0115 
0116             end
0117 
0118             alam2 = alam ;
0119             f2 = f ;
0120             fold2 = fold ;
0121             alam = max([tmplam;(0.1*alam)]) ;
0122         end
0123     end
0124 end
0125 
0126 % 01/14/01 MJ lnsearch is now a separate function
0127 % 01/12/03 MJ check for finite summ to avoid infinite loop when Jacobian
0128 %             is singular or model is denormalized