Blame view
fvn_fnlib/comp3.f
1.61 KB
38581db0c git-svn-id: https... |
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 50 51 52 53 54 55 |
subroutine comp3 (fn, dfn, xstart, xend, xincr, ystart, yend, 1 yincr, zstart, zend, zincr, fname) c august 1979 edition. w. fullerton, c3, los alamos scientific lab. double precision dfn, xstart, xend, xincr, ystart, yend, yincr, 1 zstart, zend, zincr, dx, dy, dz, df integer fname(1), fmt(23) external fn, dfn, i1mach c call s9comp (48, 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, 12x, 1hz, 69x, 4haerr, 7x, 4hrerr) c aerrmx = 0.0 rerrmx = 0.0 do 50 i=1,11 dx = xstart + dble(float(i-1))*xincr if (dx.gt.xend) go to 60 c do 40 j=1,11 dy = ystart + dble(float(j-1))*yincr if (dy.gt.yend) go to 50 c write (iwunit, 20) do 30 k=1,11 dz = zstart + dble(float(k-1))*zincr if (dz.gt.zend) go to 40 c x = dx y = dy z = dz f = fn (x, y, z) df = dfn (dble(x), dble(y), dble(z)) 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, dz) c write (iwunit, 20) x, y, z, f, df, aerr, rerr 20 format (1x, 3e13.5, e22.14, d36.28, 2x, 2e11.3) 30 continue 40 continue 50 continue c 60 write (iwunit, 70) aerrmx, rerrmx 70 format (/100x, 2e11.3) c return end |