*
* $Id: rcspln64.F,v 1.1.1.1 1996/04/01 15:02:27 mclareni Exp $
*
* $Log: rcspln64.F,v $
* Revision 1.1.1.1  1996/04/01 15:02:27  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
#if !defined(CERNLIB_DOUBLE)
      SUBROUTINE RCSPLN(N,X,M,Y,NDIM,MODE,A,B,C,D)
#endif
#if defined(CERNLIB_DOUBLE)
      SUBROUTINE DCSPLN(N,X,M,Y,NDIM,MODE,A,B,C,D)
#include "gen/imp64.inc"
#endif
      CHARACTER NAMEN*(*),NAMET*(*)
      CHARACTER*80 ERRTXT
#if !defined(CERNLIB_DOUBLE)
      PARAMETER (NAMEN = 'RCSPLN', NAMET = 'RCSPLT')
#endif
#if defined(CERNLIB_DOUBLE)
      PARAMETER (NAMEN = 'DCSPLN', NAMET = 'DCSPLT')
#endif
      LOGICAL LNT

      DIMENSION X(0:*),Y(0:NDIM,*)
      DIMENSION A(NDIM,*),B(NDIM,*),C(NDIM,*),D(NDIM,*)

      PARAMETER (Z1 = 1, C1 = 3*Z1/2, C2 = Z1/3, C3 = 2*Z1/3)
      PARAMETER (C4 = Z1/2, C5 = Z1/6, C6 = 2*Z1/15, C7 = 7*Z1/60)

      LNT=.FALSE.
      GO TO 50

#if !defined(CERNLIB_DOUBLE)
      ENTRY RCSPNT(N,X,M,Y,NDIM,MODE,A,B,C,D)
#endif
#if defined(CERNLIB_DOUBLE)
      ENTRY DCSPNT(N,X,M,Y,NDIM,MODE,A,B,C,D)
