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