This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / e / rcspln64.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:27  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10 #if !defined(CERNLIB_DOUBLE)
11       SUBROUTINE RCSPLN(N,X,M,Y,NDIM,MODE,A,B,C,D)
12 #endif
13 #if defined(CERNLIB_DOUBLE)
14       SUBROUTINE DCSPLN(N,X,M,Y,NDIM,MODE,A,B,C,D)
15 #include "gen/imp64.inc"
16 #endif
17       CHARACTER NAMEN*(*),NAMET*(*)
18       CHARACTER*80 ERRTXT
19 #if !defined(CERNLIB_DOUBLE)
20       PARAMETER (NAMEN = 'RCSPLN', NAMET = 'RCSPLT')
21 #endif
22 #if defined(CERNLIB_DOUBLE)
23       PARAMETER (NAMEN = 'DCSPLN', NAMET = 'DCSPLT')
24 #endif
25       LOGICAL LNT
26
27       DIMENSION X(0:*),Y(0:NDIM,*)
28       DIMENSION A(NDIM,*),B(NDIM,*),C(NDIM,*),D(NDIM,*)
29
30       PARAMETER (Z1 = 1, C1 = 3*Z1/2, C2 = Z1/3, C3 = 2*Z1/3)
31       PARAMETER (C4 = Z1/2, C5 = Z1/6, C6 = 2*Z1/15, C7 = 7*Z1/60)
32
33       LNT=.FALSE.
34       GO TO 50
35
36 #if !defined(CERNLIB_DOUBLE)
37       ENTRY RCSPNT(N,X,M,Y,NDIM,MODE,A,B,C,D)
38 #endif
39 #if defined(CERNLIB_DOUBLE)
40       ENTRY DCSPNT(N,X,M,Y,NDIM,MODE,A,B,C,D)
41 #endif
42
43       LNT=.TRUE.
44    50 IF(N .LT. 2 .OR. M .LT. 1 .OR. NDIM .LT. N .OR.
45      1   MODE .NE. 0. AND. MODE .NE. 1) THEN
46        IF(N .LT. 2) THEN
47         WRITE(ERRTXT,101) N
48         IF(.NOT.LNT) CALL MTLPRT(NAMEN,'E211.1',ERRTXT)
49         IF(     LNT) CALL MTLPRT(NAMET,'E211.1',ERRTXT)
50        ENDIF
51        IF(M .LT. 1) THEN
52         WRITE(ERRTXT,102) M
53         IF(.NOT.LNT) CALL MTLPRT(NAMEN,'E211.2',ERRTXT)
54         IF(     LNT) CALL MTLPRT(NAMET,'E211.2',ERRTXT)
55        ENDIF
56        IF(NDIM .LT. N) THEN
57         WRITE(ERRTXT,103) NDIM,N
58         IF(.NOT.LNT) CALL MTLPRT(NAMEN,'E211.3',ERRTXT)
59         IF(     LNT) CALL MTLPRT(NAMET,'E211.3',ERRTXT)
60        ENDIF
61        IF(MODE .NE. 0 .AND. MODE .NE. 1) THEN
62         WRITE(ERRTXT,104) MODE
63         IF(.NOT.LNT) CALL MTLPRT(NAMEN,'E211.4',ERRTXT)
64         IF(     LNT) CALL MTLPRT(NAMET,'E211.4',ERRTXT)
65        ENDIF
66        RETURN
67       ENDIF
68       DO 1 I = 1,N
69     1 D(I,1)=X(I)-X(I-1)
70       DO 2 K = 1,M
71       DO 2 I = 1,N
72     2 A(I,K)=(Y(I,K)-Y(I-1,K))/D(I,1)
73       IF(MODE .EQ. 1) THEN
74        IF(N .EQ. 2) THEN
75         T1=1/(D(1,1)+D(2,1))
76         DO 3 K = 1,M
77     3   C(2,K)=T1*(A(2,K)-A(1,K))
78        ELSE
79         DO 4 I = 2,N
80     4   B(I,1)=1/(D(I,1)+D(I-1,1))
81         DO 5 K = 1,M
82     5   C(1,K)=0
83         B(1,1)=1
84
85         DO 6 I = 2,N-1
86         T1=3*B(I,1)
87         T2=B(I,1)*D(I-1,1)
88         T3=1/(2+T2*B(I-1,1))
89         B(I,1)=(T2-1)*T3
90         DO 6 K = 1,M
91     6   C(I,K)=(T1*(A(I,K)-A(I-1,K))-T2*C(I-1,K))*T3
92
93         T1=3*B(N,1)
94         T2=B(N,1)*D(N-1,1)
95         T3=1/(3-T2*(1-B(N-1,1)))
96         DO 7 K = 1,M
97     7   C(N,K)=(T1*(A(N,K)-A(N-1,K))-T2*C(N-1,K))*T3
98        END IF
99
100        DO 8 I = N-1,2,-1
101        T1=B(I,1)
102        DO 8 K = 1,M
103     8  C(I,K)=T1*C(I+1,K)+C(I,K)
104        DO 9 K = 1,M
105     9  C(1,K)=C(2,K)
106        IF(.NOT.LNT) THEN
107         DO 10 K = M,1,-1
108         B(1,K)=A(1,K)-C(2,K)*D(1,1)
109         D(1,K)=0
110         DO 11 I = 2,N-1
111         B(I,K)=A(I,K)-C2*(C(I+1,K)+2*C(I,K))*D(I,1)
112    11   D(I,K)=(C(I+1,K)-C(I,K))/(3*D(I,1))
113         B(N,K)=A(N,K)-C(N,K)*D(N,1)
114    10   D(N,K)=0
115         DO 12 K = 1,M
116         DO 12 I = 1,N
117    12   A(I,K)=Y(I-1,K)
118        ENDIF
119       ELSE
120        IF(N .EQ. 2) THEN
121         T1=C1/(D(1,1)+D(2,1))
122         DO 23 K = 1,M
123    23   C(2,K)=T1*(A(2,K)-A(1,K))
124        ELSE
125         DO 24 I = 2,N
126    24   B(I,1)=1/(D(I,1)+D(I-1,1))
127         DO 25 K = 1,M
128    25   C(1,K)=0
129         B(1,1)=0
130
131         DO 26 I = 2,N-1
132         T1=3*B(I,1)
133         T2=B(I,1)*D(I-1,1)
134         T3=1/(2+T2*B(I-1,1))
135         B(I,1)=(T2-1)*T3
136         DO 26 K = 1,M
137    26   C(I,K)=(T1*(A(I,K)-A(I-1,K))-T2*C(I-1,K))*T3
138
139         T1=3*B(N,1)
140         T2=B(N,1)*D(N-1,1)
141         T3=1/(2+T2*B(N-1,1))
142         DO 27 K = 1,M
143    27   C(N,K)=(T1*(A(N,K)-A(N-1,K))-T2*C(N-1,K))*T3
144        END IF
145
146        DO 28 I = N-1,2,-1
147        T1=B(I,1)
148        DO 28 K = 1,M
149    28  C(I,K)=T1*C(I+1,K)+C(I,K)
150        DO 29 K = 1,M
151    29  C(1,K)=0
152        IF(.NOT.LNT) THEN
153         T1=C2*D(1,1)
154         T2=C2/D(1,1)
155         T3=C3*D(N,1)
156         T4=C2/D(N,1)
157         DO 30 K = M,1,-1
158         B(1,K)=A(1,K)-T1*C(2,K)
159         D(1,K)=T2*C(2,K)
160         DO 31 I = 2,N-1
161         B(I,K)=A(I,K)-C2*(C(I+1,K)+2*C(I,K))*D(I,1)
162    31   D(I,K)=(C(I+1,K)-C(I,K))/(3*D(I,1))
163         B(N,K)=A(N,K)-T3*C(N,K)
164    30   D(N,K)=-T4*C(N,K)
165         DO 32 K = 1,M
166         DO 32 I = 1,N
167    32   A(I,K)=Y(I-1,K)
168        ENDIF
169       ENDIF
170
171       IF(LNT) THEN
172        DO 41 K = 1,M
173        T1=D(1,1)**2
174        A(1,K)=C4*(Y(0,K)+Y(1,K)-C5*(C(1,K)+C(2,K))*T1)*D(1,1)
175        B(1,K)=C2*(Y(0,K)+C4*Y(1,K)-(C6*C(1,K)+C7*C(2,K))*T1)*T1
176        DO 42 I = 2,N-1
177        T1=D(I,1)**2
178        A(I,K)=A(I-1,K)+
179      1        C4*(Y(I-1,K)+Y(I,K)-C5*(C(I,K)+C(I+1,K))*T1)*D(I,1)
180    42  B(I,K)=B(I-1,K)+C2*(Y(I-1,K)+C4*Y(I,K)-
181      1                (C6*C(I,K)+C7*C(I+1,K))*T1)*T1+A(I-1,K)*D(I,1)
182        T1=D(N,1)**2
183        A(N,K)=A(N-1,K)+
184      1        C4*(Y(N-1,K)+Y(N,K)-C(N,K)*C5*(1+MODE)*T1)*D(N,1)
185    41  B(N,K)=B(N-1,K)+C2*(Y(N-1,K)+C4*Y(N,K)-
186      1        C(N,K)*(C6+MODE*C7)*T1)*T1+A(N-1,K)*D(N,1)
187       ENDIF
188       RETURN
189   101 FORMAT('ILLEGAL N = ',I5,' < 2')
190   102 FORMAT('ILLEGAL M = ',I5,' < 1')
191   103 FORMAT('ILLEGAL NDIM = ',I5,' < ',I5,' = N')
192   104 FORMAT('ILLEGAL MODE = ',I5)
193       END