#endif

      LNT=.TRUE.
   50 IF(N .LT. 2 .OR. M .LT. 1 .OR. NDIM .LT. N .OR.
     1   MODE .NE. 0. AND. MODE .NE. 1) THEN
       IF(N .LT. 2) THEN
        WRITE(ERRTXT,101) N
        IF(.NOT.LNT) CALL MTLPRT(NAMEN,'E211.1',ERRTXT)
        IF(     LNT) CALL MTLPRT(NAMET,'E211.1',ERRTXT)
       ENDIF
       IF(M .LT. 1) THEN
        WRITE(ERRTXT,102) M
        IF(.NOT.LNT) CALL MTLPRT(NAMEN,'E211.2',ERRTXT)
        IF(     LNT) CALL MTLPRT(NAMET,'E211.2',ERRTXT)
       ENDIF
       IF(NDIM .LT. N) THEN
        WRITE(ERRTXT,103) NDIM,N
        IF(.NOT.LNT) CALL MTLPRT(NAMEN,'E211.3',ERRTXT)
        IF(     LNT) CALL MTLPRT(NAMET,'E211.3',ERRTXT)
       ENDIF
       IF(MODE .NE. 0 .AND. MODE .NE. 1) THEN
        WRITE(ERRTXT,104) MODE
        IF(.NOT.LNT) CALL MTLPRT(NAMEN,'E211.4',ERRTXT)
        IF(     LNT) CALL MTLPRT(NAMET,'E211.4',ERRTXT)
       ENDIF
       RETURN
      ENDIF
      DO 1 I = 1,N
    1 D(I,1)=X(I)-X(I-1)
      DO 2 K = 1,M
      DO 2 I = 1,N
    2 A(I,K)=(Y(I,K)-Y(I-1,K))/D(I,1)
      IF(MODE .EQ. 1) THEN
       IF(N .EQ. 2) THEN
        T1=1/(D(1,1)+D(2,1))
        DO 3 K = 1,M
    3   C(2,K)=T1*(A(2,K)-A(1,K))
       ELSE
        DO 4 I = 2,N
    4   B(I,1)=1/(D(I,1)+D(I-1,1))
        DO 5 K = 1,M
    5   C(1,K)=0
        B(1,1)=1

        DO 6 I = 2,N-1
        T1=3*B(I,1)
        T2=B(I,1)*D(I-1,1)
        T3=1/(2+T2*B(I-1,1))
        B(I,1)=(T2-1)*T3
        DO 6 K = 1,M
    6   C(I,K)=(T1*(A(I,K)-A(I-1,K))-T2*C(I-1,K))*T3

        T1=3*B(N,1)
        T2=B(N,1)*D(N-1,1)
        T3=1/(3-T2*(1-B(N-1,1)))
        DO 7 K = 1,M
    7   C(N,K)=(T1*(A(N,K)-A(N-1,K))-T2*C(N-1,K))*T3
       END IF

       DO 8 I = N-1,2,-1
       T1=B(I,1)
       DO 8 K = 1,M
    8  C(I,K)=T1*C(I+1,K)+C(I,K)
       DO 9 K = 1,M
    9  C(1,K)=C(2,K)
       IF(.NOT.LNT) THEN
        DO 10 K = M,1,-1
        B(1,K)=A(1,K)-C(2,K)*D(1,1)
        D(1,K)=0
        DO 11 I = 2,N-1
        B(I,K)=A(I,K)-C2*(C(I+1,K)+2*C(I,K))*D(I,1)
   11   D(I,K)=(C(I+1,K)-C(I,K))/(3*D(I,1))
        B(N,K)=A(N,K)-C(N,K)*D(N,1)
   10   D(N,K)=0
        DO 12 K = 1,M
        DO 12 I = 1,N
   12   A(I,K)=Y(I-1,K)
       ENDIF
      ELSE
       IF(N .EQ. 2) THEN
        T1=C1/(D(1,1)+D(2,1))
        DO 23 K = 1,M
   23   C(2,K)=T1*(A(2,K)-A(1,K))
       ELSE
        DO 24 I = 2,N
   24   B(I,1)=1/(D(I,1)+D(I-1,1))
        DO 25 K = 1,M
   25   C(1,K)=0
        B(1,1)=0

        DO 26 I = 2,N-1
        T1=3*B(I,1)
        T2=B(I,1)*D(I-1,1)
        T3=1/(2+T2*B(I-1,1))
        B(I,1)=(T2-1)*T3
        DO 26 K = 1,M
   26   C(I,K)=(T1*(A(I,K)-A(I-1,K))-T2*C(I-1,K))*T3

        T1=3*B(N,1)
        T2=B(N,1)*D(N-1,1)
        T3=1/(2+T2*B(N-1,1))
        DO 27 K = 1,M
   27   C(N,K)=(T1*(A(N,K)-A(N-1,K))-T2*C(N-1,K))*T3
       END IF

       DO 28 I = N-1,2,-1
       T1=B(I,1)
       DO 28 K = 1,M
   28  C(I,K)=T1*C(I+1,K)+C(I,K)
       DO 29 K = 1,M
   29  C(1,K)=0
       IF(.NOT.LNT) THEN
        T1=C2*D(1,1)
        T2=C2/D(1,1)
        T3=C3*D(N,1)
        T4=C2/D(N,1)
        DO 30 K = M,1,-1
        B(1,K)=A(1,K)-T1*C(2,K)
        D(1,K)=T2*C(2,K)
        DO 31 I = 2,N-1
        B(I,K)=A(I,K)-C2*(C(I+1,K)+2*C(I,K))*D(I,1)
   31   D(I,K)=(C(I+1,K)-C(I,K))/(3*D(I,1))
        B(N,K)=A(N,K)-T3*C(N,K)
   30   D(N,K)=-T4*C(N,K)
        DO 32 K = 1,M
        DO 32 I = 1,N
   32   A(I,K)=Y(I-1,K)
       ENDIF
      ENDIF

      IF(LNT) THEN
       DO 41 K = 1,M
       T1=D(1,1)**2
       A(1,K)=C4*(Y(0,K)+Y(1,K)-C5*(C(1,K)+C(2,K))*T1)*D(1,1)
       B(1,K)=C2*(Y(0,K)+C4*Y(1,K)-(C6*C(1,K)+C7*C(2,K))*T1)*T1
       DO 42 I = 2,N-1
       T1=D(I,1)**2
       A(I,K)=A(I-1,K)+
     1        C4*(Y(I-1,K)+Y(I,K)-C5*(C(I,K)+C(I+1,K))*T1)*D(I,1)
   42  B(I,K)=B(I-1,K)+C2*(Y(I-1,K)+C4*Y(I,K)-
     1                (C6*C(I,K)+C7*C(I+1,K))*T1)*T1+A(I-1,K)*D(I,1)
       T1=D(N,1)**2
       A(N,K)=A(N-1,K)+
     1        C4*(Y(N-1,K)+Y(N,K)-C(N,K)*C5*(1+MODE)*T1)*D(N,1)
   41  B(N,K)=B(N-1,K)+C2*(Y(N-1,K)+C4*Y(N,K)-
     1        C(N,K)*(C6+MODE*C7)*T1)*T1+A(N-1,K)*D(N,1)
      ENDIF
      RETURN
  101 FORMAT('ILLEGAL N = ',I5,' < 2')
  102 FORMAT('ILLEGAL M = ',I5,' < 1')
  103 FORMAT('ILLEGAL NDIM = ',I5,' < ',I5,' = N')
  104 FORMAT('ILLEGAL MODE = ',I5)
      END