← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/optimize/minpack2/dcsrch.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/optimize/minpack2/dcsrch.f
 
1
      subroutine dcsrch(stp,f,g,ftol,gtol,xtol,task,stpmin,stpmax,
 
2
     +                  isave,dsave)
 
3
      character*(*) task
 
4
      integer isave(2)
 
5
      double precision f, g, stp, ftol, gtol, xtol, stpmin, stpmax
 
6
      double precision dsave(13)
 
7
c     **********
 
8
c
 
9
c     Subroutine dcsrch
 
10
c
 
11
c     This subroutine finds a step that satisfies a sufficient
 
12
c     decrease condition and a curvature condition.
 
13
c
 
14
c     Each call of the subroutine updates an interval with
 
15
c     endpoints stx and sty. The interval is initially chosen
 
16
c     so that it contains a minimizer of the modified function
 
17
c
 
18
c           psi(stp) = f(stp) - f(0) - ftol*stp*f'(0).
 
19
c
 
20
c     If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
 
21
c     interval is chosen so that it contains a minimizer of f.
 
22
c
 
23
c     The algorithm is designed to find a step that satisfies
 
24
c     the sufficient decrease condition
 
25
c
 
26
c           f(stp) <= f(0) + ftol*stp*f'(0),
 
27
c
 
28
c     and the curvature condition
 
29
c
 
30
c           abs(f'(stp)) <= gtol*abs(f'(0)).
 
31
c
 
32
c     If ftol is less than gtol and if, for example, the function
 
33
c     is bounded below, then there is always a step which satisfies
 
34
c     both conditions.
 
35
c
 
36
c     If no step can be found that satisfies both conditions, then
 
37
c     the algorithm stops with a warning. In this case stp only
 
38
c     satisfies the sufficient decrease condition.
 
39
c
 
40
c     A typical invocation of dcsrch has the following outline:
 
41
c
 
42
c     Evaluate the function at stp = 0.0d0; store in f.
 
43
c     Evaluate the gradient at stp = 0.0d0; store in g.
 
44
c     Choose a starting step stp.
 
45
c
 
46
c     task = 'START'
 
47
c  10 continue
 
48
c        call dcsrch(stp,f,g,ftol,gtol,xtol,task,stpmin,stpmax,
 
49
c    +               isave,dsave)
 
50
c        if (task .eq. 'FG') then
 
51
c           Evaluate the function and the gradient at stp
 
52
c           go to 10
 
53
c           end if
 
54
c
 
55
c     NOTE: The user must not alter work arrays between calls.
 
56
c
 
57
c     The subroutine statement is
 
58
c
 
59
c       subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
 
60
c                         task,isave,dsave)
 
61
c     where
 
62
c
 
63
c       stp is a double precision variable.
 
64
c         On entry stp is the current estimate of a satisfactory
 
65
c            step. On initial entry, a positive initial estimate
 
66
c            must be provided.
 
67
c         On exit stp is the current estimate of a satisfactory step
 
68
c            if task = 'FG'. If task = 'CONV' then stp satisfies
 
69
c            the sufficient decrease and curvature condition.
 
70
c
 
71
c       f is a double precision variable.
 
72
c         On initial entry f is the value of the function at 0.
 
73
c            On subsequent entries f is the value of the
 
74
c            function at stp.
 
75
c         On exit f is the value of the function at stp.
 
76
c
 
77
c       g is a double precision variable.
 
78
c         On initial entry g is the derivative of the function at 0.
 
79
c            On subsequent entries g is the derivative of the
 
80
c            function at stp.
 
81
c         On exit g is the derivative of the function at stp.
 
82
c
 
83
c       ftol is a double precision variable.
 
84
c         On entry ftol specifies a nonnegative tolerance for the
 
85
c            sufficient decrease condition.
 
86
c         On exit ftol is unchanged.
 
87
c
 
88
c       gtol is a double precision variable.
 
89
c         On entry gtol specifies a nonnegative tolerance for the
 
90
c            curvature condition.
 
91
c         On exit gtol is unchanged.
 
92
c
 
93
c       xtol is a double precision variable.
 
94
c         On entry xtol specifies a nonnegative relative tolerance
 
95
c            for an acceptable step. The subroutine exits with a
 
96
c            warning if the relative difference between sty and stx
 
97
c            is less than xtol.
 
98
c         On exit xtol is unchanged.
 
99
c
 
100
c       task is a character variable of length at least 60.
 
101
c         On initial entry task must be set to 'START'.
 
102
c         On exit task indicates the required action:
 
103
c
 
104
c            If task(1:2) = 'FG' then evaluate the function and
 
105
c            derivative at stp and call dcsrch again.
 
106
c
 
107
c            If task(1:4) = 'CONV' then the search is successful.
 
108
c
 
109
c            If task(1:4) = 'WARN' then the subroutine is not able
 
110
c            to satisfy the convergence conditions. The exit value of
 
