function [x,rc] = csolve(FUN,x,gradfun,crit,itmax,varargin)
FUN should be written so that any parametric arguments are packed in to xp,
and so that if presented with a matrix x, it produces a return value of
same dimension of x. The number of rows in x and FUN(x) are always the
same. The number of columns is the number of different input arguments
at which FUN is to be evaluated.
gradfun: string naming the function called to evaluate the gradient matrix. If this
is null (i.e. just "[]"), a numerical gradient is used instead.
crit: if the sum of absolute values that FUN returns is less than this,
the equation is solved.
itmax: the solver stops when this number of iterations is reached, with rc=4
varargin: in this position the user can place any number of additional arguments, all
of which are passed on to FUN and gradfun (when it is non-empty) as a list of
arguments following x.
rc: 0 means normal solution, 1 and 3 mean no solution despite extremely fine adjustments
in step length (very likely a numerical problem, or a discontinuity). 4 means itmax
termination.0001 function [x,rc] = csolve(FUN,x,gradfun,crit,itmax,varargin) 0002 %function [x,rc] = csolve(FUN,x,gradfun,crit,itmax,varargin) 0003 % FUN should be written so that any parametric arguments are packed in to xp, 0004 % and so that if presented with a matrix x, it produces a return value of 0005 % same dimension of x. The number of rows in x and FUN(x) are always the 0006 % same. The number of columns is the number of different input arguments 0007 % at which FUN is to be evaluated. 0008 % 0009 % gradfun: string naming the function called to evaluate the gradient matrix. If this 0010 % is null (i.e. just "[]"), a numerical gradient is used instead. 0011 % crit: if the sum of absolute values that FUN returns is less than this, 0012 % the equation is solved. 0013 % itmax: the solver stops when this number of iterations is reached, with rc=4 0014 % varargin: in this position the user can place any number of additional arguments, all 0015 % of which are passed on to FUN and gradfun (when it is non-empty) as a list of 0016 % arguments following x. 0017 % rc: 0 means normal solution, 1 and 3 mean no solution despite extremely fine adjustments 0018 % in step length (very likely a numerical problem, or a discontinuity). 4 means itmax 0019 % termination. 0020 0021 % Original file downloaded from: 0022 % http://sims.princeton.edu/yftp/optimize/mfiles/csolve.m 0023 0024 % Copyright (C) 1993-2007 Christopher Sims 0025 % Copyright (C) 2007-2011 Dynare Team 0026 % 0027 % This file is part of Dynare. 0028 % 0029 % Dynare is free software: you can redistribute it and/or modify 0030 % it under the terms of the GNU General Public License as published by 0031 % the Free Software Foundation, either version 3 of the License, or 0032 % (at your option) any later version. 0033 % 0034 % Dynare is distributed in the hope that it will be useful, 0035 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0036 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0037 % GNU General Public License for more details. 0038 % 0039 % You should have received a copy of the GNU General Public License 0040 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0041 0042 %---------- delta -------------------- 0043 % differencing interval for numerical gradient 0044 delta = 1e-6; 0045 %------------------------------------- 0046 %------------ alpha ------------------ 0047 % tolerance on rate of descent 0048 alpha=1e-3; 0049 %--------------------------------------- 0050 %------------ verbose ---------------- 0051 verbose=0;% if this is set to zero, all screen output is suppressed 0052 %------------------------------------- 0053 %------------ analyticg -------------- 0054 analyticg=1-isempty(gradfun); %if the grad argument is [], numerical derivatives are used. 0055 %------------------------------------- 0056 nv=length(x); 0057 tvec=delta*eye(nv); 0058 done=0; 0059 if isempty(varargin) 0060 f0=feval(FUN,x); 0061 else 0062 f0=feval(FUN,x,varargin{:}); 0063 end 0064 af0=sum(abs(f0)); 0065 af00=af0; 0066 itct=0; 0067 while ~done 0068 % disp([af00-af0 crit*max(1,af0)]) 0069 if itct>3 && af00-af0<crit*max(1,af0) && rem(itct,2)==1 0070 randomize=1; 0071 else 0072 if ~analyticg 0073 % $$$ if isempty(varargin) 0074 % $$$ grad = (feval(FUN,x*ones(1,nv)+tvec)-f0*ones(1,nv))/delta; 0075 % $$$ else 0076 % $$$ grad = (feval(FUN,x*ones(1,nv)+tvec,varargin{:})-f0*ones(1,nv))/delta; 0077 % $$$ end 0078 grad = zeros(nv,nv); 0079 for i=1:nv 0080 grad(:,i) = (feval(FUN,x+tvec(:,i),varargin{:})-f0)/delta; 0081 end 0082 else % use analytic gradient 0083 % grad=feval(gradfun,x,varargin{:}); 0084 [f0,grad] = feval(gradfun,x,varargin{:}); 0085 end 0086 if isreal(grad) 0087 if rcond(grad)<1e-12 0088 grad=grad+tvec; 0089 end 0090 dx0=-grad\f0; 0091 randomize=0; 0092 else 0093 if(verbose),disp('gradient imaginary'),end 0094 randomize=1; 0095 end 0096 end 0097 if randomize 0098 if(verbose),fprintf(1,'\n Random Search'),end 0099 dx0=norm(x)./randn(size(x)); 0100 end 0101 lambda=1; 0102 lambdamin=1; 0103 fmin=f0; 0104 xmin=x; 0105 afmin=af0; 0106 dxSize=norm(dx0); 0107 factor=.6; 0108 shrink=1; 0109 subDone=0; 0110 while ~subDone 0111 dx=lambda*dx0; 0112 f=feval(FUN,x+dx,varargin{:}); 0113 af=sum(abs(f)); 0114 if af<afmin 0115 afmin=af; 0116 fmin=f; 0117 lambdamin=lambda; 0118 xmin=x+dx; 0119 end 0120 if ((lambda >0) && (af0-af < alpha*lambda*af0)) || ((lambda<0) && (af0-af < 0) ) 0121 if ~shrink 0122 factor=factor^.6; 0123 shrink=1; 0124 end 0125 if abs(lambda*(1-factor))*dxSize > .1*delta; 0126 lambda = factor*lambda; 0127 elseif (lambda > 0) && (factor==.6) %i.e., we've only been shrinking 0128 lambda=-.3; 0129 else % 0130 subDone=1; 0131 if lambda > 0 0132 if factor==.6 0133 rc = 2; 0134 else 0135 rc = 1; 0136 end 0137 else 0138 rc=3; 0139 end 0140 end 0141 elseif (lambda >0) && (af-af0 > (1-alpha)*lambda*af0) 0142 if shrink 0143 factor=factor^.6; 0144 shrink=0; 0145 end 0146 lambda=lambda/factor; 0147 else % good value found 0148 subDone=1; 0149 rc=0; 0150 end 0151 end % while ~subDone 0152 itct=itct+1; 0153 if(verbose) 0154 fprintf(1,'\nitct %d, af %g, lambda %g, rc %g',itct,afmin,lambdamin,rc) 0155 fprintf(1,'\n x %10g %10g %10g %10g',xmin); 0156 fprintf(1,'\n f %10g %10g %10g %10g',fmin); 0157 end 0158 x=xmin; 0159 f0=fmin; 0160 af00=af0; 0161 af0=afmin; 0162 if itct >= itmax 0163 done=1; 0164 rc=4; 0165 elseif af0<crit; 0166 done=1; 0167 rc=0; 0168 end 0169 end