4.0: merged changesets 2078, 2089 and 2090 from trunk

* new debug messages in dynare_solve.m
* solve_algo=2 exits immediately after failure to solve one block
* new solve_algo=4 with different way of dealing with badly scaled or nearly singular Jacobian
* solve_algo=0 fails if Matlab's Optimization Toolbox not present
* dynare_solve fails if solve_algo not in [0,1,2,3,4]


git-svn-id: https://www.dynare.org/svn/dynare/branches/4.0@2095 ac1d8469-bf42-47a9-8791-bf33cf982152

matlab/dynare_solve.m

@@ -38,7 +38,11 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
   options_ = set_default_option(options_,'solve_algo',2);
   info = 0;
   if options_.solve_algo == 0
-    if ~isempty(which('fsolve'))
+    if exist('OCTAVE_VERSION') || isempty(ver('optim'))
+      % Note that fsolve() exists under Octave, but has a different syntax
+      % So we fail for the moment under Octave, until we add the corresponding code
+      error('DYNARE_SOLVE: you can''t use solve_algo=0 since you don''t have Matlab''s Optimization Toolbox')
+    end
     options=optimset('fsolve');
     options.MaxFunEvals = 50000;
     options.MaxIter = 2000;
@@ -55,16 +59,10 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
     else
       info = 1;
     end
-      return
-    else 
-      options_.solve_algo = 1;
-    end
-  end
-
-  if options_.solve_algo == 1
+  elseif options_.solve_algo == 1
     nn = size(x,1);
     [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
-  elseif options_.solve_algo == 2
+  elseif options_.solve_algo == 2 || options_.solve_algo == 4
     nn = size(x,1) ;
     tolf = options_.solve_tolf ;
 
@@ -84,8 +82,6 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
       error('exiting ...')
     end
     
-%    f = 0.5*fvec'*fvec ;
-
     if max(abs(fvec)) < tolf
       return ;
     end
@@ -102,15 +98,30 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
 
     [j1,j2,r,s] = dmperm(fjac);
     
+    if options_.debug
+      disp(['DYNARE_SOLVE (solve_algo=2|4): number of blocks = ' num2str(length(r))]);
+    end
+    
     for i=length(r)-1:-1:1
+      if options_.debug
+        disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]);
+      end
+      if options_.solve_algo == 2
         [x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag,varargin{:});
-      if info & options_.debug
-	error(sprintf('Solve block = %d check = %d\n',i,info));
+      else % solve_algo=4
+        [x,info]=solve2(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag,varargin{:});
+      end
+      if info
+        return
       end
     end
     fvec = feval(func,x,varargin{:});
     if max(abs(fvec)) > tolf
