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
|
function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin)
% function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin)
% Computes the optimal step by minimizing the residual sum of squares
%
% INPUTS
% xold: actual point
% fold: residual sum of squares at the point xold
% g: gradient
% p: Newton direction
% stpmax: maximum step
% func: name of the function
% j1: equations index to be solved
% j2: unknowns index
% varargin: list of arguments following j2
%
% OUTPUTS
% x: chosen point
% f: residual sum of squares value for a given x
% fvec: residuals vector
% check=1: problem of the looping which continues indefinitely
%
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2001-2010 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 options_
alf = 1e-4 ;
tolx = options_.solve_tolx;
alam = 1;
x = xold;
nn = length(j2);
summ = sqrt(sum(p.*p)) ;
if ~isfinite(summ)
error(['Some element of Newton direction isn''t finite. Jacobian maybe' ...
' singular or there is a problem with initial values'])
end
if summ > stpmax
p=p.*stpmax/summ ;
end
slope = g'*p ;
test = max(abs(p)'./max([abs(xold(j2))';ones(1,nn)])) ;
alamin = tolx/test ;
if alamin > 0.1
alamin = 0.1;
end
while 1
if alam < alamin
check = 1 ;
return
end
x(j2) = xold(j2) + (alam*p) ;
fvec = feval(func,x,varargin{:}) ;
fvec = fvec(j1);
f = 0.5*fvec'*fvec ;
if any(isnan(fvec))
alam = alam/2 ;
alam2 = alam ;
f2 = f ;
fold2 = fold ;
else
if f <= fold+alf*alam*slope
check = 0;
break ;
else
if alam == 1
tmplam = -slope/(2*(f-fold-slope)) ;
else
rhs1 = f-fold-alam*slope ;
rhs2 = f2-fold2-alam2*slope ;
a = (rhs1/(alam^2)-rhs2/(alam2^2))/(alam-alam2) ;
b = (-alam2*rhs1/(alam^2)+alam*rhs2/(alam2^2))/(alam-alam2) ;
if a == 0
tmplam = -slope/(2*b) ;
else
disc = (b^2)-3*a*slope ;
if disc < 0
error ('Roundoff problem in nlsearch') ;
else
tmplam = (-b+sqrt(disc))/(3*a) ;
end
end
if tmplam > 0.5*alam
tmplam = 0.5*alam;
end
end
alam2 = alam ;
f2 = f ;
fold2 = fold ;
alam = max([tmplam;(0.1*alam)]) ;
end
end
end
% 01/14/01 MJ lnsearch is now a separate function
% 01/12/03 MJ check for finite summ to avoid infinite loop when Jacobian
% is singular or model is denormalized
|
