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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
|
subroutine dcsrch(stp,f,g,ftol,gtol,xtol,task,stpmin,stpmax,
+ isave,dsave)
character*(*) task
integer isave(2)
double precision f, g, stp, ftol, gtol, xtol, stpmin, stpmax
double precision dsave(13)
c **********
c
c Subroutine dcsrch
c
c This subroutine finds a step that satisfies a sufficient
c decrease condition and a curvature condition.
c
c Each call of the subroutine updates an interval with
c endpoints stx and sty. The interval is initially chosen
c so that it contains a minimizer of the modified function
c
c psi(stp) = f(stp) - f(0) - ftol*stp*f'(0).
c
c If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
c interval is chosen so that it contains a minimizer of f.
c
c The algorithm is designed to find a step that satisfies
c the sufficient decrease condition
c
c f(stp) <= f(0) + ftol*stp*f'(0),
c
c and the curvature condition
c
c abs(f'(stp)) <= gtol*abs(f'(0)).
c
c If ftol is less than gtol and if, for example, the function
c is bounded below, then there is always a step which satisfies
c both conditions.
c
c If no step can be found that satisfies both conditions, then
c the algorithm stops with a warning. In this case stp only
c satisfies the sufficient decrease condition.
c
c A typical invocation of dcsrch has the following outline:
c
c Evaluate the function at stp = 0.0d0; store in f.
c Evaluate the gradient at stp = 0.0d0; store in g.
c Choose a starting step stp.
c
c task = 'START'
c 10 continue
c call dcsrch(stp,f,g,ftol,gtol,xtol,task,stpmin,stpmax,
c + isave,dsave)
c if (task .eq. 'FG') then
c Evaluate the function and the gradient at stp
c go to 10
c end if
c
c NOTE: The user must not alter work arrays between calls.
c
c The subroutine statement is
c
c subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
c task,isave,dsave)
c where
c
c stp is a double precision variable.
c On entry stp is the current estimate of a satisfactory
c step. On initial entry, a positive initial estimate
c must be provided.
c On exit stp is the current estimate of a satisfactory step
c if task = 'FG'. If task = 'CONV' then stp satisfies
c the sufficient decrease and curvature condition.
c
c f is a double precision variable.
c On initial entry f is the value of the function at 0.
c On subsequent entries f is the value of the
c function at stp.
c On exit f is the value of the function at stp.
c
c g is a double precision variable.
c On initial entry g is the derivative of the function at 0.
c On subsequent entries g is the derivative of the
c function at stp.
c On exit g is the derivative of the function at stp.
c
c ftol is a double precision variable.
c On entry ftol specifies a nonnegative tolerance for the
c sufficient decrease condition.
c On exit ftol is unchanged.
c
c gtol is a double precision variable.
c On entry gtol specifies a nonnegative tolerance for the
c curvature condition.
c On exit gtol is unchanged.
c
c xtol is a double precision variable.
c On entry xtol specifies a nonnegative relative tolerance
c for an acceptable step. The subroutine exits with a
c warning if the relative difference between sty and stx
c is less than xtol.
c On exit xtol is unchanged.
c
c task is a character variable of length at least 60.
c On initial entry task must be set to 'START'.
c On exit task indicates the required action:
c
c If task(1:2) = 'FG' then evaluate the function and
c derivative at stp and call dcsrch again.
c
c If task(1:4) = 'CONV' then the search is successful.
c
c If task(1:4) = 'WARN' then the subroutine is not able
c to satisfy the convergence conditions. The exit value of
c stp contains the best point found during the search.
c
c If task(1:5) = 'ERROR' then there is an error in the
c input arguments.
c
c On exit with convergence, a warning or an error, the
c variable task contains additional information.
c
c stpmin is a double precision variable.
c On entry stpmin is a nonnegative lower bound for the step.
c On exit stpmin is unchanged.
c
c stpmax is a double precision variable.
c On entry stpmax is a nonnegative upper bound for the step.
c On exit stpmax is unchanged.
c
c isave is an integer work array of dimension 2.
c
c dsave is a double precision work array of dimension 13.
c
c Subprograms called
c
c MINPACK-2 ... dcstep
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, Richard G. Carter, and Jorge J. More'.
c
c **********
double precision zero, p5, p66
parameter (zero=0.0d0,p5=0.5d0,p66=0.66d0)
double precision xtrapl, xtrapu
parameter (xtrapl=1.1d0,xtrapu=4.0d0)
logical brackt
integer stage
double precision finit, ftest, fm, fx, fxm, fy, fym, ginit, gtest,
+ gm, gx, gxm, gy, gym, stx, sty, stmin, stmax,
+ width, width1
external dcstep
c Initialization block.
if (task(1:5) .eq. 'START') then
c Check the input arguments for errors.
if (stp .lt. stpmin) task = 'ERROR: STP .LT. STPMIN'
if (stp .gt. stpmax) task = 'ERROR: STP .GT. STPMAX'
if (g .ge. zero) task = 'ERROR: INITIAL G .GE. ZERO'
if (ftol .lt. zero) task = 'ERROR: FTOL .LT. ZERO'
if (gtol .lt. zero) task = 'ERROR: GTOL .LT. ZERO'
if (xtol .lt. zero) task = 'ERROR: XTOL .LT. ZERO'
if (stpmin .lt. zero) task = 'ERROR: STPMIN .LT. ZERO'
if (stpmax .lt. stpmin) task = 'ERROR: STPMAX .LT. STPMIN'
c Exit if there are errors on input.
if (task(1:5) .eq. 'ERROR') return
c Initialize local variables.
brackt = .false.
stage = 1
finit = f
ginit = g
gtest = ftol*ginit
width = stpmax - stpmin
width1 = width/p5
c The variables stx, fx, gx contain the values of the step,
c function, and derivative at the best step.
c The variables sty, fy, gy contain the value of the step,
c function, and derivative at sty.
c The variables stp, f, g contain the values of the step,
c function, and derivative at stp.
stx = zero
fx = finit
gx = ginit
sty = zero
fy = finit
gy = ginit
stmin = zero
stmax = stp + xtrapu*stp
task = 'FG'
go to 10
else
c Restore local variables.
if (isave(1) .eq. 1) then
brackt = .true.
else
brackt = .false.
end if
stage = isave(2)
ginit = dsave(1)
gtest = dsave(2)
gx = dsave(3)
gy = dsave(4)
finit = dsave(5)
fx = dsave(6)
fy = dsave(7)
stx = dsave(8)
sty = dsave(9)
stmin = dsave(10)
stmax = dsave(11)
width = dsave(12)
width1 = dsave(13)
end if
c If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
c algorithm enters the second stage.
ftest = finit + stp*gtest
if (stage .eq. 1 .and. f .le. ftest .and. g .ge. zero) stage = 2
c Test for warnings.
if (brackt .and. (stp .le. stmin .or. stp .ge. stmax))
+ task = 'WARNING: ROUNDING ERRORS PREVENT PROGRESS'
if (brackt .and. stmax-stmin .le. xtol*stmax)
+ task = 'WARNING: XTOL TEST SATISFIED'
if (stp .eq. stpmax .and. f .le. ftest .and. g .le. gtest)
+ task = 'WARNING: STP = STPMAX'
if (stp .eq. stpmin .and. (f .gt. ftest .or. g .ge. gtest))
+ task = 'WARNING: STP = STPMIN'
c Test for convergence.
if (f .le. ftest .and. abs(g) .le. gtol*(-ginit))
+ task = 'CONVERGENCE'
c Test for termination.
if (task(1:4) .eq. 'WARN' .or. task(1:4) .eq. 'CONV') go to 10
c A modified function is used to predict the step during the
c first stage if a lower function value has been obtained but
c the decrease is not sufficient.
if (stage .eq. 1 .and. f .le. fx .and. f .gt. ftest) then
c Define the modified function and derivative values.
fm = f - stp*gtest
fxm = fx - stx*gtest
fym = fy - sty*gtest
gm = g - gtest
gxm = gx - gtest
gym = gy - gtest
c Call dcstep to update stx, sty, and to compute the new step.
call dcstep(stx,fxm,gxm,sty,fym,gym,stp,fm,gm,brackt,stmin,
+ stmax)
c Reset the function and derivative values for f.
fx = fxm + stx*gtest
fy = fym + sty*gtest
gx = gxm + gtest
gy = gym + gtest
else
c Call dcstep to update stx, sty, and to compute the new step.
call dcstep(stx,fx,gx,sty,fy,gy,stp,f,g,brackt,stmin,stmax)
end if
c Decide if a bisection step is needed.
if (brackt) then
if (abs(sty-stx) .ge. p66*width1) stp = stx + p5*(sty-stx)
width1 = width
width = abs(sty-stx)
end if
c Set the minimum and maximum steps allowed for stp.
if (brackt) then
stmin = min(stx,sty)
stmax = max(stx,sty)
else
stmin = stp + xtrapl*(stp-stx)
stmax = stp + xtrapu*(stp-stx)
end if
c Force the step to be within the bounds stpmax and stpmin.
stp = max(stp,stpmin)
stp = min(stp,stpmax)
c If further progress is not possible, let stp be the best
c point obtained during the search.
if (brackt .and. (stp .le. stmin .or. stp .ge. stmax) .or.
+ (brackt .and. stmax-stmin .le. xtol*stmax)) stp = stx
c Obtain another function and derivative.
task = 'FG'
10 continue
c Save local variables.
if (brackt) then
isave(1) = 1
else
isave(1) = 0
end if
isave(2) = stage
dsave(1) = ginit
dsave(2) = gtest
dsave(3) = gx
dsave(4) = gy
dsave(5) = finit
dsave(6) = fx
dsave(7) = fy
dsave(8) = stx
dsave(9) = sty
dsave(10) = stmin
dsave(11) = stmax
dsave(12) = width
dsave(13) = width1
end
|
