3 c ###################################################
4 c ## COPYRIGHT (C) 1990 by Jay William Ponder ##
5 c ## All Rights Reserved ##
6 c ###################################################
8 c #################################################################
10 c ## subroutine search -- perform unidimensional line search ##
12 c #################################################################
15 c "search" is a unidimensional line search based upon parabolic
16 c extrapolation and cubic interpolation using both function and
19 c variables used by the routine :
21 c f function value at the best line search point
22 c x current values of variables during line search
23 c g gradient at the current point during line search
24 c p initial search vector, unchanged by this routine
25 c s scaled search vector at current line search point
26 c angle angle between search and negative gradient vector
28 c parameters used by the routine :
30 c stpmin minimum step length in current line search direction
31 c stpmax maximum step length in current line search direction
32 c cappa stringency of line search (0=tight < cappa < 1=loose)
33 c slpmax projected gradient above which stepsize is reduced
34 c angmax maximum angle between search direction and -gradient
35 c intmax maximum number of interpolations during line search
37 c status codes upon return :
39 c Success normal termination after satisfying "cappa" test
40 c ScaleStep normal termination after a step size rescaling
41 c ReSearch normal termination after a reinterpolation
42 c WideAngle large angle between search direction and -gradient
43 c BadIntpln unsatisfied "cappa" test after two searches
44 c IntplnErr function value increase or serious gradient error
47 subroutine search (nvar,f,g,x,p,f_move,angle,ncalls,
69 real*8, allocatable :: x_0(:)
70 real*8, allocatable :: s(:)
77 c use default parameters for the line search if needed
80 if (stpmin .eq. 0.0d0) stpmin = 1.0d-16
81 if (stpmax .eq. 0.0d0) stpmax = 2.0d0
82 if (cappa .eq. 0.0d0) cappa = 0.1d0
83 if (slpmax .eq. 0.0d0) slpmax = 10000.0d0
84 if (angmax .eq. 0.0d0) angmax = 180.0d0
85 if (intmax .eq. 0) intmax = 5
87 c perform dynamic allocation of some local arrays
92 c copy the search direction into a new vector
98 c compute the length of gradient and search direction
103 g_norm = g_norm + g(i)*g(i)
104 s_norm = s_norm + s(i)*s(i)
106 g_norm = sqrt(g_norm)
107 s_norm = sqrt(s_norm)
109 c store initial function, then normalize the
110 c search vector and find projected gradient
117 sg_0 = sg_0 + s(i)*g(i)
120 c check the angle between the search direction
121 c and the negative gradient vector
123 cosang = -sg_0 / g_norm
124 cosang = min(1.0d0,max(-1.0d0,cosang))
125 angle = radian * acos(cosang)
126 if (angle .gt. angmax) then
133 c set the initial stepsize to the length of the passed
134 c search vector, or based on previous function decrease
136 step = 2.0d0 * abs(f_move/sg_0)
137 step = min(step,s_norm)
138 if (step .gt. stpmax) step = stpmax
139 if (step .lt. stpmin) step = stpmin
141 c beginning of the parabolic extrapolation procedure
149 c replace last point by latest and take another step
155 x(i) = x(i) + step*s(i)
158 c get new function and projected gradient following a step
164 sg_b = sg_b + s(i)*g(i)
167 c scale stepsize if initial gradient change is too large
169 if (abs(sg_b/sg_a).ge.slpmax .and. restart) then
179 c return if the gradient is small and function decreases
181 if (abs(sg_b/sg_0).le.cappa .and. f_b.lt.f_a) then
183 if (status .eq. blank) status = ' Success '
189 c interpolate if gradient changes sign or function increases
191 if (sg_b*sg_a.lt.0.0d0 .or. f_b.gt.f_a) goto 30
193 c if the finite difference curvature is negative double the step;
194 c or if (step < parabolic estimate < 4*step) use this estimate,
195 c otherwise truncate to step or 4*step, respectively
198 if (sg_b .gt. sg_a) then
199 parab = (f_a-f_b) / (sg_b-sg_a)
200 if (parab .gt. 2.0d0*step) parab = 2.0d0 * step
201 if (parab .lt. 0.5d0*step) parab = 0.5d0 * step
204 if (step .gt. stpmax) step = stpmax
207 c beginning of the cubic interpolation procedure
211 sss = 3.0d0*(f_b-f_a)/step - sg_a - sg_b
212 ttt = sss*sss - sg_a*sg_b
213 if (ttt .lt. 0.0d0) then
221 cube = step * (sg_b+ttt+sss)/(sg_b-sg_a+2.0d0*ttt)
222 if (cube.lt.0.0d0 .or. cube.gt.step) then
230 x(i) = x(i) - cube*s(i)
233 c get new function and gradient, then test for termination
239 sg_c = sg_c + s(i)*g(i)
241 if (abs(sg_c/sg_0) .le. cappa) then
243 if (status .eq. blank) status = ' Success '
249 c get the next pair of bracketing points by replacing one
250 c of the current brackets with the interpolated point
252 if (f_c.le.f_a .or. f_c.le.f_b) then
253 cubstp = min(abs(cube),abs(step-cube))
254 if (cubstp.ge.stpmin .and. intpln.lt.intmax) then
256 c if the current brackets have slopes of opposite sign,
257 c then substitute the interpolated point for the bracket
258 c point with slope of same sign as the interpolated point
260 if (sg_a*sg_b .lt. 0.0d0) then
261 if (sg_a*sg_c .lt. 0.0d0) then
270 x(i) = x(i) + cube*s(i)
274 c if current brackets have slope of same sign, then replace
275 c the far bracket if the interpolated point has a slope of
276 c the opposite sign or a lower function value than the near
277 c bracket, otherwise replace the near bracket point
280 if (sg_a*sg_c.lt.0.0d0 .or. f_a.le.f_c) then
289 x(i) = x(i) + cube*s(i)
297 c interpolation has failed, reset to best current point
299 f_1 = min(f_a,f_b,f_c)
300 if (f_1 .eq. f_a) then
303 x(i) = x(i) + (cube-step)*s(i)
305 else if (f_1 .eq. f_b) then
308 x(i) = x(i) + cube*s(i)
310 else if (f_1 .eq. f_c) then
314 c try to restart from best point with smaller stepsize
316 if (f_1 .gt. f_0) then
326 if (sg_1 .gt. 0.0d0) then
332 step = max(cube,step-cube) / 10.0d0
333 if (step .lt. stpmin) step = stpmin
335 c if already restarted once, then return with best point
337 if (status .eq. ' ReSearch') then