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
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
|
subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,stpmin,
+ stpmax)
logical brackt
double precision stx, fx, dx, sty, fy, dy, stp, fp, dp, stpmin,
+ stpmax
c **********
c
c Subroutine dcstep
c
c This subroutine computes a safeguarded step for a search
c procedure and updates an interval that contains a step that
c satisfies a sufficient decrease and a curvature condition.
c
c The parameter stx contains the step with the least function
c value. If brackt is set to .true. then a minimizer has
c been bracketed in an interval with endpoints stx and sty.
c The parameter stp contains the current step.
c The subroutine assumes that if brackt is set to .true. then
c
c min(stx,sty) < stp < max(stx,sty),
c
c and that the derivative at stx is negative in the direction
c of the step.
c
c The subroutine statement is
c
c subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,
c stpmin,stpmax)
c
c where
c
c stx is a double precision variable.
c On entry stx is the best step obtained so far and is an
c endpoint of the interval that contains the minimizer.
c On exit stx is the updated best step.
c
c fx is a double precision variable.
c On entry fx is the function at stx.
c On exit fx is the function at stx.
c
c dx is a double precision variable.
c On entry dx is the derivative of the function at
c stx. The derivative must be negative in the direction of
c the step, that is, dx and stp - stx must have opposite
c signs.
c On exit dx is the derivative of the function at stx.
c
c sty is a double precision variable.
c On entry sty is the second endpoint of the interval that
c contains the minimizer.
c On exit sty is the updated endpoint of the interval that
c contains the minimizer.
c
c fy is a double precision variable.
c On entry fy is the function at sty.
c On exit fy is the function at sty.
c
c dy is a double precision variable.
c On entry dy is the derivative of the function at sty.
c On exit dy is the derivative of the function at the exit sty.
c
c stp is a double precision variable.
c On entry stp is the current step. If brackt is set to .true.
c then on input stp must be between stx and sty.
c On exit stp is a new trial step.
c
c fp is a double precision variable.
c On entry fp is the function at stp
c On exit fp is unchanged.
c
c dp is a double precision variable.
c On entry dp is the the derivative of the function at stp.
c On exit dp is unchanged.
c
c brackt is an logical variable.
c On entry brackt specifies if a minimizer has been bracketed.
c Initially brackt must be set to .false.
c On exit brackt specifies if a minimizer has been bracketed.
c When a minimizer is bracketed brackt is set to .true.
c
c stpmin is a double precision variable.
c On entry stpmin is a lower bound for the step.
c On exit stpmin is unchanged.
c
c stpmax is a double precision variable.
c On entry stpmax is an upper bound for the step.
c On exit stpmax is unchanged.
c
c MINPACK-1 Project. June 1983
c Argonne National Laboratory.
c Jorge J. More' and David J. Thuente.
c
c MINPACK-2 Project. November 1993.
c Argonne National Laboratory and University of Minnesota.
c Brett M. Averick and Jorge J. More'.
c
c **********
double precision zero, p66, two, three
parameter (zero=0.0d0,p66=0.66d0,two=2.0d0,three=3.0d0)
double precision gamma, p, q, r, s, sgnd, stpc, stpf, stpq, theta
sgnd = dp*(dx/abs(dx))
c First case: A higher function value. The minimum is bracketed.
c If the cubic step is closer to stx than the quadratic step, the
c cubic step is taken, otherwise the average of the cubic and
c quadratic steps is taken.
if (fp .gt. fx) then
theta = three*(fx-fp)/(stp-stx) + dx + dp
s = max(abs(theta),abs(dx),abs(dp))
gamma = s*sqrt((theta/s)**2-(dx/s)*(dp/s))
if (stp .lt. stx) gamma = -gamma
p = (gamma-dx) + theta
q = ((gamma-dx)+gamma) + dp
r = p/q
stpc = stx + r*(stp-stx)
stpq = stx + ((dx/((fx-fp)/(stp-stx)+dx))/two)*(stp-stx)
if (abs(stpc-stx) .lt. abs(stpq-stx)) then
stpf = stpc
else
stpf = stpc + (stpq-stpc)/two
end if
brackt = .true.
c Second case: A lower function value and derivatives of opposite
c sign. The minimum is bracketed. If the cubic step is farther from
c stp than the secant step, the cubic step is taken, otherwise the
c secant step is taken.
else if (sgnd .lt. zero) then
theta = three*(fx-fp)/(stp-stx) + dx + dp
s = max(abs(theta),abs(dx),abs(dp))
gamma = s*sqrt((theta/s)**2-(dx/s)*(dp/s))
if (stp .gt. stx) gamma = -gamma
p = (gamma-dp) + theta
q = ((gamma-dp)+gamma) + dx
r = p/q
stpc = stp + r*(stx-stp)
stpq = stp + (dp/(dp-dx))*(stx-stp)
if (abs(stpc-stp) .gt. abs(stpq-stp)) then
stpf = stpc
else
stpf = stpq
end if
brackt = .true.
c Third case: A lower function value, derivatives of the same sign,
c and the magnitude of the derivative decreases.
else if (abs(dp) .lt. abs(dx)) then
c The cubic step is computed only if the cubic tends to infinity
c in the direction of the step or if the minimum of the cubic
c is beyond stp. Otherwise the cubic step is defined to be the
c secant step.
theta = three*(fx-fp)/(stp-stx) + dx + dp
s = max(abs(theta),abs(dx),abs(dp))
c The case gamma = 0 only arises if the cubic does not tend
c to infinity in the direction of the step.
gamma = s*sqrt(max(zero,(theta/s)**2-(dx/s)*(dp/s)))
if (stp .gt. stx) gamma = -gamma
p = (gamma-dp) + theta
q = (gamma+(dx-dp)) + gamma
r = p/q
if (r .lt. zero .and. gamma .ne. zero) then
stpc = stp + r*(stx-stp)
else if (stp .gt. stx) then
stpc = stpmax
else
stpc = stpmin
end if
stpq = stp + (dp/(dp-dx))*(stx-stp)
if (brackt) then
c A minimizer has been bracketed. If the cubic step is
c closer to stp than the secant step, the cubic step is
c taken, otherwise the secant step is taken.
if (abs(stpc-stp) .lt. abs(stpq-stp)) then
stpf = stpc
else
stpf = stpq
end if
if (stp .gt. stx) then
stpf = min(stp+p66*(sty-stp),stpf)
else
stpf = max(stp+p66*(sty-stp),stpf)
end if
else
c A minimizer has not been bracketed. If the cubic step is
c farther from stp than the secant step, the cubic step is
c taken, otherwise the secant step is taken.
if (abs(stpc-stp) .gt. abs(stpq-stp)) then
stpf = stpc
else
stpf = stpq
end if
stpf = min(stpmax,stpf)
stpf = max(stpmin,stpf)
end if
c Fourth case: A lower function value, derivatives of the same sign,
c and the magnitude of the derivative does not decrease. If the
c minimum is not bracketed, the step is either stpmin or stpmax,
c otherwise the cubic step is taken.
else
if (brackt) then
theta = three*(fp-fy)/(sty-stp) + dy + dp
s = max(abs(theta),abs(dy),abs(dp))
gamma = s*sqrt((theta/s)**2-(dy/s)*(dp/s))
if (stp .gt. sty) gamma = -gamma
p = (gamma-dp) + theta
q = ((gamma-dp)+gamma) + dy
r = p/q
stpc = stp + r*(sty-stp)
stpf = stpc
else if (stp .gt. stx) then
stpf = stpmax
else
stpf = stpmin
end if
end if
c Update the interval which contains a minimizer.
if (fp .gt. fx) then
sty = stp
fy = fp
dy = dp
else
if (sgnd .lt. zero) then
sty = stx
fy = fx
dy = dx
end if
stx = stp
fx = fp
dx = dp
end if
c Compute the new step.
stp = stpf
end
|
