5 * Revision 1.1.1.1 1996/04/01 15:02:25 mclareni
10 SUBROUTINE DSPIN1(K,N,XI,YI,KNOT,T,C,W,IW,NERR)
12 #include "gen/imp64.inc"
13 DIMENSION XI(*),YI(*),T(*),W(*),IW(*),C(*)
16 PARAMETER (NAME = 'DSPIN1')
18 ************************************************************************
19 * NORBAS, VERSION: 10.02.1993
20 ************************************************************************
22 * DSPIN1 COMPUTES THE COEFFICIENTS C(1),...,C(N) OF A POLYNOMIAL
23 * INTERPOLATION SPLINE S(X) IN B-SPLINE REPRESENTATION
25 * S(X) = SUMME(I=1,...,N) C(I) * B(I,K)(X)
27 * TO A USER SUPPLIED DATA SET
29 * (XI(J),YI(J)) , J = 1,2,...,N , N >= K+1
31 * OF A FUNCTION Y=F(X) , I.E.
33 * S(XI(J)) = YI(J) , J = 1,2,...,N .
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
38 * T(1),T(2), ... ,T(M) ( M = N+K+1 )
39 * (KNOTS IN ASCENDING ORDER, WITH MULTIPLICITIES NOT GREATER
41 * FOR FURTHER DETAILS TO THE ONE-DIMENSIONAL NORMALIZED B-SPLINES SEE
42 * THE COMMENTS TO DSPNB1.
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,
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.
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)) ,
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
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
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
92 * IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
93 * PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
97 ************************************************************************
99 PARAMETER (Z1 = 1 , Z2 = 2 , HALF = Z1/Z2)
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)
113 * COMPUTE KNOTS FROM INTERPOLATION POINTS (IF KNOT = 1 OR 2)
115 IF (KNOT .EQ. 1) THEN
116 CALL DSPKN1(K,M,XI(1),XI(N),T,NERR)
117 ELSEIF (KNOT .EQ. 2) THEN
122 20 T(I)=HALF*(XI(I-K-1)+XI(I))
125 * COMPUTE BAND-MATRIX W
130 IF (1 .LE. IJ .AND. IJ .LE. N)
131 + W((J-1)*L+I)=DSPNB1(K,M,J,0,XI(IJ),T,NERR)
134 * SOLVE SYSTEM OF EQUATIONS FOR COMPUTING C
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
148 101 FORMAT(1X,A5,' =',I6,' NOT IN RANGE')