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
161 c 3/14/2020 Adam Liwo: added the condition to prevent from infinite
163 if (ncalls.gt.200) then
173 c write (2,*) "ncalls",ncalls," f_a",f_a," f_b",f_b," sg_a",sg_a
176 sg_b = sg_b + s(i)*g(i)
179 c scale stepsize if initial gradient change is too large
181 if (abs(sg_b/sg_a).ge.slpmax .and. restart) then
191 c return if the gradient is small and function decreases
193 if (abs(sg_b/sg_0).le.cappa .and. f_b.lt.f_a) then
195 if (status .eq. blank) status = ' Success '
201 c interpolate if gradient changes sign or function increases
203 if (sg_b*sg_a.lt.0.0d0 .or. f_b.gt.f_a) goto 30
205 c if the finite difference curvature is negative double the step;
206 c or if (step < parabolic estimate < 4*step) use this estimate,
207 c otherwise truncate to step or 4*step, respectively
210 if (sg_b .gt. sg_a) then
211 parab = (f_a-f_b) / (sg_b-sg_a)
212 if (parab .gt. 2.0d0*step) parab = 2.0d0 * step
213 if (parab .lt. 0.5d0*step) parab = 0.5d0 * step
216 if (step .gt. stpmax) step = stpmax
219 c beginning of the cubic interpolation procedure
223 sss = 3.0d0*(f_b-f_a)/step - sg_a - sg_b
224 ttt = sss*sss - sg_a*sg_b
225 if (ttt .lt. 0.0d0) then
233 cube = step * (sg_b+ttt+sss)/(sg_b-sg_a+2.0d0*ttt)
234 if (cube.lt.0.0d0 .or. cube.gt.step) then
242 x(i) = x(i) - cube*s(i)
245 c get new function and gradient, then test for termination
251 sg_c = sg_c + s(i)*g(i)
253 if (abs(sg_c/sg_0) .le. cappa) then
255 if (status .eq. blank) status = ' Success '
261 c get the next pair of bracketing points by replacing one
262 c of the current brackets with the interpolated point
264 if (f_c.le.f_a .or. f_c.le.f_b) then
265 cubstp = min(abs(cube),abs(step-cube))
266 if (cubstp.ge.stpmin .and. intpln.lt.intmax) then
268 c if the current brackets have slopes of opposite sign,
269 c then substitute the interpolated point for the bracket
270 c point with slope of same sign as the interpolated point
272 if (sg_a*sg_b .lt. 0.0d0) then
273 if (sg_a*sg_c .lt. 0.0d0) then
282 x(i) = x(i) + cube*s(i)
286 c if current brackets have slope of same sign, then replace
287 c the far bracket if the interpolated point has a slope of
288 c the opposite sign or a lower function value than the near
289 c bracket, otherwise replace the near bracket point
292 if (sg_a*sg_c.lt.0.0d0 .or. f_a.le.f_c) then
301 x(i) = x(i) + cube*s(i)
309 c interpolation has failed, reset to best current point
311 f_1 = min(f_a,f_b,f_c)
312 if (f_1 .eq. f_a) then
315 x(i) = x(i) + (cube-step)*s(i)
317 else if (f_1 .eq. f_b) then
320 x(i) = x(i) + cube*s(i)
322 else if (f_1 .eq. f_c) then
326 c try to restart from best point with smaller stepsize
328 if (f_1 .gt. f_0) then
338 if (sg_1 .gt. 0.0d0) then
344 step = max(cube,step-cube) / 10.0d0
345 if (step .lt. stpmin) step = stpmin
347 c if already restarted once, then return with best point
349 if (status .eq. ' ReSearch') then