+      if options_.solve_algo == 2
         [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
+      else % solve_algo=4
+        [x,info]=solve2(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
+      end
     end
   elseif options_.solve_algo == 3
     if jacobian_flag
@@ -118,9 +129,6 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
     else
       [x,info] = csolve(func,x,[],1e-6,500,varargin{:});
     end
+  else
+    error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4]')
   end
-%    fvec1 = feval(func,x,varargin{:})
-
-  % 08/28/03 MJ add a final call to solve1 for solve_algo == 1 in case
-  %             initvals generates 'false' zeros in the Jacobian
-  
\ No newline at end of file

matlab/solve2.m

+function [x,check] = solve2(func,x,j1,j2,jacobian_flag,varargin)
+% function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin)
+% Solves systems of non linear equations of several variables
+%
+% This solver is used for solve_algo = ...
+%
+% This is a modified version of solve1:
+% In solve1, before proceeding to the Newton step, we first check
+% that the condition number is reasonable, and we use an alternate
+% step formula if it is not the case.
+% Here, we first do the Newton step, then we check if the left
+% matrix division returned a warning (in case of badly scale or
+% nearly singular jacobian) in which case we use the alternate
+% step formula.
+%
+% INPUTS
+%    func:            name of the function to be solved
+%    x:               guess values
+%    j1:              equations index for which the model is solved
+%    j2:              unknown variables index
+%    jacobian_flag=1: jacobian given by the 'func' function
+%    jacobian_flag=0: jacobian obtained numerically
+%    varargin:        list of arguments following jacobian_flag
+%    
+% OUTPUTS
+%    x:               results
+%    check=1:         the model can not be solved
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% Copyright (C) 2001-2008 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/>.
+
+  global M_ options_ fjac  
+
+  nn = length(j1);
+  fjac = zeros(nn,nn) ;
+  g = zeros(nn,1) ;
+
+  tolf = options_.solve_tolf ;
+  tolx = options_.solve_tolx;
+  tolmin = tolx ;
+
+  stpmx = 100 ;
+  maxit = options_.solve_maxit ;
+
+  check = 0 ;
+
+  fvec = feval(func,x,varargin{:});
+  fvec = fvec(j1);
+  
+  i = find(~isfinite(fvec));
+  
+  if ~isempty(i)
+    disp(['STEADY:  numerical initial values incompatible with the following' ...
+	  ' equations'])
+    disp(j1(i)')
+  end
+  
+  f = 0.5*fvec'*fvec ;
+
+  if max(abs(fvec)) < tolf
+    return ;
+  end
+
+  stpmax = stpmx*max([sqrt(x'*x);nn]) ;
+  first_time = 1;
+  for its = 1:maxit
+    if jacobian_flag
+      [fvec,fjac] = feval(func,x,varargin{:});
+      fvec = fvec(j1);
+      fjac = fjac(j1,j2);
+    else
+      dh = max(abs(x(j2)),options_.gstep*ones(nn,1))*eps^(1/3);
+      
+      for j = 1:nn
+	xdh = x ;
+	xdh(j2(j)) = xdh(j2(j))+dh(j) ;
+	t = feval(func,xdh,varargin{:});
+	fjac(:,j) = (t(j1) - fvec)./dh(j) ;
+	g(j) = fvec'*fjac(:,j) ;
+      end
+    end
+
+    g = (fvec'*fjac)';
+    if options_.debug
+      disp(['cond(fjac) ' num2str(cond(fjac))])
+    end
+    M_.unit_root = 0;
+    if M_.unit_root
+      if first_time
+	    first_time = 0;
+	    [q,r,e]=qr(fjac);
+	    n = sum(abs(diag(r)) < 1e-12);
+	    fvec = q'*fvec;
+	    p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
+	    disp(' ')
+	    disp('STEADY with unit roots:')
+	    disp(' ')
+	    if n > 0
+	      disp(['   The following variable(s) kept their value given in INITVAL' ...
+		    ' or ENDVAL'])
+	      disp(char(e(:,end-n+1:end)'*M_.endo_names))
+        else
+	      disp('   STEADY can''t find any unit root!')
+	    end
+      else
+	    [q,r]=qr(fjac*e);
+	    fvec = q'*fvec;
+	    p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
+      end	
+    else
+      lastwarn('');
+      p = -fjac\fvec;
+      if ~isempty(lastwarn)
+        fjac2=fjac'*fjac;
+        p=-(fjac2+sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec);
+      end
+    end
+    xold = x ;
+    fold = f ;
+
+    [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin{:});
+
+    if options_.debug
+      disp([its f])
+      disp([xold x])
+    end
+      
+    if check > 0
+      den = max([f;0.5*nn]) ;
+      if max(abs(g).*max([abs(x(j2)') ones(1,nn)])')/den < tolmin
+	    return
+      else
+	    disp (' ')
+	    disp (['SOLVE: Iteration ' num2str(its)])
+	    disp (['Spurious convergence.'])
+	    disp (x)
+	    return
+      end
+
+      if max(abs(x(j2)-xold(j2))./max([abs(x(j2)') ones(1,nn)])') < tolx
+	disp (' ')
+	disp (['SOLVE: Iteration ' num2str(its)])
+	disp (['Convergence on dX.'])
+	disp (x)
+	return
+      end
+    elseif max(abs(fvec)) < tolf
+      return
+    end
+  end
+  
+  check = 1;
+  disp(' ')
+  disp('SOLVE: maxit has been reached')