comp2.f
1.36 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
subroutine comp2 (fn, dfn, xstart, xend, xincr, ystart, yend,
1 yincr, fname)
c august 1979 edition. w. fullerton, c3, los alamos scientific lab.
double precision dfn, xstart, xend, xincr, ystart, yend, yincr,
1 dx, dy, df
integer fname(1), fmt(23)
external fn, dfn, i1mach
c
call s9comp (39, fmt, nwd)
iwunit = i1mach(2)
write (iwunit, 5)
5 format (1h1)
write (iwunit, fmt) (fname(i), i=1,nwd), (fname(i), i=1,nwd)
write (iwunit, 10)
10 format (1h+, 8x, 1hx, 12x, 1hy, 78x, 4haerr, 8x, 4hrerr)
c
aerrmx = 0.0
rerrmx = 0.0
do 40 i=1,21
dx = xstart + dble(float(i-1))*xincr
if (dx.gt.xend) go to 50
c
write (iwunit, 20)
do 30 j=1,21
dy = ystart + dble(float(j-1))*yincr
if (dy.gt.yend) go to 40
c
x = dx
y = dy
f = fn (x, y)
df = dfn (dble(x), dble(y))
aerr = df - dble(f)
rerr = 0.0
if (df.ne.0.0d0) rerr = abs (aerr/sngl(df))
c
aerrmx = amax1 (abs(aerr), aerrmx)
rerrmx = amax1 (rerr, rerrmx)
df = dfn (dx, dy)
c
write (iwunit, 20) x, y, f, df, aerr, rerr
20 format (1x, 2e13.5, e25.15, d40.30, 3x, 2e12.3)
30 continue
40 continue
c
50 write (iwunit, 60) aerrmx, rerrmx
60 format (/95x, 2e12.3)
c
return
end