Blame view

fvn_fnlib/comp3.f 1.61 KB
38581db0c   daniau   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