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