]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 |