MINICERN/mathlib/gen/e/dspin1.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:25  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE DSPIN1(K,N,XI,YI,KNOT,T,C,W,IW,NERR)
  11
  12 #include "gen/imp64.inc"
  13       DIMENSION XI(*),YI(*),T(*),W(*),IW(*),C(*)
  14       CHARACTER NAME*(*)
  15       CHARACTER*80 ERRTXT
  16       PARAMETER (NAME = 'DSPIN1')
  17
  18 ************************************************************************
  19 *   NORBAS, VERSION: 10.02.1993
  20 ************************************************************************
  21 *
  22 *   DSPIN1 COMPUTES THE COEFFICIENTS  C(1),...,C(N)  OF A POLYNOMIAL
  23 *   INTERPOLATION SPLINE  S(X)  IN B-SPLINE REPRESENTATION
  24 *
  25 *          S(X) =  SUMME(I=1,...,N)  C(I) * B(I,K)(X)
  26 *
  27 *   TO A USER SUPPLIED DATA SET
  28 *
  29 *          (XI(J),YI(J)) ,   J = 1,2,...,N  ,  N >= K+1
  30 *
  31 *   OF A FUNCTION  Y=F(X) , I.E.
  32 *
  33 *           S(XI(J)) = YI(J) , J = 1,2,...,N .
  34 *
  35 *   THE FUNCTIONS  B(I,K)(X)  ARE NORMALIZED B-SPLINES OF DEGREE  K
  36 *   ( K <= N-1  AND  0<= K <= 25)  WITH INDEX  I  (1 <= I <= N) OVER A
  37 *   SET OF SPLINE-KNOTS
  38 *              T(1),T(2), ... ,T(M)    ( M = N+K+1 )
  39 *   (KNOTS IN ASCENDING ORDER, WITH MULTIPLICITIES NOT GREATER
  40 *   THAN  K+1).
  41 *   FOR FURTHER DETAILS TO THE ONE-DIMENSIONAL NORMALIZED B-SPLINES SEE
  42 *   THE COMMENTS TO  DSPNB1.
  43 *
  44 *   PARAMETERS:
  45 *
  46 *   N     (INTEGER) NUMBER OF INTERPOLATION POINTS .
  47 *   K     (INTEGER) DEGREE (= ORDER - 1) OF B-SPLINES.
  48 *   XI    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER N .
  49 *         XI MUST CONTAIN THE INTERPOLATION POINTS IN ASCENDING ORDER,
  50 *         ON ENTRY.
  51 *   YI    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER N CONTAINING
  52 *         THE FUNCTION VALUES YI(J), J=1,...,N,  ON ENTRY.
  53 *   KNOT  (INTEGER) PARAMETER FOR STEERING THE CHOICE OF KNOTS.
  54 *         ON ENTRY:
  55 *         = 1 : KNOTS ARE COMPUTED BY  DSPIN1  IN THE FOLLOWING WAY:
  56 *               T(J) = XI(1) ,                J = 1,...,K+1
  57 *               T(J) = XI(1)+(J-K-1)*(XI(N)-XI(1)) ,
  58 *                                             J = K+2,...,N
  59 *               T(N+J) = XI(N) ,              J = 1,...,K+1
  60 *         = 2 : KNOTS ARE COMPUTED BY  DSPIN1  IN THE FOLLOWING WAY:
  61 *               T(J) = XI(1) ,                J = 1,...,K+1
  62 *               T(J) = (XI(J-K-1)+XI(J))/2 ,  J = K+2,...,N
  63 *               T(N+J) = XI(N) ,              J = 1,...,K+1
  64 *         OTHERWISE KNOTS ARE USER SUPPLIED. RECOMMENDED CHOICE :
  65 *               T(1) <= ... <= T(K+1) <= XI(1)
  66 *               XI(1) < T(K+2) < ... < T(N) < XI(N)
  67 *               XI(N) <= T(N+1) <= ... <= T(N+K+1)
  68 *   T     (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER M .
  69 *         IF THE INPUT VALUE OF THE PARAMETER KNOT IS  1 OR 2  THE
  70 *         KNOTS ARE COMPUTED BY  DSPIN1  AND THEY ARE GIVEN IN THE
  71 *         ARRAY T, ON EXIT.
  72 *         IN THE OTHER CASES THE ARRAY  T  MUST CONTAIN THE USER
  73 *         SUPPLIED KNOTS, ON ENTRY.
  74 *   W     (DOUBLE PECISION) WORKING ARRAY OF AT LEAST ORDER
  75 *         (3*K+1)*N .
  76 *   IW    (INTEGER) WORKING ARRAY OF AT LEAST ORDER  N .
  77 *   C     (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER  N . ON EXIT
  78 *         C(1),...,C(N)  CONTAIN THE COEFFICIENTS OF THE B-SPLINE
  79 *         REPRESENTATION OF  S(X).
  80 *   NERR  (INTEGER) ERROR INDICATOR. ON EXIT:
  81 *         = 0: NO ERROR DETECTED
  82 *         = 1: AT LEAST ONE OF THE CONSTANTS  K , N  IS ILLEGAL
  83 *         = 2: THE LAPACK ROUTINES  DGBTRF , DGBTRS  COULD NOT SOLVE
  84 *              THE LINEAR SYSTEM OF EQUATIONS WITH BAND-MATRIX FOR
  85 *              COMPUTING C(1),...,C(N) . IT INDICATES THAT A SOLUTION
  86 *              OF THE INTERPOLATION PROBLEM DOES NOT EXIST.
  87 *              (ESPECIALLY, THE EXISTENCE OF A SOLUTION DEPENDS ON THE
  88 *              SET OF KNOTS!)
  89 *
  90 *   ERROR MESSAGES:
  91 *
  92 *   IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
  93 *   PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
  94 *     K  < 0      OR    K > 25    OR
  95 *     N < K+1 .
  96 *
  97 ************************************************************************
  98
  99       PARAMETER (Z1 = 1 , Z2 = 2 , HALF = Z1/Z2)
 100
 101       NERR=1
 102       IF(K .LT. 0 .OR. K .GT. 25) THEN
 103        WRITE(ERRTXT,101) 'K',K
 104        CALL MTLPRT(NAME,'E210.1',ERRTXT)
 105       ELSEIF(N .LT. K+1) THEN
 106        WRITE(ERRTXT,101) 'N',N
 107        CALL MTLPRT(NAME,'E210.4',ERRTXT)
 108       ELSE
 109
 110        M=N+K+1
 111        L=3*K+1
 112 *
 113 *   COMPUTE KNOTS FROM INTERPOLATION POINTS (IF KNOT = 1 OR 2)
 114 *
 115        IF (KNOT .EQ. 1) THEN
 116         CALL DSPKN1(K,M,XI(1),XI(N),T,NERR)
 117        ELSEIF (KNOT .EQ. 2) THEN
 118         DO 10 I=1,K+1
 119         T(I)=XI(1)
 120    10   T(N+I)=XI(N)
 121         DO 20 I=K+2,N
 122    20   T(I)=HALF*(XI(I-K-1)+XI(I))
 123        ENDIF
 124 *
 125 *   COMPUTE BAND-MATRIX  W
 126 *
 127        DO 50 I=K+1,3*K+1
 128        DO 50 J=1,N
 129        IJ=I+J-2*K-1
 130        IF (1 .LE. IJ .AND. IJ .LE. N)
 131      +    W((J-1)*L+I)=DSPNB1(K,M,J,0,XI(IJ),T,NERR)
 132    50 CONTINUE
 133 *
 134 *   SOLVE SYSTEM OF EQUATIONS FOR COMPUTING  C
 135 *
 136        DO 60 J=1,N
 137    60  IW(J)=J
 138        NERR=2
 139        CALL DGBTRF(N,N,K,K,W,L,IW,INFO)
 140        IF(INFO .NE. 0) RETURN
 141        CALL DVCPY(N,YI(1),YI(2),C(1),C(2))
 142        CALL DGBTRS('N',N,K,K,1,W,L,IW,C,N,INFO)
 143        IF(INFO .NE. 0) RETURN
 144        NERR=0
 145       ENDIF
 146       RETURN
 147
 148   101 FORMAT(1X,A5,' =',I6,'   NOT IN RANGE')
 149       END
 150
 151
 152