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 = 0d0
76 DO 1 i = 1, 5
77 CALL rand48(ix,iy,y)
78 ix = iy
79 sum = sum + y
80 1 CONTINUE
81 x = 6d0 * (sum / 5d0) - 3d0
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)
|