111
c            stp contains the best point found during the search.
 
112
c
 
113
c            If task(1:5) = 'ERROR' then there is an error in the
 
114
c            input arguments.
 
115
c
 
116
c         On exit with convergence, a warning or an error, the
 
117
c            variable task contains additional information.
 
118
c
 
119
c       stpmin is a double precision variable.
 
120
c         On entry stpmin is a nonnegative lower bound for the step.
 
121
c         On exit stpmin is unchanged.
 
122
c
 
123
c       stpmax is a double precision variable.
 
124
c         On entry stpmax is a nonnegative upper bound for the step.
 
125
c         On exit stpmax is unchanged.
 
126
c
 
127
c       isave is an integer work array of dimension 2.
 
128
c
 
129
c       dsave is a double precision work array of dimension 13.
 
130
c
 
131
c     Subprograms called
 
132
c
 
133
c       MINPACK-2 ... dcstep
 
134
c
 
135
c     MINPACK-1 Project. June 1983.
 
136
c     Argonne National Laboratory.
 
137
c     Jorge J. More' and David J. Thuente.
 
138
c
 
139
c     MINPACK-2 Project. November 1993.
 
140
c     Argonne National Laboratory and University of Minnesota.
 
141
c     Brett M. Averick, Richard G. Carter, and Jorge J. More'.
 
142
c
 
143
c     **********
 
144
      double precision zero, p5, p66
 
145
      parameter (zero=0.0d0,p5=0.5d0,p66=0.66d0)
 
146
      double precision xtrapl, xtrapu
 
147
      parameter (xtrapl=1.1d0,xtrapu=4.0d0)
 
148
 
 
149
      logical brackt
 
150
      integer stage
 
151
      double precision finit, ftest, fm, fx, fxm, fy, fym, ginit, gtest,
 
152
     +                 gm, gx, gxm, gy, gym, stx, sty, stmin, stmax,
 
153
     +                 width, width1
 
154
 
 
155
      external dcstep
 
156
 
 
157
c     Initialization block.
 
158
 
 
159
      if (task(1:5) .eq. 'START') then
 
160
 
 
161
c        Check the input arguments for errors.
 
162
 
 
163
         if (stp .lt. stpmin) task = 'ERROR: STP .LT. STPMIN'
 
164
         if (stp .gt. stpmax) task = 'ERROR: STP .GT. STPMAX'
 
165
         if (g .ge. zero) task = 'ERROR: INITIAL G .GE. ZERO'
 
166
         if (ftol .lt. zero) task = 'ERROR: FTOL .LT. ZERO'
 
167
         if (gtol .lt. zero) task = 'ERROR: GTOL .LT. ZERO'
 
168
         if (xtol .lt. zero) task = 'ERROR: XTOL .LT. ZERO'
 
169
         if (stpmin .lt. zero) task = 'ERROR: STPMIN .LT. ZERO'
 
170
         if (stpmax .lt. stpmin) task = 'ERROR: STPMAX .LT. STPMIN'
 
171
 
 
172
c        Exit if there are errors on input.
 
173
 
 
174
         if (task(1:5) .eq. 'ERROR') return
 
175
 
 
176
c        Initialize local variables.
 
177
 
 
178
         brackt = .false.
 
179
         stage = 1
 
180
         finit = f
 
181
         ginit = g
 
182
         gtest = ftol*ginit
 
183
         width = stpmax - stpmin
 
184
         width1 = width/p5
 
185
 
 
186
c        The variables stx, fx, gx contain the values of the step,
 
187
c        function, and derivative at the best step.
 
188
c        The variables sty, fy, gy contain the value of the step,
 
189
c        function, and derivative at sty.
 
190
c        The variables stp, f, g contain the values of the step,
 
191
c        function, and derivative at stp.
 
192
 
 
193
         stx = zero
 
194
         fx = finit
 
195
         gx = ginit
 
196
         sty = zero
 
197
         fy = finit
 
198
         gy = ginit
 
199
         stmin = zero
 
200
         stmax = stp + xtrapu*stp
 
201
         task = 'FG'
 
202
 
 
203
         go to 10
 
204
 
 
205
      else
 
206
 
 
207
c        Restore local variables.
 
208
 
 
209
         if (isave(1) .eq. 1) then
 
210
            brackt = .true.
 
211
         else
 
212
            brackt = .false.
 
213
         end if
 
214
         stage = isave(2)
 
215
         ginit = dsave(1)
 
216
         gtest = dsave(2)
 
217
         gx = dsave(3)
 
218
         gy = dsave(4)
 
219
         finit = dsave(5)
 
220
         fx = dsave(6)
 
221
         fy = dsave(7)
 
222
         stx = dsave(8)
 
223
         sty = dsave(9)
 
224
         stmin = dsave(10)
 
225
         stmax = dsave(11)
 
226
         width = dsave(12)
 
227
         width1 = dsave(13)
 
228
 
 
229
      end if
 
