V3FIT
fmax.f
1  FUNCTION fmax (ax, bx, f, tol)
2  USE stel_kinds
3  IMPLICIT NONE
4 C-----------------------------------------------
5 C D u m m y A r g u m e n t s
6 C-----------------------------------------------
7  REAL(rprec) ax, bx, f, tol
8 C-----------------------------------------------
9 C L o c a l P a r a m e t e r s
10 C-----------------------------------------------
11  REAL(rprec), PARAMETER :: zero=0, one=1, two=2,
12  1 three = 3, five = 5, p5 = 0.5_dp
13 C-----------------------------------------------
14 C L o c a l V a r i a b l e s
15 C-----------------------------------------------
16  REAL(rprec):: a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,
17  1 u,v,w,fu,fv,fw,fx,x, fmax
18 C-----------------------------------------------
19 C E x t e r n a l F u n c t i o n s
20 C-----------------------------------------------
21  EXTERNAL f
22 C-----------------------------------------------
23 c
24 c c is the squared inverse of the golden ratio
25 c
26  c = p5*(three - sqrt(five))
27 c
28 c eps is approximately the square root of the relative machine
29 c precision.
30 c
31  eps = one
32  eps = eps/two
33  tol1 = one + eps
34  1003 CONTINUE
35  IF (tol1 .le. one) GOTO 1002
36  eps = eps/two
37  tol1 = one + eps
38  IF (tol1 .le. one) GOTO 1002
39  eps = eps/two
40  tol1 = one + eps
41  IF (tol1 .le. one) GOTO 1002
42  eps = eps/two
43  tol1 = one + eps
44  IF (tol1 .le. one) GOTO 1002
45  eps = eps/two
46  tol1 = one + eps
47  GOTO 1003
48  1002 CONTINUE
49  eps = sqrt(eps)
50 c
51 c initialization
52 c
53  a = ax
54  b = bx
55  v = a + c*(b - a)
56  w = v
57  x = v
58  e = zero
59  fx = f(x)
60  fv = fx
61  fw = fx
62 c
63 c main loop starts here
64 c
65  20 CONTINUE
66  xm = p5*(a + b)
67  tol1 = eps*abs(x) + tol/three
68  tol2 = two*tol1
69 c
70 c check STOPping criterion
71 c
72  IF (abs(x - xm) .le. tol2 - p5*(b - a)) GOTO 90
73 c
74 c is golden-section necessary
75 c
76  IF (abs(e) .le. tol1) GOTO 40
77 c
78 c fit parabola
79 c
80  r = (x - w)*(fx - fv)
81  q = (x - v)*(fx - fw)
82  p = (x - v)*q - (x - w)*r
83  q = two*(q - r)
84 cwie IF (q.gt.0.0) p = -p
85  IF (q < zero) p = -p
86  q = abs(q)
87  r = e
88  e = d
89 c
90 c is parabola acceptable
91 c
92 cwie 30 IF (ABS(p).ge.ABS(0.5*q*r)) GOTO 40
93 cwie IF (p.le.q*(a - x)) GOTO 40
94 cwie IF (p.ge.q*(b - x)) GOTO 40
95  IF (abs(p) .le. abs(p5*q*r)) GOTO 40
96  IF (p .ge. q*(a - x)) GOTO 40
97  IF (p .le. q*(b - x)) GOTO 40
98 c
99 c a parabolic interpolation step
100 c
101  d = p/q
102  u = x + d
103 c
104 c f must not be evaluated too CLOSE to ax or bx
105 c
106  IF (u - a < tol2) d = sign(tol1,xm - x)
107  IF (b - u < tol2) d = sign(tol1,xm - x)
108  GOTO 50
109 c
110 c a golden-section step
111 c
112  40 CONTINUE
113  IF (x .ge. xm) THEN
114  e = a - x
115  ELSE
116  e = b - x
117  END IF
118  d = c*e
119 c
120 c f must not be evaluated too CLOSE to x
121 c
122  50 CONTINUE
123  IF (abs(d) .ge. tol1) THEN
124  u = x + d
125  ELSE
126  u = x + sign(tol1,d)
127  END IF
128  fu = f(u)
129 c
130 c update a, b, v, w, and x
131 c
132 cwie IF (fu.gt.fx) GOTO 60
133  IF (fu < fx) GOTO 60
134  IF (u .ge. x) THEN
135  a = x
136  ELSE
137  b = x
138  END IF
139  v = w
140  fv = fw
141  w = x
142  fw = fx
143  x = u
144  fx = fu
145  GOTO 20
146  60 CONTINUE
147  IF (u < x) THEN
148  a = u
149  ELSE
150  b = u
151  END IF
152 cwie IF (fu.le.fw) GOTO 70
153  IF (fu .ge. fw) GOTO 70
154  IF (w .eq. x) GOTO 70
155 cwie IF (fu.le.fv) GOTO 80
156  IF (fu .ge. fv) GOTO 80
157  IF (v .eq. x) GOTO 80
158  IF (v .eq. w) GOTO 80
159  GOTO 20
160  70 CONTINUE
161  v = w
162  fv = fw
163  w = u
164  fw = fu
165  GOTO 20
166  80 CONTINUE
167  v = u
168  fv = fu
169  GOTO 20
170 c
171 c END of main loop
172 c
173  90 CONTINUE
174  fmax = x
175 
176  END FUNCTION fmax