MINICERN/mathlib/gen/d/deriv64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:19  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if !defined(CERNLIB_DOUBLE)
  11       SUBROUTINE RDERIV(F,X,DELTA,DFDX,RERR)
  12 #endif
  13 #if defined(CERNLIB_DOUBLE)
  14       SUBROUTINE DDERIV(F,X,DELTA,DFDX,RERR)
  15 #include "gen/imp64.inc"
  16 #endif
  17 C     Computes the derivative f'(x) of f(x) at x = X. Based on
  18 C     H. Rutishauser, Ausdehnung des Rombergschen Prinzips
  19 C     (Extension of Romberg's Principle), Numer. Math. 5 (1963) 48-54
  20
  21
  22       DIMENSION DX(0:9),W(0:9,3),T(0:9,0:9),A(0:9)
  23       LOGICAL LEV(0:9),LMT
  24
  25       CHARACTER NAME*(*)
  26       CHARACTER*80 ERRTXT
  27 #if !defined(CERNLIB_DOUBLE)
  28       PARAMETER (NAME = 'RDERIV')
  29       PARAMETER (EPS = 5D-12)
  30 #endif
  31 #if defined(CERNLIB_DOUBLE)
  32       PARAMETER (NAME = 'DDERIV')
  33       PARAMETER (EPS = 5D-14)
  34 #endif
  35       PARAMETER (Z1 = 1, S = Z1/10)
  36
  37       DATA DX /0.0256D0, 0.0192D0, 0.0128D0, 0.0096D0, 0.0064D0,
  38      1         0.0048D0, 0.0032D0, 0.0024D0, 0.0016D0, 0.0012D0/
  39
  40       DATA (LEV(K),K=0,8,2) /5*.TRUE./
  41       DATA (LEV(K),K=1,9,2) /5*.FALSE./
  42
  43       DATA W(1,1) /1.33333 33333 333333D+00/
  44       DATA W(3,1) /1.06666 66666 666667D+00/
  45       DATA W(5,1) /1.01587 30158 730159D+00/
  46       DATA W(7,1) /1.00392 15686 274510D+00/
  47
  48       DATA W(2,1) /3.33333 33333 333333D-01/
  49       DATA W(4,1) /6.66666 66666 666667D-02/
  50       DATA W(6,1) /1.58730 15873 015873D-02/
  51       DATA W(8,1) /3.92156 86274 509804D-03/
  52
  53       DATA W(0,2) /2.28571 42857 142857D+00/
  54       DATA W(2,2) /1.16363 63636 363636D+00/
  55       DATA W(4,2) /1.03643 72469 635628D+00/
  56       DATA W(6,2) /1.00886 69950 738916D+00/
  57       DATA W(8,2) /1.00220 21042 329337D+00/
  58
  59       DATA W(1,2) /1.28571 42857 142857D+00/
  60       DATA W(3,2) /1.63636 36363 636364D-01/
  61       DATA W(5,2) /3.64372 46963 562753D-02/
  62       DATA W(7,2) /8.86699 50738 916256D-03/
  63       DATA W(9,2) /2.20210 42329 336922D-03/
  64
  65       DATA W(0,3) /1.80000 00000 000000D+00/
  66       DATA W(2,3) /1.12500 00000 000000D+00/
  67       DATA W(4,3) /1.02857 14285 714286D+00/
  68       DATA W(6,3) /1.00699 30069 930070D+00/
  69       DATA W(8,3) /1.00173 91304 347826D+00/
  70
  71       DATA W(1,3) /8.00000 00000 000000D-01/
  72       DATA W(3,3) /1.25000 00000 000000D-01/
  73       DATA W(5,3) /2.85714 28571 428571D-02/
  74       DATA W(7,3) /6.99300 69930 069930D-03/
  75       DATA W(9,3) /1.73913 04347 826087D-03/
  76
  77 #if !defined(CERNLIB_DOUBLE)
  78       ENTRY DERIV(F,X,DELTA,DFDX,RERR)
  79 #endif
  80
  81       DEL=10*ABS(DELTA)
  82       IS=10
  83
  84     4 IS=IS-1
  85       DEL=S*DEL
  86       IF(IS .EQ. 0 .OR. X+DEL*DX(9) .EQ. X) THEN
  87        DELTA=DEL
  88        DFDX=0
  89        RERR=1
  90        WRITE(ERRTXT,101) X
  91        CALL MTLPRT(NAME,'D401.1',ERRTXT)
  92        RETURN
  93       ENDIF
  94       DO 1 K = 0,9
  95       H=DEL*DX(K)
  96       T(K,0)=(F(X+H)-F(X-H))/(H+H)
  97     1 A(K)=T(K,0)
  98
  99       IF(A(0) .GE. A(9)) THEN
 100        DO 5 K = 0,9
 101     5  A(K)=-A(K)
 102       ENDIF
 103
 104       LMT=.TRUE.
 105       DO 3 K = 1,9
 106       H=A(K-1)-A(K)
 107     3 LMT=LMT .AND. (H .LE. 0 .OR. ABS(H) .LE. EPS*ABS(A(K)))
 108       IF(.NOT.LMT) GO TO 4
 109
 110       DO 2 M = 1,9
 111       DO 2 K = 0,9-M
 112       IF(LEV(M)) THEN
 113        T(K,M)=W(M-1,1)*T(K+1,M-1)-W(M,1)*T(K,M-1)
 114       ELSEIF(LEV(K)) THEN
 115        T(K,M)=W(M-1,2)*T(K+1,M-1)-W(M,2)*T(K,M-1)
 116       ELSE
 117        T(K,M)=W(M-1,3)*T(K+1,M-1)-W(M,3)*T(K,M-1)
 118       ENDIF
 119     2 CONTINUE
 120       DFDX=T(0,9)
 121       RERR=0
 122       IF(DFDX .NE. 0) RERR=(DFDX-T(0,8))/DFDX
 123       DELTA=DEL
 124       RETURN
 125   101 FORMAT('FAILURE FOR X = ',1P,D15.8)
 126       END