]>
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 DSPIN2(KX,KY,NX,NY,XI,YI,ZI,NDIMZ,KNOT,TX,TY,C,NDIMC, | |
11 | + W,IW,NERR) | |
12 | ||
13 | #include "gen/imp64.inc" | |
14 | DIMENSION XI(*),YI(*),ZI(NDIMZ,*),TX(*),TY(*) | |
15 | DIMENSION W(*),IW(*),C(NDIMC,*) | |
16 | CHARACTER NAME*(*) | |
17 | CHARACTER*80 ERRTXT | |
18 | PARAMETER (NAME = 'DSPIN2') | |
19 | ||
20 | ||
21 | ************************************************************************ | |
22 | * NORBAS, VERSION: 10.02.1993 | |
23 | ************************************************************************ | |
24 | * | |
25 | * DSPIN2 COMPUTES THE COEFFICIENTS | |
26 | * C(I,J) (I=1,...,NX , J=1,...,NY) | |
27 | * OF A TWO-DIMENSIONAL POLYNOMIAL INTERPOLATION SPLINE Z = S(X,Y) IN | |
28 | * REPRESENTATION OF NORMALIZED TWO-DIMENSIONAL B-SPLINES B(I,J)(X,Y) | |
29 | * | |
30 | * S(X,Y) = SUMME(I=1,...,NX) | |
31 | * SUMME(J=1,...,NY) C(I,J) * B(I,J)(X,Y) | |
32 | * | |
33 | * TO A USER SUPPLIED DATA SET | |
34 | * | |
35 | * (XI(I),YI(J),ZI(I,J)) (I=1,...,NX , J=1,...,NY) | |
36 | * | |
37 | * OF A FUNCTION Z = F(X,Y) , I.E. | |
38 | * | |
39 | * S(XI(I),YI(J)) = Z(I,J) (I=1,...,NX , J=1,...,NY) . | |
40 | * | |
41 | * THE TWO-DIMENSIONAL B-SPLINES B(I,J)(X,Y) ARE THE PRODUCT OF TWO | |
42 | * ONE-DIMENSIONAL B-SPLINES BX , BY | |
43 | * B(I,J)(X,Y) = BX(I,KX)(X) * BY(J,KY)(Y) | |
44 | * OF DEGREE KX AND KY ( 0 <= KX , KY <= 25 ) WITH INDICES I , J | |
45 | * ( 1 <= I <= NX , 1 <= J <= NY ) OVER TWO SETS OF SPLINE-KNOTS | |
46 | * TX(1),TX(2),...,TX(MX) ( MX = NX+KX+1 ) | |
47 | * TY(1),TY(2),...,TY(MY) ( MY = NY+KY+1 ) , | |
48 | * RESPECTIVELY. | |
49 | * FOR FURTHER DETAILS TO THE ONE- AND TWO-DIMENSIONAL NORMALIZED | |
50 | * B-SPLINES SEE THE COMMENTS TO DSPNB1 AND DSPNB2. | |
51 | * | |
52 | * PARAMETERS: | |
53 | * | |
54 | * NX (INTEGER) NUMBER OF INTERPOLATION POINTS IN X-DIRECTION : | |
55 | * XI(I) , I=1,...,NX . | |
56 | * NY (INTEGER) NUMBER OF INTERPOLATION POINTS IN Y-DIRECTION : | |
57 | * YI(J) , J=1,...,NY . | |
58 | * KX (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN X-DIRECTION | |
59 | * OVER THE SET OF KNOTS TX, | |
60 | * WITH KX <= NX-1 AND 0 <= KX <= 25 . | |
61 | * KY (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN Y-DIRECTION | |
62 | * OVER THE SET OF KNOTS TY, | |
63 | * WITH KY <= NY-1 AND 0 <= KY <= 25 . | |
64 | * NDIMC (INTEGER) DECLARED FIRST DIMENSION OF ARRAY C IN THE | |
65 | * CALLING PROGRAM, WITH NDIMC >= NX . | |
66 | * NDIMZ (INTEGER) DECLARED FIRST DIMENSION OF ARRAY ZI IN THE | |
67 | * CALLING PROGRAM, WITH NDIMZ >= NX . | |
68 | * XI (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER NX . | |
69 | * XI MUST CONTAIN THE INTERPOLATION POINTS IN X-DIRECTION IN | |
70 | * ASCENDING ORDER, ON ENTRY. | |
71 | * YI (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER NY . | |
72 | * YI MUST CONTAIN THE INTERPOLATION POINTS IN Y-DIRECTION IN | |
73 | * ASCENDING ORDER, ON ENTRY. | |
74 | * ZI (DOUBLE PRECISION) ARRAY OF ORDER (NDIMZ , >= NY) . | |
75 | * ON ENTRY ZI MUST CONTAIN THE GIVEN FUNCTION VALUES | |
76 | * Z(I,J) AT THE INTERPOLATION POINTS (X(I),Y(J)) | |
77 | * ( I=1,...,NX , J=1,...,NY ). | |
78 | * KNOT (INTEGER) PARAMETER FOR STEERING THE CHOICE OF KNOTS. | |
79 | * ON ENTRY: | |
80 | * = 1 : KNOTS ARE COMPUTED BY DSPIN2 IN THE FOLLOWING WAY: | |
81 | * TX(J) = XI(1) , J = 1,...,KX+1 | |
82 | * TX(J) = XI(1)+(J-KX-1)*(XI(NX)-XI(1)) , | |
83 | * J = KX+2,...,NX | |
84 | * TX(N+J) = XI(NX) , J = 1,...,KX+1 | |
85 | * = 2 : KNOTS ARE COMPUTED BY DSPIN2 IN THE FOLLOWING WAY: | |
86 | * TX(J) = XI(1) , J = 1,...,KX+1 | |
87 | * TX(J) = (XI(J-KX-1)+XI(J))/2 , J = KX+2,...,NX | |
88 | * TX(N+J) = XI(NX) , J = 1,...,KX+1 | |
89 | * OTHERWISE KNOTS ARE USER SUPPLIED. RECOMMENDED CHOICE : | |
90 | * T(1) <= ... <= T(K+1) <= XI(1) | |
91 | * XI(1) < T(K+2) < ... < T(N) < XI(N) | |
92 | * XI(N) <= T(N+1) <= ... <= T(N+K+1) | |
93 | * OTHERWISE KNOTS ARE USER SUPPLIED. RECOMMENDED CHOICE : | |
94 | * TX(1) <= ... <= TX(KX+1) <= XI(1) | |
95 | * XI(1) < TX(KX+2) < ... < TX(NX) < XI(NX) | |
96 | * XI(NX) <= TX(NX+1) <= ... <= TX(NX+KX+1) | |
97 | * IN ALL CASES THE SAME CHOICE IS USED FOR KNOTS TY IN | |
98 | * Y-DIRECTION . | |
99 | * TX (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MX. | |
100 | * IF THE INPUT VALUE OF THE PARAMETER KNOT IS 1 OR 2 THE | |
101 | * KNOTS ARE COMPUTED BY DSPIN2 AND THEY ARE GIVEN IN THE | |
102 | * ARRAY TX, ON EXIT. | |
103 | * IN THE OTHER CASES THE ARRAY TX MUST CONTAIN THE USER | |
104 | * SUPPLIED KNOTS IN X-DIRECTION, ON ENTRY. | |
105 | * TY (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MY. | |
106 | * IF THE INPUT VALUE OF THE PARAMETER KNOT IS 1 OR 2 THE | |
107 | * KNOTS ARE COMPUTED BY DSPIN2 AND THEY ARE GIVEN IN THE | |
108 | * ARRAY TY, ON EXIT. | |
109 | * IN THE OTHER CASES THE ARRAY TY MUST CONTAIN THE USER | |
110 | * SUPPLIED KNOTS IN Y-DIRECTION, ON ENTRY. | |
111 | * W (DOUBLE PRECISION) WORKING ARRAY OF AT LEAST ORDER (L+1)*NINT, | |
112 | * WITH NINT=NX*NY , K=KX*NY , L=3*K+1 . | |
113 | * IW (INTEGER) WORKING ARRAY OF AT LEAST ORDER NX*NY. | |
114 | * C (DOUBLE PRECISION) ARRAY OF ORDER (NDIMC, >= NY). | |
115 | * ON EXIT C(I,J) CONTAINS THE (I,J)-TH COEFFICIENT OF THE | |
116 | * TWO-DIMENSIONAL B-SPLINE REPRESENTATION OF S(X,Y) . | |
117 | * NERR (INTEGER) ERROR INDICATOR. ON EXIT: | |
118 | * = 0: NO ERROR DETECTED | |
119 | * = 1: AT LEAST ONE OF THE CONSTANTS KX , KY , NX , NY | |
120 | * IS ILLEGAL | |
121 | * = 2: THE LAPACK ROUTINES DGBTRF , DGBTRS COULD NOT SOLVE | |
122 | * THE LINEAR SYSTEM OF EQUATIONS WITH BAND-MATRIX FOR | |
123 | * COMPUTING C(1,1),...,C(NX,NY) . IT INDICATES THAT | |
124 | * A SOLUTION OF THE INTERPOLATION PROBLEM DOES NOT EXIST. | |
125 | * (ESPECIALLY, THE EXISTENCE OF A SOLUTION DEPENDS ON THE | |
126 | * SET OF KNOTS!) | |
127 | * | |
128 | * ERROR MESSAGES: | |
129 | * | |
130 | * IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT- | |
131 | * PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED: | |
132 | * KX < 0 OR KX > 25 OR | |
133 | * KY < 0 OR KY > 25 OR | |
134 | * NX < KX+1 OR NY < KY+1 . | |
135 | * | |
136 | ************************************************************************ | |
137 | ||
138 | ||
139 | PARAMETER (Z1 = 1 , Z2 = 2 , HALF = Z1/Z2) | |
140 | ||
141 | NERR=1 | |
142 | IF(KX .LT. 0 .OR. KX .GT. 25) THEN | |
143 | WRITE(ERRTXT,101) 'KX',KX | |
144 | CALL MTLPRT(NAME,'E210.1',ERRTXT) | |
145 | ELSEIF(KY .LT. 0 .OR. KY .GT. 25) THEN | |
146 | WRITE(ERRTXT,101) 'KY',KY | |
147 | CALL MTLPRT(NAME,'E210.1',ERRTXT) | |
148 | ELSEIF(NX .LT. KX+1) THEN | |
149 | WRITE(ERRTXT,101) 'NX',NX | |
150 | CALL MTLPRT(NAME,'E210.4',ERRTXT) | |
151 | ELSEIF(NY .LT. KY+1) THEN | |
152 | WRITE(ERRTXT,101) 'NY',NY | |
153 | CALL MTLPRT(NAME,'E210.4',ERRTXT) | |
154 | ELSE | |
155 | ||
156 | K=KX*NY | |
157 | L=3*K+1 | |
158 | NINT=NX*NY | |
159 | MX=NX+KX+1 | |
160 | MY=NY+KY+1 | |
161 | * | |
162 | * COMPUTE KNOTS FROM INTERPOLATION POINTS ( IF KNOT=1 OR 2 ) | |
163 | * | |
164 | IF (KNOT .EQ. 1) THEN | |
165 | CALL DSPKN2(KX,KY,MX,MY,XI(1),XI(NX),YI(1),YI(NY),TX,TY,NERR) | |
166 | ELSEIF (KNOT .EQ. 2) THEN | |
167 | DO 10 I=1,KX+1 | |
168 | TX(I)=XI(1) | |
169 | 10 TX(NX+I)=XI(NX) | |
170 | DO 20 I=KX+2,NX | |
171 | 20 TX(I)=HALF*(XI(I-KX-1)+XI(I)) | |
172 | DO 30 J=1,KY+1 | |
173 | TY(J)=YI(1) | |
174 | 30 TY(NY+J)=YI(NY) | |
175 | DO 40 J=KY+2,NY | |
176 | 40 TY(J)=HALF*(YI(J-KY-1)+YI(J)) | |
177 | ENDIF | |
178 | * | |
179 | * COMPUTE BAND-MATRIX W | |
180 | * | |
181 | * NUMBER OF SUB-, SUPER-DIAGONALS: KL=KU=(KX-1)*NY+KY-1<=KX*NY | |
182 | * | |
183 | DO 90 IC=1,NX | |
184 | X=XI(IC) | |
185 | DO 90 JC=1,NY | |
186 | Y=YI(JC) | |
187 | IJ=(IC-1)*NX+JC | |
188 | DO 90 I=1,NX | |
189 | DO 90 J=1,NY | |
190 | JZ=(I-1)*NX+J | |
191 | IZ=IJ+2*K+1-JZ | |
192 | IF (K+1 .LE. IZ .AND. IZ .LE. 3*K+1) | |
193 | + W((JZ-1)*L+IZ)=DSPNB2(KX,KY,MX,MY,I,J,0,0,X,Y,TX,TY,NERR) | |
194 | 90 CONTINUE | |
195 | * | |
196 | * SOLVE SYSTEM OF EQUATIONS FOR COMPUTING C | |
197 | * | |
198 | DO 110 J=1,NINT | |
199 | 110 IW(J)=J | |
200 | NERR=2 | |
201 | CALL DGBTRF(NINT,NINT,K,K,W,L,IW,INFO) | |
202 | IF(INFO .NE. 0) RETURN | |
203 | LW=L*NINT | |
204 | DO 120 I=1,NX | |
205 | DO 120 J=1,NY | |
206 | 120 W(LW+(I-1)*NY+J)=ZI(I,J) | |
207 | CALL DGBTRS('N',NINT,K,K,1,W,L,IW,W(LW+1),NINT,INFO) | |
208 | IF(INFO .NE. 0) RETURN | |
209 | DO 130 I=1,NX | |
210 | DO 130 J=1,NY | |
211 | 130 C(I,J)=W(LW+(I-1)*NY+J) | |
212 | NERR=0 | |
213 | ENDIF | |
214 | ||
215 | RETURN | |
216 | ||
217 | 101 FORMAT(1X,A5,' =',I6,' NOT IN RANGE') | |
218 | END | |
219 | ||
220 | ||
221 |