vif-math.f
1 ! @section Synopsis
2 !
3 ! Auxiliary math functions for VIF.
4 !
5 ! @author J. Marcel van der Veer
6 !
7 ! @section copyright
8 !
9 ! This file is part of VIF - vintage fortran compiler.
10 ! Copyright 2020-2025 J. Marcel van der Veer <algol68g@xs4all.nl>.
11 !
12 ! @section license
13 !
14 ! This program is free software; you can redistribute it and/or modify it
15 ! under the terms of the gnu general public license as published by the
16 ! free software foundation; either version 3 of the license, or
17 ! (at your option) any later version.
18 !
19 ! This program is distributed in the hope that it will be useful, but
20 ! without any warranty; without even the implied warranty of merchantability
21 ! or fitness for a particular purpose. See the GNU general public license for
22 ! more details. You should have received a copy of the GNU general public
23 ! license along with this program. If not, see <http://www.gnu.org/licenses/>.
24 !
25
26 subroutine rand48(ix,iy,yfl)
27 !
28 ! Generates random numbers between 0 and 1 (SUN version).
29 !
30 ! ix - initial seed
31 ! iy - new seed value
32 ! yfl - random number double precision
33 !
34 ! Adapted from CASCADE: netlib.org/ieeecs/cascade
35 ! The CASCADE library is in the public domain.
36 !
37 double precision yfl, drand48
38 integer ix, iy, iseed
39 logical first
40 data first /.true./
41 save first, iseed
42 !
43 ! if this is the first call, initialize topmost 32 bits of the seed
44 !
45 if (first .or. iseed .ne. ix) then
46 call srand48 (ix)
47 first = .false.
48 end if
49 !
50 ! Use SUN-supplied function drand48
51 !
52 iy = ix
53 iseed = ix
54 yfl = drand48()
55 !
56 return
57 end
58
59 subroutine gauss(ix,x)
60 ! The routine gauss generates a random number from a gaussian
61 ! distribution, with a mean of zero and a standard deviation of one.
62 !
63 ! VAX/VMS version
64 !
65 ! Variables:
66 !
67 ! ix - initial seed integer; gauss updates it with a new seed value
68 !
69 ! x - the output gaussian variable
70 !
71 ! Adapted from CASCADE: netlib.org/ieeecs/cascade
72 ! The CASCADE library is in the public domain.
73 !
74 double precision x, y, sum
75 sum=0.0
76 do 10 i=1,5
77 call rand48(ix,iy,y)
78 ix = iy
79 sum=sum+y
80 10 continue
81 x=6.0d0*(sum/5.0d0)-3.0d0
82 return
83 end
84
85 real*8 function zero in (ax, bx, f, tol, ierr)
86 implicit real*8 (a - h, o - z)
87 !
88 ! Compute a zero of the function f(x) in the interval [ax, bx].
89 !
90 ! It is assumed that f(ax) and f(bx) have opposite signs.
91 ! Zeroin returns a zero in the given interval ax, bx to within a
92 ! tolerance 4 * macheps * abs(x) + tol, where macheps
93 ! is the relative machine precision.
94 !
95 ! Adapted from zeroin.f from netlib.org/fmm.
96 ! That function is a slightly modified translation of the ALGOL 60
97 ! procedure zero given in
98 !
99 ! Richard Brent, Algorithms for Minimization without Derivatives,
100 ! Prentice - Hall, inc. (1973).
101 !
102 ! Parameters:
103 !
104 ! ax left endpoint of initial interval
105 ! bx right endpoint of initial interval
106 ! f function subprogram which evaluates f(x) for any x in the interval ax, bx
107 ! tol desired length of the interval of uncertainty of the final result (>= 0.0d0)
108 !
109 ! Result:
110 !
111 ! Approximation of a zero of f in the interval ax, bx
112 ! ierr = 0 'zero in' found a root
113 ! ierr = 1 'zero in' did not find a root
114 !
115 ! Compare to SLATEC routines DFZERO and FZERO.
116 !
117
118 ! Initialization.
119 eps = d1mach(3)
120 a = ax
121 b = bx
122 fa = f(a)
123 fb = f(b)
124 ! Begin.
125 1 c = a
126 fc = fa
127 d = b - a
128 e = d
129 2 if (dabs(fc) < dabs(fb)) then
130 a = b
131 b = c
132 c = a
133 fa = fb
134 fb = fc
135 fc = fa
136 end if
137 ! Convergence check.
138 tol1 = 2.0d0 * eps * dabs(b) + 0.5d0 * tol
139 xm = .5d0 * (c - b)
140 if (dabs(xm) <= tol1 | fb == 0.0d0) go to 5
141 ! Is bisection necessary?
142 if (dabs(e) < tol1 | dabs(fa) <= dabs(fb)) go to 3
143 ! Is quadratic interpolation possible?
144 if (a == c) then
145 ! Linear interpolation.
146 s = fb / fa
147 p = 2.0d0 * xm * s
148 q = 1.0d0 - s
149 else
150 ! Inverse quadratic interpolation.
151 q = fa / fc
152 r = fb / fc
153 s = fb / fa
154 p = s * (2.0d0 * xm * q * (q - r) - (b - a) * (r - 1.0d0))
155 q = (q - 1.0d0) * (r - 1.0d0) * (s - 1.0d0)
156 end if
157 ! Adjust signs.
158 if (p > 0.0d0) q = -q
159 p = dabs(p)
160 ! Is interpolation acceptable?
161 if ((2.0d0 * p) >= (3.0d0 * xm * q - dabs(tol1 * q))) go to 3
162 if (p >= dabs(0.5d0 * e * q)) go to 3
163 e = d
164 d = p / q
165 go to 4
166 ! Bisection.
167 3 d = xm
168 e = d
169 ! Complete step.
170 4 a = b
171 fa = fb
172 if (dabs(d) > tol1) b = b + d
173 if (dabs(d) <= tol1) b = b + dsign(tol1, xm)
174 fb = f(b)
175 if ((fb * (fc / dabs(fc))) > 0.0d0) go to 1
176 go to 2
177 ! Done.
178 5 zero in = b
179 if (dabs (fb) <= tol) then
180 ierr = 0
181 else
182 ierr = 1
183 end if
184 return
185 end
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
