File: solve1.m

package info (click to toggle)
dynare 6.4-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 67,648 kB
  • sloc: cpp: 79,109; ansic: 28,917; objc: 12,430; yacc: 4,528; pascal: 1,993; lex: 1,441; sh: 1,129; python: 634; makefile: 626; lisp: 163; xml: 18
file content (193 lines) | stat: -rw-r--r-- 5,507 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
function [x, errorflag, errorcode] = solve1(func, x, j1, j2, jacobian_flag, gstep, tolf, tolx, maxit, fake, debug, varargin)

% Solves systems of non linear equations of several variables
%
% 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=true: jacobian given by the 'func' function
%    jacobian_flag=false: jacobian obtained numerically
%    gstep            increment multiplier in numercial derivative
%                     computation
%    tolf             tolerance for residuals
%    tolx             tolerance for solution variation
%    maxit            maximum number of iterations
%    fake             unused argument (compatibity with trust_region).
%    debug            debug flag
%    varargin:        list of extra arguments to the function
%
% OUTPUTS
%    x:               results
%    errorflag=1:     the model can not be solved
%
% SPECIAL REQUIREMENTS
%    none

% Copyright © 2001-2022 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 <https://www.gnu.org/licenses/>.

nn = length(j1);
g = zeros(nn,1) ;

tolmin = tolx ;

stpmx = 100 ;

errorflag = false ;

fvec = feval(func,x,varargin{:});
fvec = fvec(j1);

idInf = isinf(fvec);
idNan = isnan(fvec);
idCpx = ~isreal(fvec);

if any(idInf)
    disp('SOLVE1: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number:')
    disp(j1(idInf)')
    errorcode = 0;
    errorflag = true;
    return
end

if any(idNan)
    disp('SOLVE1: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a nan:')
    disp(j1(idNan)')
    errorcode = 0;
    errorflag = true;
    return
end

if any(idNan)
    disp('SOLVE1: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a complex number:')
    disp(j1(idCpx)')
    errorcode = 0;
    errorflag = true;
    return
end

f = 0.5*(fvec'*fvec);

if max(abs(fvec))<tolf*tolf
    % Initial guess is a solution
    errorcode = -1;
    return
end

stpmax = stpmx*max([sqrt(x'*x);nn]) ;
first_time = 1;
if ~jacobian_flag
    fjac = zeros(nn,nn);
end
for its = 1:maxit
    if jacobian_flag
        [fvec,fjac] = feval(func,x, varargin{:});
        fvec = fvec(j1);
        fjac = fjac(j1,j2);
        g = (fvec'*fjac)';
    else
        dh = max(abs(x(j2)),gstep(1)*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
    if debug
        disp(['cond(fjac) ' num2str(condest(fjac))])
    end
    if issparse(fjac)
        rcond_fjac = 1/condest(fjac);
    else
        rcond_fjac = rcond(fjac);
    end
    if rcond_fjac < sqrt(eps)
        fjac2=fjac'*fjac;
        temp=max(sum(abs(fjac2)));
        if temp>0
            p=-(fjac2+sqrt(nn*eps)*temp*eye(nn))\(fjac'*fvec);
        else
            errorflag = true;
            errorcode = 5;
            if nargout<3
                skipline()
                dprintf('SOLVE: Iteration %s', num2str(its))
                disp('Zero Jacobian.')
                skipline()
            end
            return
        end
    else
        p = -fjac\fvec ;
    end
    xold = x ;
    fold = f ;

    [x, f, fvec, lnsearchflag] = lnsrch1(xold, fold, g, p, stpmax, func, j1, j2, tolx, varargin{:});

    if debug
        disp([its f])
        disp([xold x])
    end

    if lnsearchflag
        errorflag = true;
        den = max([f;0.5*nn]) ;
        if max(abs(g).*max([abs(x(j2)') ones(1,nn)])')/den < tolmin
            if max(abs(x(j2)-xold(j2))./max([abs(x(j2)') ones(1,nn)])') < tolx
                errorcode = 3;
                if nargout<3
                    skipline()
                    dprintf('SOLVE: Iteration %s', num2str(its))
                    disp('Convergence on dX.')
                    skipline()
                end
                return
            end
        else
            errorcode = 4;
            if nargout<3
                skipline()
                dprintf('SOLVE: Iteration %s', num2str(its))
                disp('Spurious convergence.')
                disp(x)
            end
            return
        end
    elseif max(abs(fvec)) < tolf
        errorcode = 1;
        return
    end
end

errorflag = true;
errorcode = 2;

if nargout<3
    skipline()
    disp('SOLVE: maxit has been reached')
end

% 01/14/01 MJ lnsearch is now a separate function
% 01/16/01 MJ added varargin to function evaluation
% 04/13/01 MJ added test  f < tolf !!
% 05/11/01 MJ changed tests for 'check' so as to remove 'continue' which is
%             an instruction which appears only in version 6