]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 |