1      double precision function linesrch(n,x,f,g)
2*
3* $Id$
4*
5      implicit none
6      integer n
7      double precision x(n)
8      double precision f(n)
9      double precision g(n)
10      integer ia, ib, ic
11      double precision dx, tol, prec
12      double precision step
13      logical linebracket
14      external linebracket
15
16      write(6,*)
17      write(6,*)
18      do ia=max(1,(n-4)),n
19        write(6,899) ia,x(ia),f(ia),g(ia)
20 899    format(i5,2f24.12,e20.6)
21      enddo
22      if (n.eq.2) then
23        dx = (g(n) - g(n-1))/(x(n)-x(n-1))
24        write(6,884) dx
25 884    format('Curvature:',e15.4)
26        if ((dx.lt.0.d0).and.(f(n).gt.f(n-1))) then
27          linesrch = -x(n)
28          return
29        endif
30      endif
31      tol = 1.d-13
32      step = (g(n)*x(n-1) - g(n-1)*x(n))/(g(n) - g(n-1))
33c
34c  Bracketing
35c
36      if (n.gt.2) then
37        if (linebracket(n,x,f,ia,ib,ic)) then
38          write(6,900) ib,x(ib),ia,x(ia),(f(ia)-f(ib)),
39     $                          ic,x(ic),(f(ic)-f(ib))
40 900      format(  'Minimum:',10x,i5,f22.14
41     $           /,'Left Bracket:',5x,i5,f22.14,5x,e14.6,
42     $           /,'Right Bracket:',4x,i5,f22.14,5x,e14.6)
43          write(6,771) step,n,n-1
44 771      format(5x,'Step:',f20.10,' from ',i3,'-',i3)
45c
46c  If last step is not the minimum cannot use it
47c  in NR extrapolation, use current brackets and
48c  sign of gradient for extrapolation.
49c
50          if (ib.ne.n) then
51            if (g(ib).gt.0) then
52              step = (g(ia)*x(ib) - g(ib)*x(ia))/(g(ia) - g(ib))
53              write(6,771) step,ia,ib
54            else
55              step = (g(ic)*x(ib) - g(ib)*x(ic))/(g(ic) - g(ib))
56              write(6,771) step,ib,ic
57            endif
58          endif
59c
60c  Apply bracketing criterion only if brackets
61c  are greater than sqrt(relative machine-precision)*x
62c  and difference in function values is greater than machine precision
63c
64          prec = abs(sqrt(tol)*x(ib))
65          if (((step.lt.x(ia)).and.(abs(x(ib)-x(ia)).gt.prec).and.
66     $         (abs(f(ia)-f(ib)).gt.(tol*f(ib)))).or.
67     $         ((step.gt.x(ic)).and.(abs(x(ic)-x(ia)).gt.prec).and.
68     $         (abs(f(ic)-f(ib)).gt.(tol*f(ib))))) then
69            if (g(ib).gt.0) then
70              step = (x(ia)+x(ib))/2.d0
71            else
72              step = (x(ib)+x(ic))/2.d0
73            endif
74            write(6,772) step
75 772        format('Step outside bracket -- bisection step',f20.10)
76          endif
77        else
78c
79c  Cannot bracket, assume function is monotonic
80c  and use points with least function value
81c
82          write(6,774) ib,x(ib),f(ib)
83 774      format('Minimum: ',i5,2f20.10)
84          if ((ia.ne.0).and.(g(ib).lt.0.0d0)) then
85            dx = (g(ia) - g(ib))/(x(ia)-x(ib))
86            write(6,775) ia,x(ia),f(ia),dx
87 775        format('Extrapolate right:',i5,2f20.10,5x,
88     $             '2nd deriv',2x,e12.3)
89            step = x(ib) + (x(ib) - x(1))
90          else if ((ic.ne.0).and.(g(ib).gt.0.d0)) then
91            dx = (g(ic) - g(ib))/(x(ic)-x(ib))
92            write(6,776) ic,x(ic),f(ic),dx
93 776        format('Extrapolate left:',i5,2f20.10,5x,
94     $             '2nd deriv',2x,e12.3)
95            step = x(ib) - (x(1) - x(ib))
96          else
97            write(6,773)
98 773        format('Internal error....')
99          endif
100        endif
101      else if (n.eq.1) then
102        step = 1.d0
103      endif
104      linesrch = step
105      write(6,881)
106 881  format(//)
107      return
108      end
109
110
111
112
113
114
115      logical function linebracket(n,x,f,ia,ib,ic)
116      implicit none
117      integer n
118      double precision x(n),f(n)
119      integer ia,ib,ic
120      double precision xmin, ymin, del
121      integer i
122
123      ia = 0
124      ic = 0
125      ib = 1
126      xmin = f(1)
127      do i=2,n
128        if (f(i).lt.xmin) then
129          xmin = f(i)
130          ib = i
131        endif
132      enddo
133      xmin = -1.d0*abs(xmin)
134      ymin = abs(xmin)
135      do i=1,n
136        del = x(i) - x(ib)
137        if ((del.gt.0.d0).and.(del.lt.ymin)) then
138          ymin = del
139          ic = i
140        else if ((del.lt.0.d0).and.(del.gt.xmin)) then
141          xmin = del
142          ia = i
143        endif
144      enddo
145      linebracket = ((ia.ne.0).and.(ic.ne.0))
146      return
147      end
148
149
150
151
152
153
154
155
156
157
158      double precision function linesrch1(n,x,f,g)
159      implicit none
160#include "global.fh"
161      integer n
162      double precision x(n)
163      double precision f(n)
164      double precision g(n)
165      double precision a, b1, b2, c1, c2, f1, f2, x0
166      double precision e1, e2, err
167      save e1, e2
168      data err/0.d0/
169
170      if (n.eq.1) then
171        linesrch1 = 1.d0
172        return
173      else
174        if (n.gt.2) then
175          err = max(abs((f(n)-e1)/f(n)),abs((f(n)-e2)/f(n)))
176C          if (ga_nodeid().eq.0) write(6,901) x(n),e1,e2,f(n),err
177 901      format(//,'Predictions @ ',f12.6,5x,2f20.10,
178     $           /,'Obs.',f20.10,
179     $           /,'Relative error:',e12.2)
180        endif
181        a = (g(n)-g(n-1))/(2.d0*(x(n)-x(n-1)))
182        b1 = g(n-1) - 4.d0*a*x(n-1)
183        b2 = g(n) - 4.d0*a*x(n)
184        c1 = f(n-1) - a*x(n-1)**2 - b1*x(n-1)
185        c2 = f(n) - a*x(n)**2 - b2*x(n)
186        f1 = a*x(n-1)**2 + b1*x(n-1) + c1
187        f2 = a*x(n)**2 + b2*x(n) + c2
188C        if (ga_nodeid().eq.0)
189C    $     write(6,900) a,b1,b2,c1,c2,f1,f(n-1),f2,f(n)
190 900    format('A:   ',5x,f20.10,/,
191     $         'B(1):',5x,f20.10,/,
192     $         'B(2):',5x,f20.10,/,
193     $         'C(1):',5x,f20.10,/,
194     $         'C(2):',5x,f20.10,/,
195     $         'F(1):',5x,f20.10,'    f(1):',f20.10,/,
196     $         'F(2):',5x,f20.10,'    f(2):',f20.10)
197C        if (err.lt.1.d-1) then
198          x0 = (g(n)*x(n-1) - g(n-1)*x(n))/(g(n) - g(n-1))
199          e1 = a*x0**2 + b1*x0 + c1
200          e2 = a*x0**2 + b2*x0 + c2
201          linesrch1 = x0
202C        else
203
204      endif
205      return
206      end
207
208
209
210