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