[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / e / dspnb1.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:26  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      FUNCTION DSPNB1(K,M,I,NDER,X,T,NERR)

#include "gen/imp64.inc"
      DIMENSION T(*),B(27),D(27)
      CHARACTER NAME*(*)
      CHARACTER*80 ERRTXT
      PARAMETER (NAME = 'DSPNB1')

************************************************************************
*   NORBAS, VERSION: 15.03.1993
************************************************************************
*
*   DSPNB1 COMPUTES FUNCTION VALUES, VALUES OF DERIVATIVES, AND THE
*   VALUE OF INTEGRAL, RESPECTIVELY, OF NORMALIZED B-SPLINES
*                      B(I,K)(X)
*   OF DEGREE  K  ( 0<= K <= 25 )  WITH INDEX  I  ( 1 <= I <= M-K-1 )
*   OVER A SET OF SPLINE-KNOTS
*              T(1),T(2), ... ,T(M)    ( M >= 2*K+2 )
*   (KNOTS IN ASCENDING ORDER, WITH MULTIPLICITIES NOT GREATER
*   THAN  K+1).
*
*   THE FUNCTION VALUE OF THE NORMALIZED B-SPLINE  B(I,K)(X)  IS
*   IDENTICALLY ZERO OUTSIDE THE INTERVAL  T(I) <= X < T(I+K+1).
*
*   THE NORMALIZATION OF  N(X) = B(I,K)(X)  IS SUCH THAT THE INTGRAL OF
*   N(X)  OVER THE WHOLE X-RANGE EQUALS
*                  ( T(I+K+1) - T(I) ) / (K+1)  .
*
*   PARAMETERS:
*
*   K     (INTEGER) DEGREE (= ORDER - 1) OF B-SPLINES.
*   M     (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES.
*   NDER  (INTEGER) ON ENTRY, NDER MUST CONTAIN AN INTEGER VALUE .GE. -1
*         = -1: DSPNB1 COMPUTES THE INTEGRAL OF  B(I,K)(TAU)  OVER THE
*               RANGE  TAU <= X.
*         =  0: DSPNB1 COMPUTES THE FUNCTION VALUE  B(I,K)(X)  FOR
*               FOR THE SPECIFIED VALUES OF I, K, AND X.
*         >= 1: DSPNB1 COMPUTES THE VALUE OF THE NDER-TH DERIVATIVE OF
*               B(I,K)(X)  FOR THE SPECIFIED VALUES OF I, K, AND X
*               (IF  NDER > K  ZERO RETURNS).
*   X     (DOUBLE PRECISION) ON ENTRY, X MUST CONTAIN THE VALUE OF THE
*         INDEPENDENT VARIABLE X OF  B(I,K)(X)
*   T     (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER M CONTAINING THE
*         KNOTS, ON ENTRY.
*   I     (INTEGER) INDEX OF THE B-SPLINE  B(I,K)(X)
*   NERR  (INTEGER) ERROR INDICATOR. ON EXIT:
*         = 0: NO ERROR DETECTED
*         = 1: AT LEAST ONE OF THE CONSTANTS K , M , I , NDER IS ILLEGAL
*
*   USAGE:
*
*   THE FUNCTION-CALL
*         A = DSPNB1(K,M,I,NDER,X,T,NERR)
*   RETURNS
*       - THE VALUE OF THE INTEGRAL           (NDER = -1) OR
*       - THE FUNCTION VALUE                  (NDER = 0 ) OR
*       - THE VALUE OF THE NDER-TH DERIVATIVE (NDER > 0 )
*   OF THE NORMALIZED B-SPLINE  B(I,K)(X) OF DEGREE K WITH INDEX I AT X.
*
*   ERROR MESSAGES:
*
*   IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
*   PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
*     K < 0      OR    K > 25    OR
*     M < 2*K+2  OR
*     NDER < -1  OR
*     I < 1      OR    I > M-K-1
*
************************************************************************

      PARAMETER (Z0 = 0 , Z1 = 1)

      NERR=1
      IF(K .LT. 0 .OR. K .GT. 25) THEN
       WRITE(ERRTXT,101) 'K',K
       CALL MTLPRT(NAME,'E210.1',ERRTXT)
      ELSEIF(M .LT. 2*K+2) THEN
       WRITE(ERRTXT,101) 'M',M
       CALL MTLPRT(NAME,'E210.2',ERRTXT)
      ELSEIF(I . LT. 1 .OR. I .GT. M-K-1) THEN
       WRITE(ERRTXT,101) 'I',I
       CALL MTLPRT(NAME,'E210.3',ERRTXT)
      ELSEIF(NDER .LT. -1) THEN
       WRITE(ERRTXT,101) 'NDER',NDER
       CALL MTLPRT(NAME,'E210.5',ERRTXT)
      ELSE

       NERR=0
       DSPNB1=Z0
       IF(     X .LT. T(I)
     +    .OR. X .GT. T(I+K+1) .AND. NDER .GE. 0
     +    .OR. K .LT. NDER                       ) RETURN

       IF(NDER .EQ. -1) THEN
        IF(X .GE. T(I+K+1)) THEN
         R=(T(I+K+1)-T(I))/(K+1)
        ELSE
         KK=LKKSPL(X,T(I),MIN(2*(K+1),M-K-I))+I-1
         IF(K .EQ. 0) THEN
          R=X-T(KK-1)
         ELSE
          CALL DVSET(K+1,Z0,B(1),B(2))
          B(KK-I)=1/(T(KK)-T(KK-1))
          DO 10 L=1,K
          DO 10 J=MAX(1,KK-I-L),MIN(K+1-L,KK-I)
          DIF=T(I+J+L)-T(I+J-1)
          B0=Z0
          IF(DIF .NE. 0) B0=((X-T(I+J-1))*B(J)+(T(I+J+L)-X)*B(J+1))/DIF
   10     B(J)=B0
          S=Z0
          DO 20 L=1,KK-I
   20     S=S+(X-T(I+L-1))*B(L)
          R=S*(T(I+K+1)-T(I))/(K+1)
         ENDIF
        ENDIF
        DSPNB1=R
        RETURN
       ENDIF

       IF(K .EQ. 0) THEN
        R=Z1
       ELSE
        KK=LKKSPL(X,T(I),MIN(2*(K+1),M-K-I))+I-1
        I0=I+K+2-KK
        IF(I0 .EQ. 0) THEN
         R=Z0
        ELSE
         E1=X-T(KK-1)
         B(1)=Z1
         DO 30 J=2,K-NDER+1
   30    B(J)=E1*B(J-1)/(T(KK-2+J)-T(KK-1))
         IF(KK .NE. I+1 .OR. NDER .NE. 0) THEN
          E2=T(KK)-X
          DO 40 J=1,K-NDER
          E3=X-T(KK-1-J)
          B(1)=E2*B(1)/(T(KK)-T(KK-J))
          DO 40 L=2,K-NDER+1-J
   40     B(L)=E3*B(L-1)/(T(KK-2+L)-T(KK-1-J))+
     +         (T(KK-1+L)-X)*B(L)/(T(KK-1+L)-T(KK-J))
         ENDIF
         IF(NDER .EQ. 0) THEN
          R=B(I0)
         ELSE
          CALL DVSET(K+2,Z0,D(1),D(2))
          D(I+K+2-KK)=1
          DO 50 J=1,NDER
          A=K-J+1
          DO 50 L=1,K-J+2
          DIF=T(L+KK-1)-T(L+KK-K-2+J)
          D0=Z0
          IF(DIF .NE. 0) D0=A*(D(L+1)-D(L))/DIF
   50     D(L)=D0
          R=DVMPY(K-NDER+1,B(1),B(2),D(1),D(2))
         ENDIF
        ENDIF
       ENDIF
       DSPNB1=R
      ENDIF
      RETURN

  101 FORMAT(1X,A5,' =',I6,'   NOT IN RANGE')
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1996/04/01 15:02:26 mclareni
	6	* Mathlib gen
	7	*
	8	*
	9	#include "gen/pilot.h"
	10	FUNCTION DSPNB1(K,M,I,NDER,X,T,NERR)
	11
	12	#include "gen/imp64.inc"
	13	DIMENSION T(*),B(27),D(27)
	14	CHARACTER NAME()
	15	CHARACTER*80 ERRTXT
	16	PARAMETER (NAME = 'DSPNB1')
	17
	18	************************************************************************
	19	* NORBAS, VERSION: 15.03.1993
	20	************************************************************************
	21	*
	22	* DSPNB1 COMPUTES FUNCTION VALUES, VALUES OF DERIVATIVES, AND THE
	23	* VALUE OF INTEGRAL, RESPECTIVELY, OF NORMALIZED B-SPLINES
	24	* B(I,K)(X)
	25	* OF DEGREE K ( 0<= K <= 25 ) WITH INDEX I ( 1 <= I <= M-K-1 )
	26	* OVER A SET OF SPLINE-KNOTS
	27	* T(1),T(2), ... ,T(M) ( M >= 2*K+2 )
	28	* (KNOTS IN ASCENDING ORDER, WITH MULTIPLICITIES NOT GREATER
	29	* THAN K+1).
	30	*
	31	* THE FUNCTION VALUE OF THE NORMALIZED B-SPLINE B(I,K)(X) IS
	32	* IDENTICALLY ZERO OUTSIDE THE INTERVAL T(I) <= X < T(I+K+1).
	33	*
	34	* THE NORMALIZATION OF N(X) = B(I,K)(X) IS SUCH THAT THE INTGRAL OF
	35	* N(X) OVER THE WHOLE X-RANGE EQUALS
	36	* ( T(I+K+1) - T(I) ) / (K+1) .
	37	*
	38	* PARAMETERS:
	39	*
	40	* K (INTEGER) DEGREE (= ORDER - 1) OF B-SPLINES.
	41	* M (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES.
	42	* NDER (INTEGER) ON ENTRY, NDER MUST CONTAIN AN INTEGER VALUE .GE. -1
	43	* = -1: DSPNB1 COMPUTES THE INTEGRAL OF B(I,K)(TAU) OVER THE
	44	* RANGE TAU <= X.
	45	* = 0: DSPNB1 COMPUTES THE FUNCTION VALUE B(I,K)(X) FOR
	46	* FOR THE SPECIFIED VALUES OF I, K, AND X.
	47	* >= 1: DSPNB1 COMPUTES THE VALUE OF THE NDER-TH DERIVATIVE OF
	48	* B(I,K)(X) FOR THE SPECIFIED VALUES OF I, K, AND X
	49	* (IF NDER > K ZERO RETURNS).
	50	* X (DOUBLE PRECISION) ON ENTRY, X MUST CONTAIN THE VALUE OF THE
	51	* INDEPENDENT VARIABLE X OF B(I,K)(X)
	52	* T (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER M CONTAINING THE
	53	* KNOTS, ON ENTRY.
	54	* I (INTEGER) INDEX OF THE B-SPLINE B(I,K)(X)
	55	* NERR (INTEGER) ERROR INDICATOR. ON EXIT:
	56	* = 0: NO ERROR DETECTED
	57	* = 1: AT LEAST ONE OF THE CONSTANTS K , M , I , NDER IS ILLEGAL
	58	*
	59	* USAGE:
	60	*
	61	* THE FUNCTION-CALL
	62	* A = DSPNB1(K,M,I,NDER,X,T,NERR)
	63	* RETURNS
	64	* - THE VALUE OF THE INTEGRAL (NDER = -1) OR
65	* - THE FUNCTION VALUE (NDER = 0 ) OR
66	* - THE VALUE OF THE NDER-TH DERIVATIVE (NDER > 0 )
67	* OF THE NORMALIZED B-SPLINE B(I,K)(X) OF DEGREE K WITH INDEX I AT X.
68	*
69	* ERROR MESSAGES:
70	*
71	* IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
72	* PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
73	* K < 0 OR K > 25 OR
74	* M < 2*K+2 OR
75	* NDER < -1 OR
76	* I < 1 OR I > M-K-1
77	*
78	************************************************************************
79
80	PARAMETER (Z0 = 0 , Z1 = 1)
81
82	NERR=1
83	IF(K .LT. 0 .OR. K .GT. 25) THEN
84	WRITE(ERRTXT,101) 'K',K
85	CALL MTLPRT(NAME,'E210.1',ERRTXT)
86	ELSEIF(M .LT. 2*K+2) THEN
87	WRITE(ERRTXT,101) 'M',M
88	CALL MTLPRT(NAME,'E210.2',ERRTXT)
89	ELSEIF(I . LT. 1 .OR. I .GT. M-K-1) THEN
90	WRITE(ERRTXT,101) 'I',I
91	CALL MTLPRT(NAME,'E210.3',ERRTXT)
92	ELSEIF(NDER .LT. -1) THEN
93	WRITE(ERRTXT,101) 'NDER',NDER
94	CALL MTLPRT(NAME,'E210.5',ERRTXT)
95	ELSE
96
97	NERR=0
98	DSPNB1=Z0
99	IF( X .LT. T(I)
100	+ .OR. X .GT. T(I+K+1) .AND. NDER .GE. 0
101	+ .OR. K .LT. NDER ) RETURN
102
103	IF(NDER .EQ. -1) THEN
104	IF(X .GE. T(I+K+1)) THEN
105	R=(T(I+K+1)-T(I))/(K+1)
106	ELSE
107	KK=LKKSPL(X,T(I),MIN(2*(K+1),M-K-I))+I-1
108	IF(K .EQ. 0) THEN
109	R=X-T(KK-1)
110	ELSE
111	CALL DVSET(K+1,Z0,B(1),B(2))
112	B(KK-I)=1/(T(KK)-T(KK-1))
113	DO 10 L=1,K
114	DO 10 J=MAX(1,KK-I-L),MIN(K+1-L,KK-I)
115	DIF=T(I+J+L)-T(I+J-1)
116	B0=Z0
117	IF(DIF .NE. 0) B0=((X-T(I+J-1))B(J)+(T(I+J+L)-X)B(J+1))/DIF
118	10 B(J)=B0
119	S=Z0
120	DO 20 L=1,KK-I
121	20 S=S+(X-T(I+L-1))*B(L)
122	R=S*(T(I+K+1)-T(I))/(K+1)
123	ENDIF
124	ENDIF
125	DSPNB1=R
126	RETURN
127	ENDIF
128
129	IF(K .EQ. 0) THEN
130	R=Z1
131	ELSE
132	KK=LKKSPL(X,T(I),MIN(2*(K+1),M-K-I))+I-1
133	I0=I+K+2-KK
134	IF(I0 .EQ. 0) THEN
135	R=Z0
136	ELSE
137	E1=X-T(KK-1)
138	B(1)=Z1
139	DO 30 J=2,K-NDER+1
140	30 B(J)=E1*B(J-1)/(T(KK-2+J)-T(KK-1))
141	IF(KK .NE. I+1 .OR. NDER .NE. 0) THEN
142	E2=T(KK)-X
143	DO 40 J=1,K-NDER
144	E3=X-T(KK-1-J)
145	B(1)=E2*B(1)/(T(KK)-T(KK-J))
146	DO 40 L=2,K-NDER+1-J
147	40 B(L)=E3*B(L-1)/(T(KK-2+L)-T(KK-1-J))+
148	+ (T(KK-1+L)-X)*B(L)/(T(KK-1+L)-T(KK-J))
149	ENDIF
150	IF(NDER .EQ. 0) THEN
151	R=B(I0)
152	ELSE
153	CALL DVSET(K+2,Z0,D(1),D(2))
154	D(I+K+2-KK)=1
155	DO 50 J=1,NDER
156	A=K-J+1
157	DO 50 L=1,K-J+2
158	DIF=T(L+KK-1)-T(L+KK-K-2+J)
159	D0=Z0
160	IF(DIF .NE. 0) D0=A*(D(L+1)-D(L))/DIF
161	50 D(L)=D0
162	R=DVMPY(K-NDER+1,B(1),B(2),D(1),D(2))
163	ENDIF
164	ENDIF
165	ENDIF
166	DSPNB1=R
167	ENDIF
168	RETURN
169
170	101 FORMAT(1X,A5,' =',I6,' NOT IN RANGE')
171	END
172
173