Line data Source code
1 : ! CS2TST
2 : ! 11/20/98
3 :
4 : ! This program computes interpolation errors using the
5 : ! scattered data package CSHEP2D for each of ten test
6 : ! functions and a 33 by 33 uniform grid of interpolation
7 : ! points in the unit square.
8 :
9 : ! This program uses subroutines TESTDT and TSTFN1 from
10 : ! ACM Algorithm SURVEY to generate a node set and and the
11 : ! test function values.
12 :
13 0 : subroutine test_renka_db
14 : use interp_2d_lib_db
15 :
16 : integer :: NMAX, NRMAX, NI
17 : parameter (NMAX=100, NRMAX=10, NI=33)
18 :
19 : ! Array storage:
20 :
21 0 : double precision :: X(NMAX), Y(NMAX), W(NMAX), RW(NMAX), A(9,NMAX), P(NI), FT(NI,NI)
22 : integer :: LCELL(NRMAX,NRMAX), LNEXT(NMAX)
23 :
24 0 : double precision :: DEL,DUM,DX,DY,ERMAX,ERMEAN,PW,RMAX,SSA,SSE,SSM,SUM,XMIN,YMIN
25 : ! double precision CS2VAL_db
26 : integer :: I, IER, J, K, KF, KFF, KFL, KS, LOUT, N, NC, NFUN, NP, NR, NSET, NW
27 :
28 : ! Data:
29 :
30 : ! LOUT = Logical unit number for output.
31 : ! NSET = Number of node sets.
32 : ! NFUN = Number of test functions.
33 :
34 : DATA LOUT/6/, NSET/5/, NFUN/10/
35 :
36 : ! Specify test parameters
37 :
38 0 : KS = 2
39 0 : KFF = 1
40 0 : KFL = 10
41 0 : NC = 17
42 0 : NW = 30
43 0 : NR = 6
44 0 : call testdt_db (KS, N,X,Y)
45 :
46 : ! Set up uniform grid points.
47 :
48 0 : DEL = 1./dble(NI-1)
49 0 : do I = 1,NI
50 0 : P(I) = dble(I-1)*DEL
51 : end do
52 :
53 : ! Initialize the average SSE/SSM value to zero.
54 :
55 0 : SSA = 0.
56 :
57 : ! Print a heading and loop on test functions.
58 :
59 0 : write (LOUT,200) KS, N, NI, NC, NW, NR
60 0 : do KF = KFF,KFL
61 :
62 : ! Compute true function values at the nodes.
63 :
64 0 : do K = 1,N
65 0 : call TSTFN1_db (KF,X(K),Y(K),0, W(K),DUM,DUM)
66 : end do
67 :
68 : ! Compute true function values FT on the uniform grid, and
69 : ! accumulate the sum of values SUM and sum of squared
70 : ! values SSM.
71 :
72 0 : SUM = 0.
73 0 : SSM = 0.
74 0 : do I = 1,NI
75 0 : do J = 1,NI
76 0 : call TSTFN1_db (KF,P(I),P(J),0, FT(I,J),DUM,DUM)
77 0 : SUM = SUM + FT(I,J)
78 0 : SSM = SSM + FT(I,J)**2
79 : end do
80 : end do
81 :
82 : ! Compute the sum of squared deviations from the mean SSM.
83 :
84 0 : SSM = SSM - SUM*SUM/dble(NI*NI)
85 :
86 : ! Compute parameters A and RW defining the interpolant.
87 :
88 0 : call interp_CSHEP2_db (N,X,Y,W,NC,NW,NR, LCELL,LNEXT,XMIN,YMIN,DX,DY,RMAX,RW,A,IER)
89 0 : if (IER /= 0) GOTO 21
90 :
91 : ! Compute interpolation errors.
92 :
93 0 : ERMEAN = 0.
94 0 : ERMAX = 0.
95 0 : SSE = 0.
96 0 : do I = 1,NI
97 0 : do J = 1,NI
98 0 : IER = 0
99 0 : PW = interp_CS2VAL_db (P(I),P(J),N,X,Y,W,NR,LCELL,LNEXT,XMIN,YMIN,DX,DY,RMAX,RW,A,IER) -FT(I,J)
100 0 : if (IER /= 0) then
101 0 : write(LOUT,*) 'IER nonzero from CS2VAL_db'
102 0 : stop 1
103 : end IF
104 0 : ERMEAN = ERMEAN + ABS(PW)
105 0 : ERMAX = MAX(ERMAX,ABS(PW))
106 0 : SSE = SSE + PW*PW
107 : end do
108 : end do
109 0 : NP = NI*NI
110 0 : ERMEAN = ERMEAN/dble(NP)
111 0 : SSE = SSE/SSM
112 0 : SSA = SSA + SSE
113 0 : write (LOUT,210) KF, ERMAX, ERMEAN, SSE
114 : end do
115 :
116 : ! Print the average SSE/SSM value (averaged over the test
117 : ! functions).
118 :
119 0 : return
120 :
121 : ! N is outside its valid range.
122 :
123 : 20 write (LOUT,500) N, NMAX
124 : stop 1
125 :
126 : ! Error in CSHEP2.
127 :
128 0 : 21 if (IER == 2) write (LOUT,510)
129 0 : if (IER == 3) write (LOUT,520)
130 0 : stop 1
131 :
132 : ! Print formats:
133 :
134 : 200 format ('RENKA790_db: Node set ',I2,4X,'N =',I4,4X,'NI = ',I2,4X,'NC = ',I2,4X,'NW = ',I2,4X,'NR = ',I2/1X,16X,
135 : & 'Function',4X,'Max Error',4X,'Mean Error',4X,'SSE/SSM')
136 : 210 format (1X,19X,I2,9X,F7.4,6X,F8.5,2X,F9.6)
137 :
138 : ! Error message formats:
139 :
140 : 500 format (///1X,10X,'*** Error in data -- N = ',I4,', Maximum value =',I4,' ***')
141 : 510 format (///1X,14X,'*** Error in CSHEP2 -- duplicate nodes encountered ***')
142 : 520 format (///1X,14X,'*** Error in CSHEP2 -- all nodes are collinear ***')
143 : end subroutine test_renka_db
144 :
145 :
146 0 : subroutine testdt_db (K, N,X,Y)
147 : double precision :: X(100), Y(100)
148 : integer :: K, N
149 :
150 : ! **********************************************************
151 :
152 : ! Robert J. Renka
153 : ! Dept. of Computer Science
154 : ! Univ. of North Texas
155 : ! renka@cs.unt.edu
156 : ! 03/28/97
157 :
158 : ! This subroutine returns one of five sets of nodes used
159 : ! for testing scattered data fitting methods. All five sets
160 : ! approximately cover the unit square [0,1] X [0,1]: the
161 : ! convex hulls of sets 1 and 3 extend slightly outside the
162 : ! square but do not completely cover it, those of sets 2 and
163 : ! 5 coincide with the unit square, and the convex hull of
164 : ! set 4 is a large subset of the unit square.
165 :
166 : ! On input:
167 :
168 : ! K = integer in the range 1 to 5 which determines the
169 : ! choice of data set as follows:
170 :
171 : ! K = 1 - Franke's 100-node set
172 : ! K = 2 - Franke's 33-node set
173 : ! K = 3 - Lawson's 25-node set
174 : ! K = 4 - Random 100-node set
175 : ! K = 5 - Gridded 81-node set
176 :
177 : ! X,Y = Arrays of length at least N(K), where
178 : ! N(1) = 100, N(2) = 33, N(3) = 25,
179 : ! N(4) = 100, and N(5) = 81.
180 :
181 : ! Input parameters are not altered by this routine.
182 :
183 : ! On output:
184 :
185 : ! N = Number of nodes in set K, or 0 if K is outside
186 : ! its valid range.
187 :
188 : ! X,Y = Nodal coordinates of node set K.
189 :
190 : ! Subprograms required by TESTDT: None
191 :
192 : ! **********************************************************
193 :
194 : double precision :: X1(100),Y1(100),X2(33),Y2(33),X3(25),Y3(25),X4(100),Y4(100),X5(81),Y5(81)
195 : integer :: I
196 :
197 : ! Node set 1: Franke's 100-node set.
198 :
199 : DATA (X1(I),Y1(I), I = 1,20)/
200 : & 0.0227035, -0.0310206, 0.0539888, 0.1586742,
201 : & 0.0217008, 0.2576924, 0.0175129, 0.3414014,
202 : & 0.0019029, 0.4943596, -0.0509685, 0.5782854,
203 : & 0.0395408, 0.6993418, -0.0487061, 0.7470194,
204 : & 0.0315828, 0.9107649, -0.0418785, 0.9962890,
205 : & 0.1324189, 0.0501330, 0.1090271, 0.0918555,
206 : & 0.1254439, 0.2592973, 0.0934540, 0.3381592,
207 : & 0.0767578, 0.4171125, 0.1451874, 0.5615563,
208 : & 0.0626494, 0.6552235, 0.1452734, 0.7524066,
209 : & 0.0958668, 0.9146523, 0.0695559, 0.9632421/
210 : DATA (X1(I),Y1(I), I = 21,40)/
211 : & 0.2645602, 0.0292939, 0.2391645, 0.0602303,
212 : & 0.2088990, 0.2668783, 0.2767329, 0.3696044,
213 : & 0.1714726, 0.4801738, 0.2266781, 0.5940595,
214 : & 0.1909212, 0.6878797, 0.1867647, 0.8185576,
215 : & 0.2304634, 0.9046507, 0.2426219, 0.9805412,
216 : & 0.3663168, 0.0396955, 0.3857662, 0.0684484,
217 : & 0.3832392, 0.2389548, 0.3179087, 0.3124129,
218 : & 0.3466321, 0.4902989, 0.3776591, 0.5199303,
219 : & 0.3873159, 0.6445227, 0.3812917, 0.8203789,
220 : & 0.3795364, 0.8938079, 0.2803515, 0.9711719/
221 : DATA (X1(I),Y1(I), I = 41,60)/
222 : & 0.4149771, -0.0284618, 0.4277679, 0.1560965,
223 : & 0.4200010, 0.2262471, 0.4663631, 0.3175094,
224 : & 0.4855658, 0.3891417, 0.4092026, 0.5084949,
225 : & 0.4792578, 0.6324247, 0.4812279, 0.7511007,
226 : & 0.3977761, 0.8489712, 0.4027321, 0.9978728,
227 : & 0.5848691, -0.0271948, 0.5730076, 0.1272430,
228 : & 0.6063893, 0.2709269, 0.5013894, 0.3477728,
229 : & 0.5741311, 0.4259422, 0.6106955, 0.6084711,
230 : & 0.5990105, 0.6733781, 0.5380621, 0.7235242,
231 : & 0.6096967, 0.9242411, 0.5026188, 1.0308762/
232 : DATA (X1(I),Y1(I), I = 61,80)/
233 : & 0.6616928, 0.0255959, 0.6427836, 0.0707835,
234 : & 0.6396475, 0.2008336, 0.6703963, 0.3259843,
235 : & 0.7001181, 0.4890704, 0.6333590, 0.5096324,
236 : & 0.6908947, 0.6697880, 0.6895638, 0.7759569,
237 : & 0.6718889, 0.9366096, 0.6837675, 1.0064516,
238 : & 0.7736939, 0.0285374, 0.7635332, 0.1021403,
239 : & 0.7410424, 0.1936581, 0.8258981, 0.3235775,
240 : & 0.7306034, 0.4714228, 0.8086609, 0.6091595,
241 : & 0.8214531, 0.6685053, 0.7290640, 0.8022808,
242 : & 0.8076643, 0.8476790, 0.8170951, 1.0512371/
243 : DATA (X1(I),Y1(I), I = 81,100)/
244 : & 0.8424572, 0.0380499, 0.8684053, 0.0902048,
245 : & 0.8366923, 0.2083092, 0.9418461, 0.3318491,
246 : & 0.8478122, 0.4335632, 0.8599583, 0.5910139,
247 : & 0.9175700, 0.6307383, 0.8596328, 0.8144841,
248 : & 0.9279871, 0.9042310, 0.8512805, 0.9696030,
249 : & 1.0449820, -0.0120900, 0.9670631, 0.1334114,
250 : & 0.9857884, 0.2695844, 0.9676313, 0.3795281,
251 : & 1.0129299, 0.4396054, 0.9657040, 0.5044425,
252 : & 1.0019855, 0.6941519, 1.0359297, 0.7459923,
253 : & 1.0414677, 0.8682081, 0.9471506, 0.9801409/
254 :
255 : ! Node set 2: Franke's 33-node set.
256 :
257 : DATA (X2(I),Y2(I), I = 1,33)/
258 : & 0.05, 0.45, 0.00, 0.50,
259 : & 0.00, 1.00, 0.00, 0.00,
260 : & 0.10, 0.15, 0.10, 0.75,
261 : & 0.15, 0.30, 0.20, 0.10,
262 : & 0.25, 0.20, 0.30, 0.35,
263 : & 0.35, 0.85, 0.50, 0.00,
264 : & 0.50, 1.00, 0.55, 0.95,
265 : & 0.60, 0.25, 0.60, 0.65,
266 : & 0.60, 0.85, 0.65, 0.70,
267 : & 0.70, 0.20, 0.70, 0.65,
268 : & 0.70, 0.90, 0.75, 0.10,
269 : & 0.75, 0.35, 0.75, 0.85,
270 : & 0.80, 0.40, 0.80, 0.65,
271 : & 0.85, 0.25, 0.90, 0.35,
272 : & 0.90, 0.80, 0.95, 0.90,
273 : & 1.00, 0.00, 1.00, 0.50,
274 : & 1.00, 1.00/
275 :
276 : ! Node set 3: Lawson's 25-node set.
277 :
278 : DATA (X3(I),Y3(I), I = 1,25)/
279 : & 0.13750, 0.97500, 0.91250, 0.98750,
280 : & 0.71250, 0.76250, 0.22500, 0.83750,
281 : & -0.05000, 0.41250, 0.47500, 0.63750,
282 : & 0.05000, -0.05000, 0.45000, 1.03750,
283 : & 1.08750, 0.55000, 0.53750, 0.80000,
284 : & -0.03750, 0.75000, 0.18750, 0.57500,
285 : & 0.71250, 0.55000, 0.85000, 0.43750,
286 : & 0.70000, 0.31250, 0.27500, 0.42500,
287 : & 0.45000, 0.28750, 0.81250, 0.18750,
288 : & 0.45000, -0.03750, 1.00000, 0.26250,
289 : & 0.50000, 0.46250, 0.18750, 0.26250,
290 : & 0.58750, 0.12500, 1.05000, -0.06125,
291 : & 0.10000, 0.11250/
292 :
293 : ! Node set 4: Random 100-node set.
294 :
295 : DATA (X4(I),Y4(I), I = 1,20)/
296 : & 0.0096326, 0.3083158, 0.0216348, 0.2450434,
297 : & 0.0298360, 0.8613847, 0.0417447, 0.0977864,
298 : & 0.0470462, 0.3648355, 0.0562965, 0.7156339,
299 : & 0.0646857, 0.5311312, 0.0740377, 0.9755672,
300 : & 0.0873907, 0.1781117, 0.0934832, 0.5452797,
301 : & 0.1032216, 0.1603881, 0.1110176, 0.7837139,
302 : & 0.1181193, 0.9982015, 0.1251704, 0.6910589,
303 : & 0.1327330, 0.1049580, 0.1439536, 0.8184662,
304 : & 0.1564861, 0.7086405, 0.1651043, 0.4456593,
305 : & 0.1786039, 0.1178342, 0.1886405, 0.3189021/
306 : DATA (X4(I),Y4(I), I = 21,40)/
307 : & 0.2016706, 0.9668446, 0.2099886, 0.7571834,
308 : & 0.2147003, 0.2016598, 0.2204141, 0.3232444,
309 : & 0.2343715, 0.4368583, 0.2409660, 0.8907869,
310 : & 0.2527740, 0.0647260, 0.2570839, 0.5692618,
311 : & 0.2733365, 0.2947027, 0.2853833, 0.4332426,
312 : & 0.2901755, 0.3347464, 0.2964854, 0.7436284,
313 : & 0.3019725, 0.1066265, 0.3125695, 0.8845357,
314 : & 0.3307163, 0.5158730, 0.3378504, 0.9425637,
315 : & 0.3439061, 0.4799701, 0.3529922, 0.1783069,
316 : & 0.3635507, 0.1146760, 0.3766172, 0.8225797/
317 : DATA (X4(I),Y4(I), I = 41,60)/
318 : & 0.3822429, 0.2270688, 0.3869838, 0.4073598,
319 : & 0.3973137, 0.8875080, 0.4170708, 0.7631616,
320 : & 0.4255588, 0.9972804, 0.4299218, 0.4959884,
321 : & 0.4372839, 0.3410421, 0.4705033, 0.2498120,
322 : & 0.4736655, 0.6409007, 0.4879299, 0.1058690,
323 : & 0.4940260, 0.5411969, 0.5055324, 0.0089792,
324 : & 0.5162593, 0.8784268, 0.5219219, 0.5515874,
325 : & 0.5348529, 0.4038952, 0.5483213, 0.1654023,
326 : & 0.5569571, 0.2965158, 0.5638611, 0.3660356,
327 : & 0.5784908, 0.0366554, 0.5863950, 0.9502420/
328 : DATA (X4(I),Y4(I), I = 61,80)/
329 : & 0.5929148, 0.2638101, 0.5987839, 0.9277386,
330 : & 0.6117561, 0.5377694, 0.6252296, 0.7374676,
331 : & 0.6331381, 0.4674627, 0.6399048, 0.9186109,
332 : & 0.6488972, 0.0416884, 0.6558537, 0.1291029,
333 : & 0.6677405, 0.6763676, 0.6814074, 0.8444238,
334 : & 0.6887812, 0.3273328, 0.6940896, 0.1893879,
335 : & 0.7061687, 0.0645923, 0.7160957, 0.0180147,
336 : & 0.7317445, 0.8904992, 0.7370798, 0.4160648,
337 : & 0.7462030, 0.4688995, 0.7566957, 0.2174508,
338 : & 0.7699998, 0.5734231, 0.7879347, 0.8853319/
339 : DATA (X4(I),Y4(I), I = 81,100)/
340 : & 0.7944014, 0.8018436, 0.8164468, 0.6388941,
341 : & 0.8192794, 0.8931002, 0.8368405, 0.1000558,
342 : & 0.8500993, 0.2789506, 0.8588255, 0.9082948,
343 : & 0.8646496, 0.3259159, 0.8792329, 0.8318747,
344 : & 0.8837536, 0.0508513, 0.8900077, 0.9708450,
345 : & 0.8969894, 0.5120548, 0.9044917, 0.2859716,
346 : & 0.9083947, 0.9581641, 0.9203972, 0.6183429,
347 : & 0.9347906, 0.3779934, 0.9434519, 0.4010423,
348 : & 0.9490328, 0.9478657, 0.9569571, 0.7425486,
349 : & 0.9772067, 0.8883287, 0.9983493, 0.5496750/
350 :
351 : ! Node set 5: 9 by 9 uniform grid.
352 :
353 : DATA (X5(I),Y5(I), I = 1,20)/
354 : & 0.125, 0.000, 0.000, 0.125,
355 : & 0.000, 0.250, 0.000, 0.375,
356 : & 0.000, 0.500, 0.000, 0.625,
357 : & 0.000, 0.750, 0.000, 0.875,
358 : & 0.000, 1.000, 0.000, 0.000,
359 : & 0.125, 0.125, 0.125, 0.250,
360 : & 0.125, 0.375, 0.125, 0.500,
361 : & 0.125, 0.625, 0.125, 0.750,
362 : & 0.125, 0.875, 0.125, 1.000,
363 : & 0.250, 0.000, 0.250, 0.125/
364 : DATA (X5(I),Y5(I), I = 21,40)/
365 : & 0.250, 0.250, 0.250, 0.375,
366 : & 0.250, 0.500, 0.250, 0.625,
367 : & 0.250, 0.750, 0.250, 0.875,
368 : & 0.250, 1.000, 0.375, 0.000,
369 : & 0.375, 0.125, 0.375, 0.250,
370 : & 0.375, 0.375, 0.375, 0.500,
371 : & 0.375, 0.625, 0.375, 0.750,
372 : & 0.375, 0.875, 0.375, 1.000,
373 : & 0.500, 0.000, 0.500, 0.125,
374 : & 0.500, 0.250, 0.500, 0.375/
375 : DATA (X5(I),Y5(I), I = 41,60)/
376 : & 0.500, 0.500, 0.500, 0.625,
377 : & 0.500, 0.750, 0.500, 0.875,
378 : & 0.500, 1.000, 0.625, 0.000,
379 : & 0.625, 0.125, 0.625, 0.250,
380 : & 0.625, 0.375, 0.625, 0.500,
381 : & 0.625, 0.625, 0.625, 0.750,
382 : & 0.625, 0.875, 0.625, 1.000,
383 : & 0.750, 0.000, 0.750, 0.125,
384 : & 0.750, 0.250, 0.750, 0.375,
385 : & 0.750, 0.500, 0.750, 0.625/
386 : DATA (X5(I),Y5(I), I = 61,81)/
387 : & 0.750, 0.750, 0.750, 0.875,
388 : & 0.750, 1.000, 0.875, 0.000,
389 : & 0.875, 0.125, 0.875, 0.250,
390 : & 0.875, 0.375, 0.875, 0.500,
391 : & 0.875, 0.625, 0.875, 0.750,
392 : & 0.875, 0.875, 0.875, 1.000,
393 : & 1.000, 0.000, 1.000, 0.125,
394 : & 1.000, 0.250, 1.000, 0.375,
395 : & 1.000, 0.500, 1.000, 0.625,
396 : & 1.000, 0.750, 1.000, 0.875,
397 : & 1.000, 1.000/
398 :
399 : ! Store node set K in (X,Y).
400 :
401 0 : if (K == 1) then
402 0 : do I = 1,100
403 0 : X(I) = X1(I)
404 0 : Y(I) = Y1(I)
405 : end do
406 0 : N = 100
407 0 : else if (K == 2) then
408 0 : do I = 1,33
409 0 : X(I) = X2(I)
410 0 : Y(I) = Y2(I)
411 : end do
412 0 : N = 33
413 0 : else if (K == 3) then
414 0 : do I = 1,25
415 0 : X(I) = X3(I)
416 0 : Y(I) = Y3(I)
417 : end do
418 0 : N = 25
419 0 : else if (K == 4) then
420 0 : do I = 1,100
421 0 : X(I) = X4(I)
422 0 : Y(I) = Y4(I)
423 : end do
424 0 : N = 100
425 0 : else if (K == 5) then
426 0 : do I = 1,81
427 0 : X(I) = X5(I)
428 0 : Y(I) = Y5(I)
429 : end do
430 0 : N = 81
431 : else
432 0 : N = 0
433 : end if
434 0 : return
435 : end subroutine testdt_db
436 :
437 0 : subroutine TSTFN1_db (K,X,Y,IFLAG, F,FX,FY)
438 : integer :: K, IFLAG
439 : double precision :: X, Y, F, FX, FY
440 :
441 : ! **********************************************************
442 :
443 : ! Robert J. Renka
444 : ! Dept. of Computer Science
445 : ! Univ. of North Texas
446 : ! renka@cs.unt.edu
447 : ! 10/14/98
448 :
449 : ! This subroutine computes the value and, optionally, the
450 : ! first partial derivatives of one of ten bivariate test
451 : ! functions. The first six functions were chosen by Richard
452 : ! Franke to test interpolation software (See the reference
453 : ! below). The last four functions represent more chal-
454 : ! lenging surface fitting problems.
455 :
456 : ! On input:
457 :
458 : ! K = integer in the range 1 to 10 which determines
459 : ! the choice of function as follows:
460 :
461 : ! K = 1 - Exponential
462 : ! K = 2 - Cliff
463 : ! K = 3 - Saddle
464 : ! K = 4 - Gentle
465 : ! K = 5 - Steep
466 : ! K = 6 - Sphere
467 : ! K = 7 - Trig
468 : ! K = 8 - Synergistic Gaussian
469 : ! K = 9 - Cloverleaf Asymmetric Peak/Valley
470 : ! K = 10 - Cosine Peak
471 :
472 : ! Note that function 6 is only defined inside a circle of
473 : ! radius 8/9 centered at (.5,.5). Thus, if (X-.5)**2 +
474 : ! (Y-.5)**2 >= 64/81, the value (and partials if IFLAG=1)
475 : ! are set to 0 for this function. Also, the first partial
476 : ! derivatives of function 10 are not defined at (.5,.5) --
477 : ! again, zeros are returned.
478 :
479 : ! X,Y = Coordinates of the point at which the selected
480 : ! function is to be evaluated.
481 :
482 : ! IFLAG = Derivative option indicator:
483 : ! IFLAG = 0 if only a function value is
484 : ! required.
485 : ! IFLAG = 1 if both the function and its first
486 : ! partial derivatives are to be
487 : ! evaluated.
488 :
489 : ! Input parameters are not altered by this routine.
490 :
491 : ! On output:
492 :
493 : ! F = Value of function K at (X,Y).
494 :
495 : ! FX,FY = First partial derivatives of function K at
496 : ! (X,Y) if IFLAG = 1, unaltered otherwise.
497 :
498 : ! Intrinsic functions called by TSTFN1: COS, EXP, SIN,
499 : ! SQRT, TANH
500 :
501 : ! Reference: R. Franke, A Critical Comparison of Some
502 : ! Methods for Interpolation of Scattered Data,
503 : ! Naval Postgraduate School Technical Report,
504 : ! NPS-53-79-003, 1979.
505 :
506 : ! **********************************************************
507 :
508 0 : double precision :: T1, T2, T3, T4
509 0 : if (K < 1 .or. K > 10) return
510 0 : GOTO (1,2,3,4,5,6,7,8,9,10), K
511 :
512 : ! Exponential:
513 :
514 : 1 F = .75*exp(-((9.*X-2.)**2 + (9.*Y-2.)**2)/4.) + .75*exp(-((9.*X+1.)**2)/49. - (9.*Y+1.)/10.)
515 0 : & + .5*exp(-((9.*X-7.)**2 + (9.*Y-3.)**2)/4.) - .2*exp(-(9.*X-4.)**2 - (9.*Y-7.)**2)
516 0 : if (IFLAG /= 1) return
517 0 : T1 = exp(-((9.*X-2.)**2 + (9.*Y-2.)**2)/4.)
518 0 : T2 = exp(-((9.*X+1.)**2)/49. - (9.*Y+1.)/10.)
519 0 : T3 = exp(-((9.*X-7.)**2 + (9.*Y-3.)**2)/4.)
520 0 : T4 = exp(-(9.*X-4.)**2 - (9.*Y-7.)**2)
521 0 : FX = -3.375*(9.*X-2.)*T1 - (27./98.)*(9.*X+1.)*T2 -2.25*(9.*X-7.)*T3 + 3.6*(9.*X-4.)*T4
522 0 : FY = -3.375*(9.*Y-2.)*T1 - .675*T2 -2.25*(9.*Y-3.)*T3 + 3.6*(9.*Y-7.)*T4
523 0 : return
524 :
525 : ! Cliff:
526 :
527 0 : 2 F = (tanh(9.0*(Y-X)) + 1.0)/9.0
528 0 : if (IFLAG /= 1) return
529 0 : T1 = 18.0*(Y-X)
530 0 : FX = -4.0/(exp(T1) + 2.0 + exp(-T1))
531 0 : FY = -FX
532 0 : return
533 :
534 : ! Saddle:
535 :
536 0 : 3 F = (1.25 + cos(5.4*Y))/(6.0 + 6.0*(3.0*X-1.0)**2)
537 0 : if (IFLAG /= 1) return
538 0 : T1 = 5.4*Y
539 0 : T2 = 1.0 + (3.0*X-1.)**2
540 0 : FX = -(3.0*X-1.0)*(1.25 + cos(T1))/(T2**2)
541 0 : FY = -.9*sin(T1)/T2
542 0 : return
543 :
544 : ! Gentle:
545 :
546 0 : 4 F = exp(-5.0625*((X-.5)**2 + (Y-.5)**2))/3.0
547 0 : if (IFLAG /= 1) return
548 0 : T1 = X - .5
549 0 : T2 = Y - .5
550 0 : T3 = -3.375*exp(-5.0625*(T1**2 + T2**2))
551 0 : FX = T1*T3
552 0 : FY = T2*T3
553 0 : return
554 :
555 : ! Steep:
556 :
557 0 : 5 F = exp(-20.25*((X-.5)**2 + (Y-.5)**2))/3.0
558 0 : if (IFLAG /= 1) return
559 0 : T1 = X - .5
560 0 : T2 = Y - .5
561 0 : T3 = -13.5*exp(-20.25*(T1**2 + T2**2))
562 0 : FX = T1*T3
563 0 : FY = T2*T3
564 0 : return
565 :
566 : ! Sphere:
567 :
568 0 : 6 T4 = 64.0 - 81.0*((X-.5)**2 + (Y-.5)**2)
569 0 : F = 0.
570 0 : if (T4 >= 0.) F = sqrt(T4)/9.0 - .5
571 0 : if (IFLAG /= 1) return
572 0 : T1 = X - .5
573 0 : T2 = Y - .5
574 0 : T3 = 0.
575 0 : if (T4 > 0.) T3 = -9.0/sqrt(T4)
576 0 : FX = T1*T3
577 0 : FY = T2*T3
578 0 : return
579 :
580 : ! Trig:
581 :
582 0 : 7 F = 2.0*cos(10.0*X)*sin(10.0*Y) + sin(10.0*X*Y)
583 0 : if (IFLAG /= 1) return
584 0 : T1 = 10.0*X
585 0 : T2 = 10.0*Y
586 0 : T3 = 10.0*cos(10.0*X*Y)
587 0 : FX = -20.0*sin(T1)*sin(T2) + T3*Y
588 0 : FY = 20.0*cos(T1)*cos(T2) + T3*X
589 0 : return
590 :
591 : ! Gaussx(1,.5,.1) + Gaussy(.75,.5,.1) + Gaussx(1,.5,.1)*
592 : ! Gaussy(.75,.5,.1), where Gaussx(a,b,c) is the Gaussian
593 : ! function of x with amplitude a, center (mean) b, and
594 : ! width (standard deviation) c.
595 :
596 0 : 8 T1 = 5.0 - 10.0*X
597 0 : T2 = 5.0 - 10.0*Y
598 0 : T3 = exp(-.5*T1*T1)
599 0 : T4 = exp(-.5*T2*T2)
600 0 : F = T3 + .75*T4*(1.0+T3)
601 0 : if (IFLAG /= 1) return
602 0 : FX = T1*T3*(10.0 + 7.5*T4)
603 0 : FY = T2*T4*(7.5 + 7.5*T3)
604 0 : return
605 :
606 : ! Cloverleaf Asymmetric Hill/Valley:
607 :
608 0 : 9 T1 = exp((10.0 - 20.0*X)/3.0)
609 0 : T2 = exp((10.0 - 20.0*Y)/3.0)
610 0 : T3 = 1.0/(1.0 + T1)
611 0 : T4 = 1.0/(1.0 + T2)
612 0 : F = ((20.0/3.0)**3 * T1*T2)**2 * (T3*T4)**5 * (T1-2.0*T3)*(T2-2.0*T4)
613 0 : if (IFLAG /= 1) return
614 0 : FX = ((20.0/3.0)*T1)**2 * ((20.0/3.0)*T3)**5 * (2.0*T1-3.0*T3-5.0+12.0*T3*T3)*T2*T2*T4**5 * (T2-2.0*T4)
615 0 : FY = ((20.0/3.0)*T1)**2 * ((20.0/3.0)*T3)**5 * (2.0*T2-3.0*T4-5.0+12.0*T4*T4)*T2*T2*T4**5 * (T1-2.0*T3)
616 0 : return
617 :
618 : ! Cosine Peak:
619 :
620 0 : 10 T1 = sqrt( (80.0*X - 40.0)**2 + (90.0*Y - 45.0)**2 )
621 0 : T2 = exp(-.04*T1)
622 0 : T3 = cos(.15*T1)
623 0 : F = T2*T3
624 0 : if (IFLAG /= 1) return
625 0 : T4 = sin(.15*T1)
626 0 : FX = 0.
627 0 : FY = 0.
628 0 : if (T1 == 0.) return
629 0 : T4 = sin(.15*T1)
630 0 : FX = -T2*(12.0*T4 + 3.2*T3)*(80.0*X - 40.0)/T1
631 0 : FY = -T2*(13.5*T4 + 3.6*T3)*(90.0*Y - 45.0)/T1
632 0 : return
633 : end subroutine TSTFN1_db
634 :
635 :
636 0 : subroutine TSTFN2_db (K,X,Y,IFLAG, F,FX,FY,FXX,FXY,FYY)
637 : integer :: K, IFLAG
638 : double precision :: X, Y, F, FX, FY, FXX, FXY, FYY
639 :
640 : ! **********************************************************
641 :
642 : ! Robert J. Renka
643 : ! Dept. of Computer Science
644 : ! Univ. of North Texas
645 : ! renka@cs.unt.edu
646 : ! 10/14/98
647 :
648 : ! This subroutine computes the value and, optionally, the
649 : ! first and/or second partial derivatives of one of ten
650 : ! bivariate test functions. The first six functions were
651 : ! chosen by Richard Franke to test interpolation software.
652 : ! (See the reference below.) The last four functions repre-
653 : ! sent more challenging surface fitting problems.
654 :
655 : ! On input:
656 :
657 : ! K = integer in the range 1 to 10 which determines
658 : ! the choice of function as follows:
659 :
660 : ! K = 1 - Exponential
661 : ! K = 2 - Cliff
662 : ! K = 3 - Saddle
663 : ! K = 4 - Gentle
664 : ! K = 5 - Steep
665 : ! K = 6 - Sphere
666 : ! K = 7 - Trig
667 : ! K = 8 - Synergistic Gaussian
668 : ! K = 9 - Cloverleaf Asymmetric Peak/Valley
669 : ! K = 10 - Cosine Peak
670 :
671 : ! Note that function 6 is only defined inside a circle of
672 : ! radius 8/9 centered at (.5,.5). Thus, if (X-.5)**2 +
673 : ! (Y-.5)**2 >= 64/81, the value (and partials if IFLAG=1)
674 : ! are set to 0 for this function. Also, the first partial
675 : ! derivatives of function 10 are not defined at (.5,.5) --
676 : ! again, zeros are returned.
677 :
678 : ! X,Y = Coordinates of the point at which the selected
679 : ! function is to be evaluated.
680 :
681 : ! IFLAG = Derivative option indicator:
682 : ! IFLAG = 0 if only a function value is
683 : ! required.
684 : ! IFLAG = 1 if both the function and its first
685 : ! partial derivatives are to be
686 : ! evaluated.
687 : ! IFLAG = 2 if the function, its first partial
688 : ! derivatives, and its second partial
689 : ! derivatives are to be evaluated.
690 :
691 : ! Input parameters are not altered by this routine.
692 :
693 : ! On output:
694 :
695 : ! F = Value of function K at (X,Y).
696 :
697 : ! FX,FY = First partial derivatives of function K at
698 : ! (X,Y) if IFLAG >= 1, unaltered otherwise.
699 :
700 : ! FXX,FXY,FYY = Second partial derivatives of function
701 : ! K at (X,Y) if IFLAG >= 2, unaltered
702 : ! otherwise.
703 :
704 : ! Intrinsic functions called by TSTFN2: COS, EXP, SIN,
705 : ! SQRT, TANH
706 :
707 : ! Reference: R. Franke, A Critical Comparison of Some
708 : ! Methods for Interpolation of Scattered Data,
709 : ! Naval Postgraduate School Technical Report,
710 : ! NPS-53-79-003, 1979.
711 :
712 : ! **********************************************************
713 :
714 0 : double precision :: T1, T2, T3, T4, T5, T6
715 0 : if (K < 1 .or. K > 10) return
716 0 : GOTO (1,2,3,4,5,6,7,8,9,10), K
717 :
718 : ! Exponential:
719 :
720 : 1 F = .75*exp(-((9.*X-2.)**2 + (9.*Y-2.)**2)/4.) + .75*exp(-((9.*X+1.)**2)/49. - (9.*Y+1.)/10.)
721 0 : & + .5*exp(-((9.*X-7.)**2 + (9.*Y-3.)**2)/4.) - .2*exp(-(9.*X-4.)**2 - (9.*Y-7.)**2)
722 0 : if (IFLAG < 1) return
723 0 : T1 = exp(-((9.*X-2.)**2 + (9.*Y-2.)**2)/4.)
724 0 : T2 = exp(-((9.*X+1.)**2)/49. - (9.*Y+1.)/10.)
725 0 : T3 = exp(-((9.*X-7.)**2 + (9.*Y-3.)**2)/4.)
726 0 : T4 = exp(-(9.*X-4.)**2 - (9.*Y-7.)**2)
727 0 : FX = -3.375*(9.*X-2.)*T1 - (27./98.)*(9.*X+1.)*T2 -2.25*(9.*X-7.)*T3 + 3.6*(9.*X-4.)*T4
728 0 : FY = -3.375*(9.*Y-2.)*T1 - .675*T2 -2.25*(9.*Y-3.)*T3 + 3.6*(9.*Y-7.)*T4
729 0 : if (IFLAG < 2) return
730 0 : FXX = 15.1875*((9.*X-2.)**2-2.)*T1 + 60.75*((9.*X+1.)**2-24.5)*T2 + 10.125*((9.*X-7.)**2-2.)*T3 - 64.8*((9.*X-4.)**2-.5)*T4
731 0 : FXY = 15.1875*(9.*X-2.)*(9.*Y-2.)*T1+(243./980.)*(9.*X+1.)*T2 + 10.125*(9.*X-7.)*(9.*Y-3.)*T3 - 64.8*(9.*X-4.)*(9.*Y-7.)*T4
732 0 : FYY = 15.1875*((9.*Y-2.)**2 - 2.)*T1 + .6075*T2 + 10.125*((9.*Y-3.)**2 - 2.)*T3 - 64.8*((9.*Y-7.)**2 - .5)*T4
733 0 : return
734 :
735 : ! Cliff:
736 :
737 0 : 2 F = (tanh(9.0*(Y-X)) + 1.0)/9.0
738 0 : if (IFLAG < 1) return
739 0 : T1 = 18.0*(Y-X)
740 0 : FX = -4.0/(exp(T1) + 2.0 + exp(-T1))
741 0 : FY = -FX
742 0 : if (IFLAG < 2) return
743 0 : FXX = 18.0*tanh(0.5*T1)*FX
744 0 : FXY = -FXX
745 0 : FYY = FXX
746 0 : return
747 :
748 : ! Saddle:
749 :
750 0 : 3 F = (1.25 + cos(5.4*Y))/(6.0 + 6.0*(3.0*X-1.0)**2)
751 0 : if (IFLAG < 1) return
752 0 : T1 = 5.4*Y
753 0 : T2 = 1.0 + (3.0*X-1.)**2
754 0 : FX = -(3.0*X-1.0)*(1.25 + cos(T1))/(T2**2)
755 0 : FY = -.9*sin(T1)/T2
756 0 : if (IFLAG < 2) return
757 0 : FXX = 3.0*(1.25 + cos(T1))*(3.0*T2-4.0)/(T2**3)
758 0 : FXY = 5.4*(3.0*X-1.0)*sin(T1)/(T2**2)
759 0 : FYY = -4.86*cos(T1)/T2
760 0 : return
761 :
762 : ! Gentle:
763 :
764 0 : 4 F = exp(-5.0625*((X-.5)**2 + (Y-.5)**2))/3.0
765 0 : if (IFLAG < 1) return
766 0 : T1 = X - .5
767 0 : T2 = Y - .5
768 0 : T3 = -3.375*exp(-5.0625*(T1**2 + T2**2))
769 0 : FX = T1*T3
770 0 : FY = T2*T3
771 0 : if (IFLAG < 2) return
772 0 : FXX = (1.0 - 10.125*T1*T1)*T3
773 0 : FXY = -10.125*T1*T2*T3
774 0 : FYY = (1.0 - 10.125*T2*T2)*T3
775 0 : return
776 :
777 : ! Steep:
778 :
779 0 : 5 F = exp(-20.25*((X-.5)**2 + (Y-.5)**2))/3.0
780 0 : if (IFLAG < 1) return
781 0 : T1 = X - .5
782 0 : T2 = Y - .5
783 0 : T3 = -13.5*exp(-20.25*(T1**2 + T2**2))
784 0 : FX = T1*T3
785 0 : FY = T2*T3
786 0 : if (IFLAG < 2) return
787 0 : FXX = (1.0 - 40.5*T1*T1)*T3
788 0 : FXY = -40.5*T1*T2*T3
789 0 : FYY = (1.0 - 40.5*T2*T2)*T3
790 0 : return
791 :
792 : ! Sphere:
793 :
794 0 : 6 T4 = 64.0 - 81.0*((X-.5)**2 + (Y-.5)**2)
795 0 : F = 0.
796 0 : if (T4 >= 0.) F = sqrt(T4)/9.0 - .5
797 0 : if (IFLAG < 1) return
798 0 : T1 = X - .5
799 0 : T2 = Y - .5
800 0 : T3 = 0.
801 0 : if (T4 > 0.) T3 = -9.0/sqrt(T4)
802 0 : FX = T1*T3
803 0 : FY = T2*T3
804 0 : if (IFLAG < 2) return
805 0 : FXX = (1.0 + FX*FX)*T3
806 0 : FXY = FX*FY
807 0 : FYY = (1.0 + FY*FY)*T3
808 0 : return
809 :
810 : ! Trig:
811 :
812 0 : 7 F = 2.0*cos(10.0*X)*sin(10.0*Y) + sin(10.0*X*Y)
813 0 : if (IFLAG < 1) return
814 0 : T1 = 10.0*X
815 0 : T2 = 10.0*Y
816 0 : T3 = 10.0*cos(10.0*X*Y)
817 0 : FX = -20.0*sin(T1)*sin(T2) + T3*Y
818 0 : FY = 20.0*cos(T1)*cos(T2) + T3*X
819 0 : if (IFLAG < 2) return
820 0 : T4 = 100.0*sin(10.0*X*Y)
821 0 : FXX = -200.0*cos(T1)*sin(T2) - T4*Y*Y
822 0 : FXY = -200.0*sin(T1)*cos(T2) + T3 - T4*X*Y
823 0 : FYY = -200.0*cos(T1)*sin(T2) - T4*X*X
824 0 : return
825 :
826 : ! Gaussx(1,.5,.1) + Gaussy(.75,.5,.1) + Gaussx(1,.5,.1)*
827 : ! Gaussy(.75,.5,.1), where Gaussx(a,b,c) is the Gaussian
828 : ! function of x with amplitude a, center (mean) b, and
829 : ! width (standard deviation) c.
830 :
831 0 : 8 T1 = 5.0 - 10.0*X
832 0 : T2 = 5.0 - 10.0*Y
833 0 : T3 = exp(-.5*T1*T1)
834 0 : T4 = exp(-.5*T2*T2)
835 0 : F = T3 + .75*T4*(1.0+T3)
836 0 : if (IFLAG < 1) return
837 0 : FX = T1*T3*(10.0 + 7.5*T4)
838 0 : FY = T2*T4*(7.5 + 7.5*T3)
839 0 : if (IFLAG < 2) return
840 0 : FXX = T3*(T1*T1-1.0)*(100.0 + 75.0*T4)
841 0 : FXY = 75.0*T1*T2*T3*T4
842 0 : FYY = T4*(T2*T2-1.0)*(75.0 + 75.0*T3)
843 0 : return
844 :
845 : ! Cloverleaf Asymmetric Hill/Valley:
846 :
847 0 : 9 T1 = exp((10.0 - 20.0*X)/3.0)
848 0 : T2 = exp((10.0 - 20.0*Y)/3.0)
849 0 : T3 = 1.0/(1.0 + T1)
850 0 : T4 = 1.0/(1.0 + T2)
851 0 : T5 = 20.0/3.0
852 0 : F = (T5**3 * T1*T2)**2 * (T3*T4)**5 * (T1-2.0*T3)*(T2-2.0*T4)
853 0 : if (IFLAG < 1) return
854 0 : T6 = (T5*T1*T2)**2 * (T5*T3*T4)**5
855 0 : FX = T6 * (T2-2.0*T4) * ((12.0*T3 - 3.0)*T3 + 2.0*T1 - 5.0)
856 0 : FY = T6 * (T1-2.0*T3) * ((12.0*T4 - 3.0)*T4 + 2.0*T2 - 5.0)
857 0 : if (IFLAG < 2) return
858 0 : FXX = T5*T6 * (T2-2.0*T4) * (((-84.0*T3 + 78.0)*T3 + 23.0)*T3 + 4.0*T1-25.0)
859 0 : FXY = T5*T6 * ((12.0*T4 - 3.0)*T4 + 2.0*T2 - 5.0) * ((12.0*T3 - 3.0)*T3 + 2.0*T1 - 5.0)
860 0 : FYY = T5*T6 * (T1-2.0*T3) * (((-84.0*T4 + 78.0)*T4 + 23.0)*T4 + 4.0*T2-25.0)
861 0 : return
862 :
863 : ! Cosine Peak:
864 :
865 0 : 10 T1 = sqrt( (80.0*X - 40.0)**2 + (90.0*Y - 45.0)**2 )
866 0 : T2 = exp(-.04*T1)
867 0 : T3 = cos(.15*T1)
868 0 : F = T2*T3
869 0 : if (IFLAG < 1) return
870 0 : T4 = sin(.15*T1)
871 0 : FX = 0.
872 0 : FY = 0.
873 0 : if (T1 == 0.) return
874 0 : T4 = sin(.15*T1)
875 0 : FX = -T2*(12.0*T4 + 3.2*T3)*(80.0*X - 40.0)/T1
876 0 : FY = -T2*(13.5*T4 + 3.6*T3)*(90.0*Y - 45.0)/T1
877 0 : if (IFLAG < 2) return
878 0 : FXX = 0.
879 0 : FXY = 0.
880 0 : FYY = 0.
881 0 : if (T1 == 0.) return
882 0 : T5 = T2/(T1**3)
883 0 : FXX = T5*(T1*(76.8*T4 - 133.76*T3)*(80.0*X-40.0)**2 - (960.0*T4 + 256.0*T3)*(90.0*Y-45.0)**2 )
884 0 : FXY = T5*(T1*(86.4*T4 - 150.48*T3) + 1080.0*T4 + 288.0*T3)*(80.0*X-40.0)*(90.0*Y-45.0)
885 0 : FYY = T5*(T1*(97.2*T4 - 169.29*T3)*(90.0*Y-45.0)**2 - (1215.0*T4 + 324.0*T3)*(80.0*X-40.0)**2)
886 0 : return
887 : end subroutine TSTFN2_db
|
