1function [x,rc] = csolve(FUN,x,gradfun,crit,itmax,varargin) 2%function [x,rc] = csolve(FUN,x,gradfun,crit,itmax,varargin) 3% FUN should be written so that any parametric arguments are packed in to xp, 4% and so that if presented with a matrix x, it produces a return value of 5% same dimension of x. The number of rows in x and FUN(x) are always the 6% same. The number of columns is the number of different input arguments 7% at which FUN is to be evaluated. 8% 9% gradfun: string naming the function called to evaluate the gradient matrix. If this 10% is null (i.e. just "[]"), a numerical gradient is used instead. 11% crit: if the sum of absolute values that FUN returns is less than this, 12% the equation is solved. 13% itmax: the solver stops when this number of iterations is reached, with rc=4 14% varargin: in this position the user can place any number of additional arguments, all 15% of which are passed on to FUN and gradfun (when it is non-empty) as a list of 16% arguments following x. 17% rc: 0 means normal solution, 1 and 3 mean no solution despite extremely fine adjustments 18% in step length (very likely a numerical problem, or a discontinuity). 4 means itmax 19% termination. 20 21% Original file downloaded from: 22% http://sims.princeton.edu/yftp/optimize/mfiles/csolve.m 23 24% Copyright (C) 1993-2007 Christopher Sims 25% Copyright (C) 2007-2017 Dynare Team 26% 27% This file is part of Dynare. 28% 29% Dynare is free software: you can redistribute it and/or modify 30% it under the terms of the GNU General Public License as published by 31% the Free Software Foundation, either version 3 of the License, or 32% (at your option) any later version. 33% 34% Dynare is distributed in the hope that it will be useful, 35% but WITHOUT ANY WARRANTY; without even the implied warranty of 36% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 37% GNU General Public License for more details. 38% 39% You should have received a copy of the GNU General Public License 40% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 41 42%---------- delta -------------------- 43% differencing interval for numerical gradient 44delta = 1e-6; 45%------------------------------------- 46%------------ alpha ------------------ 47% tolerance on rate of descent 48alpha=1e-3; 49%--------------------------------------- 50%------------ verbose ---------------- 51verbose=0;% if this is set to zero, all screen output is suppressed 52 %------------------------------------- 53 %------------ analyticg -------------- 54analyticg=1-isempty(gradfun); %if the grad argument is [], numerical derivatives are used. 55 %------------------------------------- 56nv=length(x); 57tvec=delta*eye(nv); 58done=0; 59if isempty(varargin) 60 f0=feval(FUN,x); 61else 62 f0=feval(FUN,x,varargin{:}); 63end 64af0=sum(abs(f0)); 65af00=af0; 66itct=0; 67while ~done 68 % disp([af00-af0 crit*max(1,af0)]) 69 if itct>3 && af00-af0<crit*max(1,af0) && rem(itct,2)==1 70 randomize=1; 71 else 72 if ~analyticg 73% $$$ if isempty(varargin) 74% $$$ grad = (feval(FUN,x*ones(1,nv)+tvec)-f0*ones(1,nv))/delta; 75% $$$ else 76% $$$ grad = (feval(FUN,x*ones(1,nv)+tvec,varargin{:})-f0*ones(1,nv))/delta; 77% $$$ end 78 grad = zeros(nv,nv); 79 for i=1:nv 80 grad(:,i) = (feval(FUN,x+tvec(:,i),varargin{:})-f0)/delta; 81 end 82 else % use analytic gradient 83 % grad=feval(gradfun,x,varargin{:}); 84 [f0,grad] = feval(gradfun,x,varargin{:}); 85 end 86 if isreal(grad) 87 if rcond(full(grad))<1e-12 88 grad=grad+tvec; 89 end 90 dx0=-grad\f0; 91 randomize=0; 92 else 93 if(verbose),disp('gradient imaginary'),end 94 randomize=1; 95 end 96 end 97 if randomize 98 if(verbose),fprintf(1,'\n Random Search'),end 99 dx0=norm(x)./randn(size(x)); 100 end 101 lambda=1; 102 lambdamin=1; 103 fmin=f0; 104 xmin=x; 105 afmin=af0; 106 dxSize=norm(dx0); 107 factor=.6; 108 shrink=1; 109 subDone=0; 110 while ~subDone 111 dx=lambda*dx0; 112 f=feval(FUN,x+dx,varargin{:}); 113 af=sum(abs(f)); 114 if af<afmin 115 afmin=af; 116 fmin=f; 117 lambdamin=lambda; 118 xmin=x+dx; 119 end 120 if ((lambda >0) && (af0-af < alpha*lambda*af0)) || ((lambda<0) && (af0-af < 0) ) 121 if ~shrink 122 factor=factor^.6; 123 shrink=1; 124 end 125 if abs(lambda*(1-factor))*dxSize > .1*delta 126 lambda = factor*lambda; 127 elseif (lambda > 0) && (factor==.6) %i.e., we've only been shrinking 128 lambda=-.3; 129 else % 130 subDone=1; 131 if lambda > 0 132 if factor==.6 133 rc = 2; 134 else 135 rc = 1; 136 end 137 else 138 rc=3; 139 end 140 end 141 elseif (lambda >0) && (af-af0 > (1-alpha)*lambda*af0) 142 if shrink 143 factor=factor^.6; 144 shrink=0; 145 end 146 lambda=lambda/factor; 147 else % good value found 148 subDone=1; 149 rc=0; 150 end 151 end % while ~subDone 152 itct=itct+1; 153 if(verbose) 154 fprintf(1,'\nitct %d, af %g, lambda %g, rc %g',itct,afmin,lambdamin,rc) 155 fprintf(1,'\n x %10g %10g %10g %10g',xmin); 156 fprintf(1,'\n f %10g %10g %10g %10g',fmin); 157 end 158 x=xmin; 159 f0=fmin; 160 af00=af0; 161 af0=afmin; 162 if itct >= itmax 163 done=1; 164 rc=4; 165 elseif af0<crit 166 done=1; 167 rc=0; 168 end 169end 170