MINICERN/mathlib/gen/e/dspnb1.F

   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