230
 
 
231
c     If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
 
232
c     algorithm enters the second stage.
 
233
 
 
234
      ftest = finit + stp*gtest
 
235
      if (stage .eq. 1 .and. f .le. ftest .and. g .ge. zero) stage = 2
 
236
 
 
237
c     Test for warnings.
 
238
 
 
239
      if (brackt .and. (stp .le. stmin .or. stp .ge. stmax))
 
240
     +    task = 'WARNING: ROUNDING ERRORS PREVENT PROGRESS'
 
241
      if (brackt .and. stmax-stmin .le. xtol*stmax)
 
242
     +    task = 'WARNING: XTOL TEST SATISFIED'
 
243
      if (stp .eq. stpmax .and. f .le. ftest .and. g .le. gtest)
 
244
     +    task = 'WARNING: STP = STPMAX'
 
245
      if (stp .eq. stpmin .and. (f .gt. ftest .or. g .ge. gtest))
 
246
     +    task = 'WARNING: STP = STPMIN'
 
247
 
 
248
c     Test for convergence.
 
249
 
 
250
      if (f .le. ftest .and. abs(g) .le. gtol*(-ginit))
 
251
     +    task = 'CONVERGENCE'
 
252
 
 
253
c     Test for termination.
 
254
 
 
255
      if (task(1:4) .eq. 'WARN' .or. task(1:4) .eq. 'CONV') go to 10
 
256
 
 
257
c     A modified function is used to predict the step during the
 
258
c     first stage if a lower function value has been obtained but
 
259
c     the decrease is not sufficient.
 
260
 
 
261
      if (stage .eq. 1 .and. f .le. fx .and. f .gt. ftest) then
 
262
 
 
263
c        Define the modified function and derivative values.
 
264
 
 
265
         fm = f - stp*gtest
 
266
         fxm = fx - stx*gtest
 
267
         fym = fy - sty*gtest
 
268
         gm = g - gtest
 
269
         gxm = gx - gtest
 
270
         gym = gy - gtest
 
271
 
 
272
c        Call dcstep to update stx, sty, and to compute the new step.
 
273
 
 
274
         call dcstep(stx,fxm,gxm,sty,fym,gym,stp,fm,gm,brackt,stmin,
 
275
     +               stmax)
 
276
 
 
277
c        Reset the function and derivative values for f.
 
278
 
 
279
         fx = fxm + stx*gtest
 
280
         fy = fym + sty*gtest
 
281
         gx = gxm + gtest
 
282
         gy = gym + gtest
 
283
 
 
284
      else
 
285
 
 
286
c       Call dcstep to update stx, sty, and to compute the new step.
 
287
 
 
288
         call dcstep(stx,fx,gx,sty,fy,gy,stp,f,g,brackt,stmin,stmax)
 
289
 
 
290
      end if
 
291
 
 
292
c     Decide if a bisection step is needed.
 
293
 
 
294
      if (brackt) then
 
295
         if (abs(sty-stx) .ge. p66*width1) stp = stx + p5*(sty-stx)
 
296
         width1 = width
 
297
         width = abs(sty-stx)
 
298
      end if
 
299
 
 
300
c     Set the minimum and maximum steps allowed for stp.
 
301
 
 
302
      if (brackt) then
 
303
         stmin = min(stx,sty)
 
304
         stmax = max(stx,sty)
 
305
      else
 
306
         stmin = stp + xtrapl*(stp-stx)
 
307
         stmax = stp + xtrapu*(stp-stx)
 
308
      end if
 
309
 
 
310
c     Force the step to be within the bounds stpmax and stpmin.
 
311
 
 
312
      stp = max(stp,stpmin)
 
313
      stp = min(stp,stpmax)
 
314
 
 
315
c     If further progress is not possible, let stp be the best
 
316
c     point obtained during the search.
 
317
 
 
318
      if (brackt .and. (stp .le. stmin .or. stp .ge. stmax) .or.
 
319
     +    (brackt .and. stmax-stmin .le. xtol*stmax)) stp = stx
 
320
 
 
321
c     Obtain another function and derivative.
 
322
 
 
323
      task = 'FG'
 
324
 
 
325
   10 continue
 
326
 
 
327
c     Save local variables.
 
328
 
 
329
      if (brackt) then
 
330
         isave(1) = 1
 
331
      else
 
332
         isave(1) = 0
 
333
      end if
 
334
      isave(2) = stage
 
335
      dsave(1) = ginit
 
336
      dsave(2) = gtest
 
337
      dsave(3) = gx
 
338
      dsave(4) = gy
 
339
      dsave(5) = finit
 
340
      dsave(6) = fx
 
341
      dsave(7) = fy
 
342
      dsave(8) = stx
 
343
      dsave(9) = sty
 
344
      dsave(10) = stmin
 
345
      dsave(11) = stmax
 
346
      dsave(12) = width
 
347
      dsave(13) = width1
 
348
 
 
349
      end