1 /////////////////////////////////////////////////////////////////////////////
2 // einspline: a library for creating and evaluating B-splines //
3 // Copyright (C) 2007 Kenneth P. Esler, Jr. //
4 // Released under the BSD-3-clause license //
5 /////////////////////////////////////////////////////////////////////////////
6
7 #include "multi_bspline.h"
8 #include "bspline.h"
9 #include <stdio.h>
10 #include <stdlib.h>
11 #include <math.h>
12 #include <time.h>
13
14
diff(double a,double b,double tol)15 inline double diff (double a, double b, double tol)
16 {
17 if (fabs(a-b) > tol)
18 return 1;
19 else
20 return 0;
21 }
22
23
24 //////////////////////////////////////////
25 // Single-precision real test functions //
26 //////////////////////////////////////////
27 int
test_1d_float_all()28 test_1d_float_all()
29 {
30 int Nx=73;
31 int num_splines = 21;
32
33 Ugrid x_grid;
34 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
35
36 BCtype_s xBC;
37 xBC.lCode = xBC.rCode = PERIODIC;
38
39 // First, create splines the normal way
40 UBspline_1d_s* norm_splines[num_splines];
41 multi_UBspline_1d_s *multi_spline;
42
43 // First, create multispline
44 multi_spline = create_multi_UBspline_1d_s (x_grid, xBC, num_splines);
45
46 float data[Nx];
47 // Now, create normal splines and set multispline data
48 for (int i=0; i<num_splines; i++) {
49 for (int j=0; j<Nx; j++)
50 data[j] = (drand48()-0.5);
51 norm_splines[i] = create_UBspline_1d_s (x_grid, xBC, data);
52 set_multi_UBspline_1d_s (multi_spline, i, data);
53 }
54
55 // fprintf (stderr, "\nnorm coef = %1.14e\n",
56 // norm_splines[19]->coefs[27]);
57 // fprintf (stderr, "multi coef = %1.14e\n",
58 // multi_spline->coefs[19+27*multi_spline->x_stride]);
59
60 // Now, test random values
61 int num_vals = 100;
62 float multi_vals[num_splines], norm_vals [num_splines];
63 float multi_grads[num_splines], norm_grads[num_splines];
64 float multi_lapl[num_splines], norm_lapl [num_splines];
65 for (int i=0; i<num_vals; i++) {
66 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
67
68 //////////////////////////
69 // Check value routine //
70 //////////////////////////
71 eval_multi_UBspline_1d_s (multi_spline, x, multi_vals);
72 for (int j=0; j<num_splines; j++)
73 eval_UBspline_1d_s (norm_splines[j], x, &(norm_vals[j]));
74 for (int j=0; j<num_splines; j++) {
75 // Check value
76 if (diff(norm_vals[j], multi_vals[j], 1.0e-6)) {
77 fprintf (stderr, " norm_vals[j] = %1.8e\n", norm_vals[j]);
78 fprintf (stderr, "multi_vals[j] = %1.8e\n", multi_vals[j]);
79 return -1;
80 }
81 }
82
83 ///////////////////////
84 // Check VG routine //
85 ///////////////////////
86 eval_multi_UBspline_1d_s_vg (multi_spline, x,
87 multi_vals, multi_grads);
88 for (int j=0; j<num_splines; j++)
89 eval_UBspline_1d_s_vg (norm_splines[j], x, &(norm_vals[j]),
90 &(norm_grads[j]));
91 for (int j=0; j<num_splines; j++) {
92 // Check value
93 if (diff(norm_vals[j], multi_vals[j], 1.0e-6))
94 return -2;
95
96 // Check gradients
97 if (diff (norm_grads[j], multi_grads[j], 1.0e-5))
98 return -3;
99 }
100
101
102 ///////////////////////
103 // Check VGL routine //
104 ///////////////////////
105 eval_multi_UBspline_1d_s_vgl (multi_spline, x, multi_vals, multi_grads, multi_lapl);
106 for (int j=0; j<num_splines; j++)
107 eval_UBspline_1d_s_vgl (norm_splines[j], x, &(norm_vals[j]),
108 &(norm_grads[j]), &(norm_lapl[j]));
109 for (int j=0; j<num_splines; j++) {
110 // Check value
111 if (diff(norm_vals[j], multi_vals[j], 1.0e-6))
112 return -4;
113
114 // Check gradients
115 if (diff (norm_grads[j], multi_grads[j], 1.0e-5))
116 return -5;
117
118 // Check laplacian
119 if (diff (norm_lapl[j], multi_lapl[j], 1.0e-3))
120 return -6;
121 }
122 }
123 return 0;
124 }
125
126
127
128 int
test_2d_float_all()129 test_2d_float_all()
130 {
131 int Nx=73; int Ny=91;
132 int num_splines = 21;
133
134 Ugrid x_grid, y_grid;
135 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
136 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
137
138 BCtype_s xBC, yBC;
139 xBC.lCode = xBC.rCode = PERIODIC;
140 yBC.lCode = yBC.rCode = PERIODIC;
141
142 // First, create splines the normal way
143 UBspline_2d_s* norm_splines[num_splines];
144 multi_UBspline_2d_s *multi_spline;
145
146 // First, create multispline
147 multi_spline = create_multi_UBspline_2d_s (x_grid, y_grid, xBC, yBC,
148 num_splines);
149
150 float data[Nx*Ny];
151 // Now, create normal splines and set multispline data
152 for (int i=0; i<num_splines; i++) {
153 for (int j=0; j<Nx*Ny; j++)
154 data[j] = (drand48()-0.5);
155 norm_splines[i] = create_UBspline_2d_s (x_grid, y_grid, xBC, yBC, data);
156 set_multi_UBspline_2d_s (multi_spline, i, data);
157 }
158
159 // fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
160 // real(norm_splines[19]->coefs[227]),
161 // imag(norm_splines[19]->coefs[227]));
162 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
163 // real(multi_spline->coefs[19+227*multi_spline->z_stride]),
164 // imag(multi_spline->coefs[19+227*multi_spline->z_stride]));
165
166 // Now, test random values
167 int num_vals = 100;
168 float multi_vals[num_splines], norm_vals[num_splines];
169 float multi_grads[2*num_splines], norm_grads[2*num_splines];
170 float multi_lapl[num_splines], norm_lapl[num_splines];
171 float multi_hess[4*num_splines], norm_hess[4*num_splines];
172 for (int i=0; i<num_vals; i++) {
173 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
174 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
175
176
177 //////////////////////////
178 // Check value routine //
179 //////////////////////////
180 eval_multi_UBspline_2d_s (multi_spline, x, y, multi_vals);
181 for (int j=0; j<num_splines; j++)
182 eval_UBspline_2d_s (norm_splines[j], x, y, &(norm_vals[j]));
183 for (int j=0; j<num_splines; j++) {
184 // Check value
185 if (diff(norm_vals[j], multi_vals[j], 1.0e-6))
186 return -1;
187 }
188
189 ///////////////////////
190 // Check VG routine //
191 ///////////////////////
192 eval_multi_UBspline_2d_s_vg (multi_spline, x, y,
193 multi_vals, multi_grads);
194 for (int j=0; j<num_splines; j++)
195 eval_UBspline_2d_s_vg (norm_splines[j], x, y, &(norm_vals[j]),
196 &(norm_grads[2*j]));
197 for (int j=0; j<num_splines; j++) {
198 // Check value
199 if (diff(norm_vals[j], multi_vals[j], 1.0e-6))
200 return -1;
201
202 // Check gradients
203 for (int n=0; n<2; n++)
204 if (diff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-5))
205 return -2;
206 }
207
208
209 ///////////////////////
210 // Check VGL routine //
211 ///////////////////////
212 eval_multi_UBspline_2d_s_vgl (multi_spline, x, y,
213 multi_vals, multi_grads, multi_lapl);
214 for (int j=0; j<num_splines; j++)
215 eval_UBspline_2d_s_vgl (norm_splines[j], x, y, &(norm_vals[j]),
216 &(norm_grads[2*j]), &(norm_lapl[j]));
217 for (int j=0; j<num_splines; j++) {
218 // Check value
219 if (diff(norm_vals[j], multi_vals[j], 1.0e-6))
220 return -3;
221
222 // Check gradients
223 for (int n=0; n<2; n++)
224 if (diff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-5))
225 return -4;
226
227 // Check laplacian
228 if (diff (norm_lapl[j], multi_lapl[j], 1.0e-3))
229 return -5;
230 }
231
232
233 ///////////////////////
234 // Check VGH routine //
235 ///////////////////////
236 eval_multi_UBspline_2d_s_vgh (multi_spline, x, y,
237 multi_vals, multi_grads, multi_hess);
238 for (int j=0; j<num_splines; j++)
239 eval_UBspline_2d_s_vgh (norm_splines[j], x, y, &(norm_vals[j]),
240 &(norm_grads[2*j]), &(norm_hess[4*j]));
241 for (int j=0; j<num_splines; j++) {
242 // Check value
243 if (diff(norm_vals[j], multi_vals[j], 1.0e-6)) {
244 fprintf (stderr, "j = %d\n", j);
245 fprintf (stderr, "norm_vals[j] = %1.14e\n", norm_vals[j]);
246 fprintf (stderr, "multi_vals[j] = %1.14e\n", multi_vals[j]);
247 //return -6;
248 }
249
250 // Check gradients
251 for (int n=0; n<2; n++)
252 if (diff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-5))
253 return -7;
254
255 // Check hessian
256 for (int n=0; n<4; n++)
257 if (diff (norm_hess[4*j+n], multi_hess[4*j+n], 1.0e-3)) {
258 fprintf (stderr, "j = %d n = %d \n", j, n);
259 fprintf (stderr, "norm_hess[j] = %1.14e\n", norm_hess[4*j+n]);
260 fprintf (stderr, "multi_hess[j] = %1.14e\n", multi_hess[4*j+n]);
261 //return -8;
262 }
263 }
264 }
265 return 0;
266 }
267
268
269 int
test_3d_float_all()270 test_3d_float_all()
271 {
272 int Nx=73; int Ny=91; int Nz = 29;
273 int num_splines = 23;
274
275 Ugrid x_grid, y_grid, z_grid;
276 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
277 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
278 z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
279
280 BCtype_s xBC, yBC, zBC;
281 xBC.lCode = xBC.rCode = PERIODIC;
282 yBC.lCode = yBC.rCode = PERIODIC;
283 zBC.lCode = zBC.rCode = PERIODIC;
284
285 // First, create splines the normal way
286 UBspline_3d_s* norm_splines[num_splines];
287 multi_UBspline_3d_s *multi_spline;
288
289 // First, create multispline
290 multi_spline = create_multi_UBspline_3d_s (x_grid, y_grid, z_grid, xBC, yBC, zBC,
291 num_splines);
292
293 float data[Nx*Ny*Nz];
294 // Now, create normal splines and set multispline data
295 for (int i=0; i<num_splines; i++) {
296 for (int j=0; j<Nx*Ny*Nz; j++)
297 data[j] = (drand48()-0.5);
298 norm_splines[i] = create_UBspline_3d_s (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
299 set_multi_UBspline_3d_s (multi_spline, i, data);
300 }
301
302 // fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
303 // real(norm_splines[19]->coefs[227]),
304 // imag(norm_splines[19]->coefs[227]));
305 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
306 // real(multi_spline->coefs[19+227*multi_spline->z_stride]),
307 // imag(multi_spline->coefs[19+227*multi_spline->z_stride]));
308
309 // Now, test random values
310 int num_vals = 100;
311 float multi_vals[num_splines], norm_vals[num_splines];
312 float multi_grads[3*num_splines], norm_grads[3*num_splines];
313 float multi_lapl[num_splines], norm_lapl[num_splines];
314 float multi_hess[9*num_splines], norm_hess[9*num_splines];
315 for (int i=0; i<num_vals; i++) {
316 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
317 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
318 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
319
320 //////////////////////////
321 // Check value routine //
322 /////////////////////////
323 eval_multi_UBspline_3d_s (multi_spline, x, y, z, multi_vals);
324 for (int j=0; j<num_splines; j++)
325 eval_UBspline_3d_s (norm_splines[j], x, y, z, &(norm_vals[j]));
326 for (int j=0; j<num_splines; j++) {
327 // Check value
328 if (diff(norm_vals[j], multi_vals[j], 1.0e-6))
329 return -1;
330 }
331
332 ///////////////////////
333 // Check VG routine //
334 ///////////////////////
335 eval_multi_UBspline_3d_s_vg (multi_spline, x, y, z,
336 multi_vals, multi_grads);
337 for (int j=0; j<num_splines; j++)
338 eval_UBspline_3d_s_vg (norm_splines[j], x, y, z, &(norm_vals[j]),
339 &(norm_grads[3*j]));
340 for (int j=0; j<num_splines; j++) {
341 // Check value
342 if (diff(norm_vals[j], multi_vals[j], 1.0e-6))
343 return -1;
344
345 // Check gradients
346 for (int n=0; n<3; n++)
347 if (diff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-4))
348 return -2;
349 }
350
351
352 ///////////////////////
353 // Check VGL routine //
354 ///////////////////////
355 eval_multi_UBspline_3d_s_vgl (multi_spline, x, y, z,
356 multi_vals, multi_grads, multi_lapl);
357 for (int j=0; j<num_splines; j++)
358 eval_UBspline_3d_s_vgl (norm_splines[j], x, y, z, &(norm_vals[j]),
359 &(norm_grads[3*j]), &(norm_lapl[j]));
360 for (int j=0; j<num_splines; j++) {
361 // Check value
362 if (diff(norm_vals[j], multi_vals[j], 1.0e-6))
363 return -3;
364
365 // Check gradients
366 for (int n=0; n<3; n++)
367 if (diff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-4))
368 return -4;
369
370 // Check laplacian
371 if (diff (norm_lapl[j], multi_lapl[j], 1.0e-3))
372 return -5;
373 }
374
375
376 ///////////////////////
377 // Check VGH routine //
378 ///////////////////////
379 eval_multi_UBspline_3d_s_vgh (multi_spline, x, y, z,
380 multi_vals, multi_grads, multi_hess);
381 for (int j=0; j<num_splines; j++)
382 eval_UBspline_3d_s_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
383 &(norm_grads[3*j]), &(norm_hess[9*j]));
384 for (int j=0; j<num_splines; j++) {
385 // Check value
386 if (diff(norm_vals[j], multi_vals[j], 1.0e-6))
387 return -6;
388
389 // Check gradients
390 for (int n=0; n<3; n++)
391 if (diff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-4)) {
392 fprintf (stderr, "n=%d j=%d\n", n, j);
393 fprintf (stderr, " norm_grads[3*j+n] = %1.8e\n",
394 norm_grads[3*j+n]);
395 fprintf (stderr, "multi_grads[3*j+n] = %1.8e\n",
396 multi_grads[3*j+n]);
397 //return -7;
398 }
399 // Check hessian
400 for (int n=0; n<9; n++)
401 if (diff (norm_hess[9*j+n], multi_hess[9*j+n], 1.0e-3))
402 return -8;
403 }
404 }
405
406
407 // num_vals = 100000;
408
409 // // Now do timing
410 // clock_t norm_start, norm_end, multi_start, multi_end, rand_start, rand_end;
411 // rand_start = clock();
412 // for (int i=0; i<num_vals; i++) {
413 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
414 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
415 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
416 // }
417 // rand_end = clock();
418
419 // norm_start = clock();
420 // for (int i=0; i<num_vals; i++) {
421 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
422 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
423 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
424
425 // for (int j=0; j<num_splines; j++)
426 // eval_UBspline_3d_s_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
427 // &(norm_grads[3*j]), &norm_hess[9*j]);
428 // }
429 // norm_end = clock();
430
431 // multi_start = clock();
432 // for (int i=0; i<num_vals; i++) {
433 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
434 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
435 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
436 // eval_multi_UBspline_3d_s_vgh (multi_spline, x, y, z, multi_vals,
437 // multi_grads, multi_hess);
438 // }
439 // multi_end = clock();
440
441 // fprintf (stderr, "Normal spline time = %1.5f\n",
442 // (double)(norm_end-norm_start+rand_start-rand_end)/CLOCKS_PER_SEC);
443 // fprintf (stderr, "Multi spline time = %1.5f\n",
444 // (double)(multi_end-multi_start+rand_start-rand_end)/CLOCKS_PER_SEC);
445
446 return 0;
447 }
448
449
450
451
452 //////////////////////////////////////////
453 // Double-precision real test functions //
454 //////////////////////////////////////////
455 int
test_1d_double_all()456 test_1d_double_all()
457 {
458 int Nx=73;
459 int num_splines = 21;
460
461 Ugrid x_grid;
462 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
463
464 BCtype_d xBC;
465 xBC.lCode = xBC.rCode = PERIODIC;
466
467 // First, create splines the normal way
468 UBspline_1d_d* norm_splines[num_splines];
469 multi_UBspline_1d_d *multi_spline;
470
471 // First, create multispline
472 multi_spline = create_multi_UBspline_1d_d (x_grid, xBC, num_splines);
473
474 double data[Nx];
475 // Now, create normal splines and set multispline data
476 for (int i=0; i<num_splines; i++) {
477 for (int j=0; j<Nx; j++)
478 data[j] = (drand48()-0.5);
479 norm_splines[i] = create_UBspline_1d_d (x_grid, xBC, data);
480 set_multi_UBspline_1d_d (multi_spline, i, data);
481 }
482
483 // Now, test random values
484 int num_vals = 100;
485 double multi_vals[num_splines], norm_vals [num_splines];
486 double multi_grads[num_splines], norm_grads[num_splines];
487 double multi_lapl[num_splines], norm_lapl [num_splines];
488 for (int i=0; i<num_vals; i++) {
489 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
490
491 //////////////////////////
492 // Check value routine //
493 //////////////////////////
494 eval_multi_UBspline_1d_d (multi_spline, x, multi_vals);
495 for (int j=0; j<num_splines; j++)
496 eval_UBspline_1d_d (norm_splines[j], x, &(norm_vals[j]));
497 for (int j=0; j<num_splines; j++) {
498 // Check value
499 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
500 return -1;
501 }
502
503 ///////////////////////
504 // Check VG routine //
505 ///////////////////////
506 eval_multi_UBspline_1d_d_vg (multi_spline, x,
507 multi_vals, multi_grads);
508 for (int j=0; j<num_splines; j++)
509 eval_UBspline_1d_d_vg (norm_splines[j], x, &(norm_vals[j]),
510 &(norm_grads[j]));
511 for (int j=0; j<num_splines; j++) {
512 // Check value
513 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
514 return -1;
515
516 // Check gradients
517 if (diff (norm_grads[j], multi_grads[j], 1.0e-12))
518 return -2;
519 }
520
521
522 ///////////////////////
523 // Check VGL routine //
524 ///////////////////////
525 eval_multi_UBspline_1d_d_vgl (multi_spline, x, multi_vals, multi_grads, multi_lapl);
526 for (int j=0; j<num_splines; j++)
527 eval_UBspline_1d_d_vgl (norm_splines[j], x, &(norm_vals[j]),
528 &(norm_grads[j]), &(norm_lapl[j]));
529 for (int j=0; j<num_splines; j++) {
530 // Check value
531 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
532 return -3;
533
534 // Check gradients
535 if (diff (norm_grads[j], multi_grads[j], 1.0e-10))
536 return -4;
537
538 // Check laplacian
539 if (diff (norm_lapl[j], multi_lapl[j], 1.0e-10))
540 return -5;
541 }
542 }
543 return 0;
544 }
545
546
547
548 int
test_2d_double_all()549 test_2d_double_all()
550 {
551 int Nx=73; int Ny=91;
552 int num_splines = 21;
553
554 Ugrid x_grid, y_grid;
555 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
556 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
557
558 BCtype_d xBC, yBC;
559 xBC.lCode = xBC.rCode = PERIODIC;
560 yBC.lCode = yBC.rCode = PERIODIC;
561
562 // First, create splines the normal way
563 UBspline_2d_d* norm_splines[num_splines];
564 multi_UBspline_2d_d *multi_spline;
565
566 // First, create multispline
567 multi_spline = create_multi_UBspline_2d_d (x_grid, y_grid, xBC, yBC,
568 num_splines);
569
570 double data[Nx*Ny];
571 // Now, create normal splines and set multispline data
572 for (int i=0; i<num_splines; i++) {
573 for (int j=0; j<Nx*Ny; j++)
574 data[j] = (drand48()-0.5);
575 norm_splines[i] = create_UBspline_2d_d (x_grid, y_grid, xBC, yBC, data);
576 set_multi_UBspline_2d_d (multi_spline, i, data);
577 }
578
579 // fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
580 // real(norm_splines[19]->coefs[227]),
581 // imag(norm_splines[19]->coefs[227]));
582 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
583 // real(multi_spline->coefs[19+227*multi_spline->z_stride]),
584 // imag(multi_spline->coefs[19+227*multi_spline->z_stride]));
585
586 // Now, test random values
587 int num_vals = 100;
588 double multi_vals[num_splines], norm_vals[num_splines];
589 double multi_grads[2*num_splines], norm_grads[2*num_splines];
590 double multi_lapl[num_splines], norm_lapl[num_splines];
591 double multi_hess[4*num_splines], norm_hess[4*num_splines];
592 for (int i=0; i<num_vals; i++) {
593 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
594 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
595
596
597 //////////////////////////
598 // Check value routine //
599 //////////////////////////
600 eval_multi_UBspline_2d_d (multi_spline, x, y, multi_vals);
601 for (int j=0; j<num_splines; j++)
602 eval_UBspline_2d_d (norm_splines[j], x, y, &(norm_vals[j]));
603 for (int j=0; j<num_splines; j++) {
604 // Check value
605 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
606 return -1;
607 }
608
609 ///////////////////////
610 // Check VG routine //
611 ///////////////////////
612 eval_multi_UBspline_2d_d_vg (multi_spline, x, y,
613 multi_vals, multi_grads);
614 for (int j=0; j<num_splines; j++)
615 eval_UBspline_2d_d_vg (norm_splines[j], x, y, &(norm_vals[j]),
616 &(norm_grads[2*j]));
617 for (int j=0; j<num_splines; j++) {
618 // Check value
619 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
620 return -1;
621
622 // Check gradients
623 for (int n=0; n<2; n++)
624 if (diff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-12))
625 return -2;
626 }
627
628
629 ///////////////////////
630 // Check VGL routine //
631 ///////////////////////
632 eval_multi_UBspline_2d_d_vgl (multi_spline, x, y,
633 multi_vals, multi_grads, multi_lapl);
634 for (int j=0; j<num_splines; j++)
635 eval_UBspline_2d_d_vgl (norm_splines[j], x, y, &(norm_vals[j]),
636 &(norm_grads[2*j]), &(norm_lapl[j]));
637 for (int j=0; j<num_splines; j++) {
638 // Check value
639 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
640 return -3;
641
642 // Check gradients
643 for (int n=0; n<2; n++)
644 if (diff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-10))
645 return -4;
646
647 // Check laplacian
648 if (diff (norm_lapl[j], multi_lapl[j], 1.0e-10))
649 return -5;
650 }
651
652
653 ///////////////////////
654 // Check VGH routine //
655 ///////////////////////
656 eval_multi_UBspline_2d_d_vgh (multi_spline, x, y,
657 multi_vals, multi_grads, multi_hess);
658 for (int j=0; j<num_splines; j++)
659 eval_UBspline_2d_d_vgh (norm_splines[j], x, y, &(norm_vals[j]),
660 &(norm_grads[2*j]), &(norm_hess[4*j]));
661 for (int j=0; j<num_splines; j++) {
662 // Check value
663 if (diff(norm_vals[j], multi_vals[j], 1.0e-12)) {
664 fprintf (stderr, "j = %d\n", j);
665 fprintf (stderr, "norm_vals[j] = %1.14e\n", norm_vals[j]);
666 fprintf (stderr, "multi_vals[j] = %1.14e\n", multi_vals[j]);
667 //return -6;
668 }
669
670 // Check gradients
671 for (int n=0; n<2; n++)
672 if (diff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-12))
673 return -7;
674
675 // Check hessian
676 for (int n=0; n<4; n++)
677 if (diff (norm_hess[4*j+n], multi_hess[4*j+n], 1.0e-10)) {
678 fprintf (stderr, "j = %d n = %d \n", j, n);
679 fprintf (stderr, "norm_hess[j] = %1.14e\n", norm_hess[4*j+n]);
680 fprintf (stderr, "multi_hess[j] = %1.14e\n", multi_hess[4*j+n]);
681 //return -8;
682 }
683 }
684 }
685 return 0;
686 }
687
688
689 int
test_3d_double_all()690 test_3d_double_all()
691 {
692 int Nx=73; int Ny=91; int Nz = 29;
693 int num_splines = 21;
694
695 Ugrid x_grid, y_grid, z_grid;
696 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
697 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
698 z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
699
700 BCtype_d xBC, yBC, zBC;
701 xBC.lCode = xBC.rCode = PERIODIC;
702 yBC.lCode = yBC.rCode = PERIODIC;
703 zBC.lCode = zBC.rCode = PERIODIC;
704
705 // First, create splines the normal way
706 UBspline_3d_d* norm_splines[num_splines];
707 multi_UBspline_3d_d *multi_spline;
708
709 // First, create multispline
710 multi_spline = create_multi_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC,
711 num_splines);
712
713 double data[Nx*Ny*Nz];
714 // Now, create normal splines and set multispline data
715 for (int i=0; i<num_splines; i++) {
716 for (int j=0; j<Nx*Ny*Nz; j++)
717 data[j] = (drand48()-0.5);
718 norm_splines[i] = create_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
719 set_multi_UBspline_3d_d (multi_spline, i, data);
720 }
721
722 // fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
723 // real(norm_splines[19]->coefs[227]),
724 // imag(norm_splines[19]->coefs[227]));
725 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
726 // real(multi_spline->coefs[19+227*multi_spline->z_stride]),
727 // imag(multi_spline->coefs[19+227*multi_spline->z_stride]));
728
729 // Now, test random values
730 int num_vals = 100;
731 double multi_vals[num_splines], norm_vals[num_splines];
732 double multi_grads[3*num_splines], norm_grads[3*num_splines];
733 double multi_lapl[num_splines], norm_lapl[num_splines];
734 double multi_hess[9*num_splines], norm_hess[9*num_splines];
735 for (int i=0; i<num_vals; i++) {
736 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
737 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
738 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
739
740
741 ///////////////////////
742 // Check VG routine //
743 ///////////////////////
744 eval_multi_UBspline_3d_d_vg (multi_spline, x, y, z,
745 multi_vals, multi_grads);
746 for (int j=0; j<num_splines; j++)
747 eval_UBspline_3d_d_vg (norm_splines[j], x, y, z, &(norm_vals[j]),
748 &(norm_grads[3*j]));
749 for (int j=0; j<num_splines; j++) {
750 // Check value
751 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
752 return -1;
753
754 // Check gradients
755 for (int n=0; n<3; n++)
756 if (diff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-12))
757 return -2;
758 }
759
760
761 ///////////////////////
762 // Check VGL routine //
763 ///////////////////////
764 eval_multi_UBspline_3d_d_vgl (multi_spline, x, y, z,
765 multi_vals, multi_grads, multi_lapl);
766 for (int j=0; j<num_splines; j++)
767 eval_UBspline_3d_d_vgl (norm_splines[j], x, y, z, &(norm_vals[j]),
768 &(norm_grads[3*j]), &(norm_lapl[j]));
769 for (int j=0; j<num_splines; j++) {
770 // Check value
771 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
772 return -3;
773
774 // Check gradients
775 for (int n=0; n<3; n++)
776 if (diff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-10))
777 return -4;
778
779 // Check laplacian
780 if (diff (norm_lapl[j], multi_lapl[j], 1.0e-10))
781 return -5;
782 }
783
784
785 ///////////////////////
786 // Check VGH routine //
787 ///////////////////////
788 eval_multi_UBspline_3d_d_vgh (multi_spline, x, y, z,
789 multi_vals, multi_grads, multi_hess);
790 for (int j=0; j<num_splines; j++)
791 eval_UBspline_3d_d_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
792 &(norm_grads[3*j]), &(norm_hess[9*j]));
793 for (int j=0; j<num_splines; j++) {
794 // Check value
795 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
796 return -6;
797
798 // Check gradients
799 for (int n=0; n<3; n++)
800 if (diff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-12))
801 return -7;
802
803 // Check hessian
804 for (int n=0; n<9; n++)
805 if (diff (norm_hess[9*j+n], multi_hess[9*j+n], 1.0e-10))
806 return -8;
807 }
808 }
809 return 0;
810 }
811
812
813
814
815 /////////////////////////////////////////////
816 // Single-precision complex test functions //
817 /////////////////////////////////////////////
818 inline int
cdiff(complex_float a,complex_float b,double tol)819 cdiff (complex_float a, complex_float b, double tol)
820 {
821 double rdiff = fabs(real(a) - real(b));
822 double idiff = fabs(imag(a) - imag(b));
823 if (rdiff > tol || idiff > tol)
824 return 1;
825 else
826 return 0;
827 }
828
829 int
test_1d_complex_float_all()830 test_1d_complex_float_all()
831 {
832 int Nx=73;
833 int num_splines = 21;
834
835 Ugrid x_grid;
836 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
837
838 BCtype_c xBC;
839 xBC.lCode = xBC.rCode = PERIODIC;
840
841 // First, create splines the normal way
842 UBspline_1d_c* norm_splines[num_splines];
843 multi_UBspline_1d_c *multi_spline;
844
845 // First, create multispline
846 multi_spline = create_multi_UBspline_1d_c (x_grid, xBC, num_splines);
847
848 complex_float data[Nx];
849 // Now, create normal splines and set multispline data
850 for (int i=0; i<num_splines; i++) {
851 for (int j=0; j<Nx; j++)
852 data[j] = complex<float>((drand48()-0.5),(drand48()-0.5));
853 norm_splines[i] = create_UBspline_1d_c (x_grid, xBC, data);
854 set_multi_UBspline_1d_c (multi_spline, i, data);
855 }
856
857 // fprintf (stderr, "\nnorm coef = %1.14e + %1.14ei\n",
858 // real(norm_splines[19]->coefs[27]),
859 // imag(norm_splines[19]->coefs[27]));
860 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
861 // real(multi_spline->coefs[19+27*multi_spline->x_stride]),
862 // imag(multi_spline->coefs[19+27*multi_spline->x_stride]));
863
864
865 // Now, test random values
866 int num_vals = 100;
867 complex_float multi_vals[num_splines], norm_vals [num_splines];
868 complex_float multi_grads[num_splines], norm_grads[num_splines];
869 complex_float multi_lapl[num_splines], norm_lapl [num_splines];
870 for (int i=0; i<num_vals; i++) {
871 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
872
873 //////////////////////////
874 // Check value routine //
875 //////////////////////////
876 eval_multi_UBspline_1d_c (multi_spline, x, multi_vals);
877 for (int j=0; j<num_splines; j++)
878 eval_UBspline_1d_c (norm_splines[j], x, &(norm_vals[j]));
879 for (int j=0; j<num_splines; j++) {
880 // Check value
881 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-6)) {
882 fprintf (stderr, " j = %d\n", j);
883 fprintf (stderr, " norm_vals[j] = %1.14e + %1.14ei\n",
884 real (norm_vals[j]), imag(norm_vals[j]));
885 fprintf (stderr, "multi_vals[j] = %1.14e + %1.14ei\n",
886 real (multi_vals[j]), imag(multi_vals[j]));
887
888 return -1;
889 }
890 }
891
892 ///////////////////////
893 // Check VG routine //
894 ///////////////////////
895 eval_multi_UBspline_1d_c_vg (multi_spline, x,
896 multi_vals, multi_grads);
897 for (int j=0; j<num_splines; j++)
898 eval_UBspline_1d_c_vg (norm_splines[j], x, &(norm_vals[j]),
899 &(norm_grads[j]));
900 for (int j=0; j<num_splines; j++) {
901 // Check value
902 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-6))
903 return -1;
904
905 // Check gradients
906 if (cdiff (norm_grads[j], multi_grads[j], 1.0e-5))
907 return -2;
908 }
909
910
911 ///////////////////////
912 // Check VGL routine //
913 ///////////////////////
914 eval_multi_UBspline_1d_c_vgl (multi_spline, x, multi_vals, multi_grads, multi_lapl);
915 for (int j=0; j<num_splines; j++)
916 eval_UBspline_1d_c_vgl (norm_splines[j], x, &(norm_vals[j]),
917 &(norm_grads[j]), &(norm_lapl[j]));
918 for (int j=0; j<num_splines; j++) {
919 // Check value
920 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-6))
921 return -3;
922
923 // Check gradients
924 if (cdiff (norm_grads[j], multi_grads[j], 1.0e-5))
925 return -4;
926
927 // Check laplacian
928 if (cdiff (norm_lapl[j], multi_lapl[j], 1.0e-3))
929 return -5;
930 }
931 }
932 return 0;
933 }
934
935
936 int
test_2d_complex_float_all()937 test_2d_complex_float_all()
938 {
939 int Nx=73; int Ny=91;
940 int num_splines = 20;
941
942 Ugrid x_grid, y_grid;
943 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
944 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
945
946 BCtype_c xBC, yBC;
947 xBC.lCode = xBC.rCode = PERIODIC;
948 yBC.lCode = yBC.rCode = PERIODIC;
949
950 // First, create splines the normal way
951 UBspline_2d_c* norm_splines[num_splines];
952 multi_UBspline_2d_c *multi_spline;
953
954 // First, create multispline
955 multi_spline = create_multi_UBspline_2d_c (x_grid, y_grid, xBC, yBC,
956 num_splines);
957
958 complex_float data[Nx*Ny];
959 // Now, create normal splines and set multispline data
960 for (int i=0; i<num_splines; i++) {
961 for (int j=0; j<Nx*Ny; j++)
962 data[j] = complex<float>((drand48()-0.5),(drand48()-0.5));
963 norm_splines[i] = create_UBspline_2d_c (x_grid, y_grid, xBC, yBC, data);
964 set_multi_UBspline_2d_c (multi_spline, i, data);
965 }
966
967 // fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
968 // real(norm_splines[19]->coefs[2127]),
969 // imag(norm_splines[19]->coefs[2127]));
970 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
971 // real(multi_spline->coefs[19+2127*multi_spline->y_stride]),
972 // imag(multi_spline->coefs[19+2127*multi_spline->y_stride]));
973
974 // Now, test random values
975 int num_vals = 100;
976 complex_float multi_vals[num_splines], norm_vals[num_splines];
977 complex_float multi_grads[2*num_splines], norm_grads[2*num_splines];
978 complex_float multi_lapl[num_splines], norm_lapl[num_splines];
979 complex_float multi_hess[4*num_splines], norm_hess[4*num_splines];
980 for (int i=0; i<num_vals; i++) {
981 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
982 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
983
984
985 //////////////////////////
986 // Check value routine //
987 //////////////////////////
988 eval_multi_UBspline_2d_c (multi_spline, x, y, multi_vals);
989 for (int j=0; j<num_splines; j++)
990 eval_UBspline_2d_c (norm_splines[j], x, y, &(norm_vals[j]));
991 for (int j=0; j<num_splines; j++) {
992 // Check value
993 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-5))
994 return -1;
995 }
996
997 ///////////////////////
998 // Check VG routine //
999 ///////////////////////
1000 eval_multi_UBspline_2d_c_vg (multi_spline, x, y,
1001 multi_vals, multi_grads);
1002 for (int j=0; j<num_splines; j++)
1003 eval_UBspline_2d_c_vg (norm_splines[j], x, y, &(norm_vals[j]),
1004 &(norm_grads[2*j]));
1005 for (int j=0; j<num_splines; j++) {
1006 // Check value
1007 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-5)) {
1008 fprintf (stderr, " norm_vals[j] = %1.8f + %1.8fi\n",
1009 real(norm_vals[j]), imag(norm_vals[j]));
1010 fprintf (stderr, "multi_vals[j] = %1.8f + %1.8fi\n",
1011 real(multi_vals[j]), imag(multi_vals[j]));
1012 return -2;
1013 }
1014
1015 // Check gradients
1016 for (int n=0; n<2; n++)
1017 if (cdiff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-3)) {
1018 fprintf (stderr, "norm_grads[j] = %1.14e + %1.14ei\n",
1019 real(norm_grads[2*j+n]), imag(norm_grads[2*j+n]));
1020 fprintf (stderr, "multi_grads[j] = %1.14e + %1.14ei\n",
1021 real(multi_grads[2*j+n]), imag(multi_grads[2*j+n]));
1022 return -3;
1023 }
1024 }
1025
1026
1027 ///////////////////////
1028 // Check VGL routine //
1029 ///////////////////////
1030 eval_multi_UBspline_2d_c_vgl (multi_spline, x, y,
1031 multi_vals, multi_grads, multi_lapl);
1032 for (int j=0; j<num_splines; j++)
1033 eval_UBspline_2d_c_vgl (norm_splines[j], x, y, &(norm_vals[j]),
1034 &(norm_grads[2*j]), &(norm_lapl[j]));
1035 for (int j=0; j<num_splines; j++) {
1036 // Check value
1037 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-5))
1038 return -4;
1039
1040 // Check gradients
1041 for (int n=0; n<2; n++)
1042 if (cdiff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-3))
1043 return -5;
1044
1045 // Check laplacian
1046 if (cdiff (norm_lapl[j], multi_lapl[j], 1.0e-2)) {
1047 fprintf (stderr, "norm_lapl[j] = %1.6f + %1.6fi\n",
1048 real(norm_lapl[j]), imag(norm_lapl[j]));
1049 fprintf (stderr, "multi_lapl[j] = %1.6f + %1.6fi\n",
1050 real(multi_lapl[j]), imag(multi_lapl[j]));
1051 return -6;
1052 }
1053 }
1054
1055
1056 ///////////////////////
1057 // Check VGH routine //
1058 ///////////////////////
1059 eval_multi_UBspline_2d_c_vgh (multi_spline, x, y,
1060 multi_vals, multi_grads, multi_hess);
1061 for (int j=0; j<num_splines; j++)
1062 eval_UBspline_2d_c_vgh (norm_splines[j], x, y, &(norm_vals[j]),
1063 &(norm_grads[2*j]), &(norm_hess[4*j]));
1064 for (int j=0; j<num_splines; j++) {
1065 // Check value
1066 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-5)) {
1067 fprintf (stderr, "j = %d\n", j);
1068 fprintf (stderr, "norm_vals[j] = %1.14e + %1.14ei\n",
1069 real(norm_vals[j]), imag(norm_vals[j]));
1070 fprintf (stderr, "multi_vals[j] = %1.14e + %1.14ei\n",
1071 real(multi_vals[j]), imag(multi_vals[j]));
1072 return -7;
1073 }
1074
1075 // Check gradients
1076 for (int n=0; n<2; n++)
1077 if (cdiff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-3)) {
1078 fprintf (stderr, "j = %d\n", j);
1079 fprintf (stderr, "norm_grads[j] = %1.14e + %1.14ei\n",
1080 real(norm_grads[2*j+n]), imag(norm_grads[2*j+n]));
1081 fprintf (stderr, "multi_grads[j] = %1.14e + %1.14ei\n",
1082 real(multi_grads[2*j+n]), imag(multi_grads[2*j+n]));
1083 return -8;
1084 }
1085
1086
1087 // Check hessian
1088 for (int n=0; n<4; n++)
1089 if (cdiff (norm_hess[4*j+n], multi_hess[4*j+n], 1.0e-2)) {
1090 fprintf (stderr, "\nj = %d n = %d \n", j, n);
1091 fprintf (stderr, "norm_hess[j] = %1.6f + %1.6fi\n",
1092 real(norm_hess[4*j+n]), imag(norm_hess[4*j+n]));
1093 fprintf (stderr, "multi_hess[j] = %1.6f + %1.6fi\n",
1094 real(multi_hess[4*j+n]), imag(multi_hess[4*j+n]));
1095 return -9;
1096 }
1097 }
1098 }
1099 return 0;
1100 }
1101
1102
1103 int
test_3d_complex_float_all()1104 test_3d_complex_float_all()
1105 {
1106 int Nx=73; int Ny=91; int Nz = 29;
1107 int num_splines = 21;
1108
1109 Ugrid x_grid, y_grid, z_grid;
1110 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
1111 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
1112 z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
1113
1114 BCtype_c xBC, yBC, zBC;
1115 xBC.lCode = xBC.rCode = PERIODIC;
1116 yBC.lCode = yBC.rCode = PERIODIC;
1117 zBC.lCode = zBC.rCode = PERIODIC;
1118
1119 // First, create splines the normal way
1120 UBspline_3d_c* norm_splines[num_splines];
1121 multi_UBspline_3d_c *multi_spline;
1122
1123 // First, create multispline
1124 multi_spline = create_multi_UBspline_3d_c (x_grid, y_grid, z_grid, xBC, yBC, zBC,
1125 num_splines);
1126
1127 complex_float data[Nx*Ny*Nz];
1128 // Now, create normal splines and set multispline data
1129 for (int i=0; i<num_splines; i++) {
1130 for (int j=0; j<Nx*Ny*Nz; j++)
1131 data[j] = complex<float>((drand48()-0.5), (drand48()-0.5));
1132 norm_splines[i] = create_UBspline_3d_c (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
1133 set_multi_UBspline_3d_c (multi_spline, i, data);
1134 }
1135
1136 // fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
1137 // real(norm_splines[19]->coefs[227]),
1138 // imag(norm_splines[19]->coefs[227]));
1139 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
1140 // real(multi_spline->coefs[19+227*multi_spline->z_stride]),
1141 // imag(multi_spline->coefs[19+227*multi_spline->z_stride]));
1142
1143 // Now, test random values
1144 int num_vals = 100;
1145 complex_float multi_vals[num_splines], norm_vals[num_splines];
1146 complex_float multi_grads[3*num_splines], norm_grads[3*num_splines];
1147 complex_float multi_lapl[num_splines], norm_lapl[num_splines];
1148 complex_float multi_hess[9*num_splines], norm_hess[9*num_splines];
1149 for (int i=0; i<num_vals; i++) {
1150 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1151 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1152 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1153
1154 /////////////////////////
1155 // Check value routine //
1156 /////////////////////////
1157 eval_multi_UBspline_3d_c (multi_spline, x, y, z, multi_vals);
1158 for (int j=0; j<num_splines; j++)
1159 eval_UBspline_3d_c (norm_splines[j], x, y, z, &(norm_vals[j]));
1160 for (int j=0; j<num_splines; j++) {
1161 // Check value
1162 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-6))
1163 return -1;
1164 }
1165
1166 ///////////////////////
1167 // Check VG routine //
1168 ///////////////////////
1169 eval_multi_UBspline_3d_c_vg (multi_spline, x, y, z,
1170 multi_vals, multi_grads);
1171 for (int j=0; j<num_splines; j++)
1172 eval_UBspline_3d_c_vg (norm_splines[j], x, y, z, &(norm_vals[j]),
1173 &(norm_grads[3*j]));
1174 for (int j=0; j<num_splines; j++) {
1175 // Check value
1176 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-6))
1177 return -2;
1178
1179 // Check gradients
1180 for (int n=0; n<3; n++)
1181 if (cdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-4))
1182 return -3;
1183 }
1184
1185
1186 ///////////////////////
1187 // Check VGL routine //
1188 ///////////////////////
1189 eval_multi_UBspline_3d_c_vgl (multi_spline, x, y, z,
1190 multi_vals, multi_grads, multi_lapl);
1191 for (int j=0; j<num_splines; j++)
1192 eval_UBspline_3d_c_vgl (norm_splines[j], x, y, z, &(norm_vals[j]),
1193 &(norm_grads[3*j]), &(norm_lapl[j]));
1194 for (int j=0; j<num_splines; j++) {
1195 // Check value
1196 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-6))
1197 return -4;
1198
1199 // Check gradients
1200 for (int n=0; n<3; n++)
1201 if (cdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-4))
1202 return -5;
1203
1204 // Check laplacian
1205 if (cdiff (norm_lapl[j], multi_lapl[j], 1.0e-2))
1206 return -6;
1207 }
1208
1209
1210 ///////////////////////
1211 // Check VGH routine //
1212 ///////////////////////
1213 eval_multi_UBspline_3d_c_vgh (multi_spline, x, y, z,
1214 multi_vals, multi_grads, multi_hess);
1215 for (int j=0; j<num_splines; j++)
1216 eval_UBspline_3d_c_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
1217 &(norm_grads[3*j]), &(norm_hess[9*j]));
1218 for (int j=0; j<num_splines; j++) {
1219 // Check value
1220 if (cdiff(norm_vals[j], multi_vals[j], 1.0e-6))
1221 return -7;
1222
1223 // Check gradients
1224 for (int n=0; n<3; n++)
1225 if (cdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-4))
1226 return -8;
1227
1228 // Check hessian
1229 for (int n=0; n<9; n++)
1230 if (cdiff (norm_hess[9*j+n], multi_hess[9*j+n], 1.0e-2))
1231 return -9;
1232 }
1233 }
1234 return 0;
1235 }
1236
1237
1238
1239 /////////////////////////////////////////////
1240 // Double-precision complex test functions //
1241 /////////////////////////////////////////////
test_complex_double()1242 void test_complex_double()
1243 {
1244 int Nx=73; int Ny=91; int Nz = 29;
1245 int num_splines = 128;
1246
1247 Ugrid x_grid, y_grid, z_grid;
1248 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
1249 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
1250 z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
1251
1252 BCtype_z xBC, yBC, zBC;
1253 xBC.lCode = xBC.rCode = PERIODIC;
1254 yBC.lCode = yBC.rCode = PERIODIC;
1255 zBC.lCode = zBC.rCode = PERIODIC;
1256
1257 // First, create splines the normal way
1258 UBspline_3d_z* norm_splines[num_splines];
1259 multi_UBspline_3d_z *multi_spline;
1260
1261 // First, create multispline
1262 multi_spline = create_multi_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC,
1263 num_splines);
1264
1265 complex_double data[Nx*Ny*Nz];
1266 // Now, create normal splines and set multispline data
1267 for (int i=0; i<num_splines; i++) {
1268 for (int j=0; j<Nx*Ny*Nz; j++)
1269 data[j] = complex<double>((drand48()-0.5),(drand48()-0.5));
1270 norm_splines[i] = create_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
1271 set_multi_UBspline_3d_z (multi_spline, i, data);
1272 }
1273
1274 fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
1275 real(norm_splines[19]->coefs[227]),
1276 imag(norm_splines[19]->coefs[227]));
1277 fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
1278 real(multi_spline->coefs[19+227*multi_spline->z_stride]),
1279 imag(multi_spline->coefs[19+227*multi_spline->z_stride]));
1280 //return;
1281
1282 // Now, test random values
1283 int num_vals = 100;
1284 complex_double multi_vals[num_splines], norm_vals[num_splines];
1285 for (int i=0; i<num_vals; i++) {
1286 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1287 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1288 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1289
1290 eval_multi_UBspline_3d_z (multi_spline, x, y, z, multi_vals);
1291 for (int j=0; j<num_splines; j++)
1292 eval_UBspline_3d_z (norm_splines[j], x, y, z, &(norm_vals[j]));
1293 for (int j=0; j<num_splines; j++) {
1294 double rdiff = real(norm_vals[j]) - real(multi_vals[j]);
1295 double idiff = imag(norm_vals[j]) - imag(multi_vals[j]);
1296 if (fabs(rdiff) > 1.0e-12 || fabs(idiff) > 1.0e-12) {
1297 fprintf (stderr, "Error! norm_vals[j] = %1.14e + %1.14ei\n",
1298 real(norm_vals[j]), imag(norm_vals[j]));
1299 fprintf (stderr, " multi_vals[j] = %1.14e + %1.14ei\n",
1300 real(multi_vals[j]), imag(multi_vals[j]));
1301 }
1302 }
1303 }
1304
1305 num_vals = 100000;
1306
1307 // Now do timing
1308 clock_t norm_start, norm_end, multi_start, multi_end, rand_start, rand_end;
1309 rand_start = clock();
1310 for (int i=0; i<num_vals; i++) {
1311 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1312 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1313 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1314 }
1315 rand_end = clock();
1316
1317 norm_start = clock();
1318 for (int i=0; i<num_vals; i++) {
1319 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1320 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1321 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1322
1323 for (int j=0; j<num_splines; j++)
1324 eval_UBspline_3d_z (norm_splines[j], x, y, z, &(norm_vals[j]));
1325 }
1326 norm_end = clock();
1327
1328 multi_start = clock();
1329 for (int i=0; i<num_vals; i++) {
1330 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1331 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1332 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1333 eval_multi_UBspline_3d_z (multi_spline, x, y, z, multi_vals);
1334 }
1335 multi_end = clock();
1336
1337 fprintf (stderr, "Normal spline time = %1.5f\n",
1338 (double)(norm_end-norm_start+rand_start-rand_end)/CLOCKS_PER_SEC);
1339 fprintf (stderr, "Multi spline time = %1.5f\n",
1340 (double)(multi_end-multi_start+rand_start-rand_end)/CLOCKS_PER_SEC);
1341
1342 }
1343
1344
1345 inline int
zdiff(complex_double a,complex_double b,double tol)1346 zdiff (complex_double a, complex_double b, double tol)
1347 {
1348 double rdiff = fabs(real(a) - real(b));
1349 double idiff = fabs(imag(a) - imag(b));
1350 if (rdiff > tol || idiff > tol)
1351 return 1;
1352 else
1353 return 0;
1354 }
1355
1356
1357 int
test_1d_complex_double_all()1358 test_1d_complex_double_all()
1359 {
1360 int Nx=73;
1361 int num_splines = 21;
1362
1363 Ugrid x_grid;
1364 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
1365
1366 BCtype_z xBC;
1367 xBC.lCode = xBC.rCode = PERIODIC;
1368
1369 // First, create splines the normal way
1370 UBspline_1d_z* norm_splines[num_splines];
1371 multi_UBspline_1d_z *multi_spline;
1372
1373 // First, create multispline
1374 multi_spline = create_multi_UBspline_1d_z (x_grid, xBC, num_splines);
1375
1376 complex_double data[Nx];
1377 // Now, create normal splines and set multispline data
1378 for (int i=0; i<num_splines; i++) {
1379 for (int j=0; j<Nx; j++)
1380 data[j] = complex<double>((drand48()-0.5), (drand48()-0.5));
1381 norm_splines[i] = create_UBspline_1d_z (x_grid, xBC, data);
1382 set_multi_UBspline_1d_z (multi_spline, i, data);
1383 }
1384
1385 // fprintf (stderr, "\nnorm coef = %1.14e + %1.14ei\n",
1386 // real(norm_splines[19]->coefs[27]),
1387 // imag(norm_splines[19]->coefs[27]));
1388 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
1389 // real(multi_spline->coefs[19+27*multi_spline->x_stride]),
1390 // imag(multi_spline->coefs[19+27*multi_spline->x_stride]));
1391
1392
1393 // Now, test random values
1394 int num_vals = 100;
1395 complex_double multi_vals[num_splines], norm_vals [num_splines];
1396 complex_double multi_grads[num_splines], norm_grads[num_splines];
1397 complex_double multi_lapl[num_splines], norm_lapl [num_splines];
1398 for (int i=0; i<num_vals; i++) {
1399 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1400
1401 //////////////////////////
1402 // Check value routine //
1403 //////////////////////////
1404 eval_multi_UBspline_1d_z (multi_spline, x, multi_vals);
1405 for (int j=0; j<num_splines; j++)
1406 eval_UBspline_1d_z (norm_splines[j], x, &(norm_vals[j]));
1407 for (int j=0; j<num_splines; j++) {
1408 // Check value
1409 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12)) {
1410 fprintf (stderr, " norm_vals[j] = %1.14e + %1.14ei\n",
1411 real (norm_vals[j]), imag(norm_vals[j]));
1412 fprintf (stderr, "multi_vals[j] = %1.14e + %1.14ei\n",
1413 real (multi_vals[j]), imag(multi_vals[j]));
1414
1415 return -1;
1416 }
1417 }
1418
1419 ///////////////////////
1420 // Check VG routine //
1421 ///////////////////////
1422 eval_multi_UBspline_1d_z_vg (multi_spline, x,
1423 multi_vals, multi_grads);
1424 for (int j=0; j<num_splines; j++)
1425 eval_UBspline_1d_z_vg (norm_splines[j], x, &(norm_vals[j]),
1426 &(norm_grads[j]));
1427 for (int j=0; j<num_splines; j++) {
1428 // Check value
1429 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
1430 return -1;
1431
1432 // Check gradients
1433 if (zdiff (norm_grads[j], multi_grads[j], 1.0e-12))
1434 return -2;
1435 }
1436
1437
1438 ///////////////////////
1439 // Check VGL routine //
1440 ///////////////////////
1441 eval_multi_UBspline_1d_z_vgl (multi_spline, x, multi_vals, multi_grads, multi_lapl);
1442 for (int j=0; j<num_splines; j++)
1443 eval_UBspline_1d_z_vgl (norm_splines[j], x, &(norm_vals[j]),
1444 &(norm_grads[j]), &(norm_lapl[j]));
1445 for (int j=0; j<num_splines; j++) {
1446 // Check value
1447 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
1448 return -3;
1449
1450 // Check gradients
1451 if (zdiff (norm_grads[j], multi_grads[j], 1.0e-10))
1452 return -4;
1453
1454 // Check laplacian
1455 if (zdiff (norm_lapl[j], multi_lapl[j], 1.0e-10))
1456 return -5;
1457 }
1458 }
1459 return 0;
1460 }
1461
1462
1463 int
test_2d_complex_double_all()1464 test_2d_complex_double_all()
1465 {
1466 int Nx=73; int Ny=91;
1467 int num_splines = 21;
1468
1469 Ugrid x_grid, y_grid;
1470 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
1471 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
1472
1473 BCtype_z xBC, yBC;
1474 xBC.lCode = xBC.rCode = PERIODIC;
1475 yBC.lCode = yBC.rCode = PERIODIC;
1476
1477 // First, create splines the normal way
1478 UBspline_2d_z* norm_splines[num_splines];
1479 multi_UBspline_2d_z *multi_spline;
1480
1481 // First, create multispline
1482 multi_spline = create_multi_UBspline_2d_z (x_grid, y_grid, xBC, yBC,
1483 num_splines);
1484
1485 complex_double data[Nx*Ny];
1486 // Now, create normal splines and set multispline data
1487 for (int i=0; i<num_splines; i++) {
1488 for (int j=0; j<Nx*Ny; j++)
1489 data[j] = complex<double>((drand48()-0.5),(drand48()-0.5));
1490 norm_splines[i] = create_UBspline_2d_z (x_grid, y_grid, xBC, yBC, data);
1491 set_multi_UBspline_2d_z (multi_spline, i, data);
1492 }
1493
1494 // fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
1495 // real(norm_splines[19]->coefs[227]),
1496 // imag(norm_splines[19]->coefs[227]));
1497 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
1498 // real(multi_spline->coefs[19+227*multi_spline->y_stride]),
1499 // imag(multi_spline->coefs[19+227*multi_spline->y_stride]));
1500
1501 // Now, test random values
1502 int num_vals = 100;
1503 complex_double multi_vals[num_splines], norm_vals[num_splines];
1504 complex_double multi_grads[2*num_splines], norm_grads[2*num_splines];
1505 complex_double multi_lapl[num_splines], norm_lapl[num_splines];
1506 complex_double multi_hess[4*num_splines], norm_hess[4*num_splines];
1507 for (int i=0; i<num_vals; i++) {
1508 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1509 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1510
1511
1512 //////////////////////////
1513 // Check value routine //
1514 //////////////////////////
1515 eval_multi_UBspline_2d_z (multi_spline, x, y, multi_vals);
1516 for (int j=0; j<num_splines; j++)
1517 eval_UBspline_2d_z (norm_splines[j], x, y, &(norm_vals[j]));
1518 for (int j=0; j<num_splines; j++) {
1519 // Check value
1520 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
1521 return -1;
1522 }
1523
1524 ///////////////////////
1525 // Check VG routine //
1526 ///////////////////////
1527 eval_multi_UBspline_2d_z_vg (multi_spline, x, y,
1528 multi_vals, multi_grads);
1529 for (int j=0; j<num_splines; j++)
1530 eval_UBspline_2d_z_vg (norm_splines[j], x, y, &(norm_vals[j]),
1531 &(norm_grads[2*j]));
1532 for (int j=0; j<num_splines; j++) {
1533 // Check value
1534 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
1535 return -1;
1536
1537 // Check gradients
1538 for (int n=0; n<2; n++)
1539 if (zdiff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-12))
1540 return -2;
1541 }
1542
1543
1544 ///////////////////////
1545 // Check VGL routine //
1546 ///////////////////////
1547 eval_multi_UBspline_2d_z_vgl (multi_spline, x, y,
1548 multi_vals, multi_grads, multi_lapl);
1549 for (int j=0; j<num_splines; j++)
1550 eval_UBspline_2d_z_vgl (norm_splines[j], x, y, &(norm_vals[j]),
1551 &(norm_grads[2*j]), &(norm_lapl[j]));
1552 for (int j=0; j<num_splines; j++) {
1553 // Check value
1554 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
1555 return -3;
1556
1557 // Check gradients
1558 for (int n=0; n<2; n++)
1559 if (zdiff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-10))
1560 return -4;
1561
1562 // Check laplacian
1563 if (zdiff (norm_lapl[j], multi_lapl[j], 1.0e-9)) {
1564 fprintf (stderr, "norm_lapl[j] = %1.14e + %1.14ei\n",
1565 real(norm_lapl[j]), imag(norm_lapl[j]));
1566 fprintf (stderr, "multi_lapl[j] = %1.14e + %1.14ei\n",
1567 real(multi_lapl[j]), imag(multi_lapl[j]));
1568 return -5;
1569 }
1570 }
1571
1572
1573 ///////////////////////
1574 // Check VGH routine //
1575 ///////////////////////
1576 eval_multi_UBspline_2d_z_vgh (multi_spline, x, y,
1577 multi_vals, multi_grads, multi_hess);
1578 for (int j=0; j<num_splines; j++)
1579 eval_UBspline_2d_z_vgh (norm_splines[j], x, y, &(norm_vals[j]),
1580 &(norm_grads[2*j]), &(norm_hess[4*j]));
1581 for (int j=0; j<num_splines; j++) {
1582 // Check value
1583 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12)) {
1584 fprintf (stderr, "j = %d\n", j);
1585 fprintf (stderr, "norm_vals[j] = %1.14e + %1.14ei\n",
1586 real(norm_vals[j]), imag(norm_vals[j]));
1587 fprintf (stderr, "multi_vals[j] = %1.14e + %1.14ei\n",
1588 real(multi_vals[j]), imag(multi_vals[j]));
1589 return -6;
1590 }
1591
1592 // Check gradients
1593 for (int n=0; n<2; n++)
1594 if (zdiff (norm_grads[2*j+n], multi_grads[2*j+n], 1.0e-12))
1595 return -7;
1596
1597 // Check hessian
1598 for (int n=0; n<4; n++)
1599 if (zdiff (norm_hess[4*j+n], multi_hess[4*j+n], 1.0e-10)) {
1600 fprintf (stderr, "j = %d n = %d \n", j, n);
1601 fprintf (stderr, "norm_hess[j] = %1.14e + %1.14ei\n",
1602 real(norm_hess[4*j+n]), imag(norm_hess[4*j+n]));
1603 fprintf (stderr, "multi_hess[j] = %1.14e + %1.15ei\n",
1604 real(multi_hess[4*j+n]), imag(multi_hess[4*j+n]));
1605 return -8;
1606 }
1607 }
1608 }
1609 return 0;
1610 }
1611
1612
1613 int
test_3d_complex_double_all()1614 test_3d_complex_double_all()
1615 {
1616 int Nx=73; int Ny=91; int Nz = 29;
1617 int num_splines = 21;
1618
1619 Ugrid x_grid, y_grid, z_grid;
1620 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
1621 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
1622 z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
1623
1624 BCtype_z xBC, yBC, zBC;
1625 xBC.lCode = xBC.rCode = PERIODIC;
1626 yBC.lCode = yBC.rCode = PERIODIC;
1627 zBC.lCode = zBC.rCode = PERIODIC;
1628
1629 // First, create splines the normal way
1630 UBspline_3d_z* norm_splines[num_splines];
1631 multi_UBspline_3d_z *multi_spline;
1632
1633 // First, create multispline
1634 multi_spline = create_multi_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC,
1635 num_splines);
1636
1637 complex_double data[Nx*Ny*Nz];
1638 // Now, create normal splines and set multispline data
1639 for (int i=0; i<num_splines; i++) {
1640 for (int j=0; j<Nx*Ny*Nz; j++)
1641 data[j] = complex<double>((drand48()-0.5),(drand48()-0.5));
1642 norm_splines[i] = create_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
1643 set_multi_UBspline_3d_z (multi_spline, i, data);
1644 }
1645
1646 // fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
1647 // real(norm_splines[19]->coefs[227]),
1648 // imag(norm_splines[19]->coefs[227]));
1649 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
1650 // real(multi_spline->coefs[19+227*multi_spline->z_stride]),
1651 // imag(multi_spline->coefs[19+227*multi_spline->z_stride]));
1652
1653 // Now, test random values
1654 int num_vals = 100;
1655 complex_double multi_vals[num_splines], norm_vals[num_splines];
1656 complex_double multi_grads[3*num_splines], norm_grads[3*num_splines];
1657 complex_double multi_lapl[num_splines], norm_lapl[num_splines];
1658 complex_double multi_hess[9*num_splines], norm_hess[9*num_splines];
1659 for (int i=0; i<num_vals; i++) {
1660 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1661 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1662 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1663
1664 ///////////////////////
1665 // Check value only //
1666 ///////////////////////
1667 eval_multi_UBspline_3d_z (multi_spline, x, y, z, multi_vals);
1668 for (int j=0; j<num_splines; j++)
1669 eval_UBspline_3d_z (norm_splines[j], x, y, z, &(norm_vals[j]));
1670 for (int j=0; j<num_splines; j++)
1671 // Check value
1672 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
1673 return -1;
1674
1675 ///////////////////////
1676 // Check VG routine //
1677 ///////////////////////
1678 eval_multi_UBspline_3d_z_vg (multi_spline, x, y, z,
1679 multi_vals, multi_grads);
1680 for (int j=0; j<num_splines; j++)
1681 eval_UBspline_3d_z_vg (norm_splines[j], x, y, z, &(norm_vals[j]),
1682 &(norm_grads[3*j]));
1683 for (int j=0; j<num_splines; j++) {
1684 // Check value
1685 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
1686 return -2;
1687
1688 // Check gradients
1689 for (int n=0; n<3; n++)
1690 if (zdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-12))
1691 return -3;
1692 }
1693
1694
1695 ///////////////////////
1696 // Check VGL routine //
1697 ///////////////////////
1698 eval_multi_UBspline_3d_z_vgl (multi_spline, x, y, z,
1699 multi_vals, multi_grads, multi_lapl);
1700 for (int j=0; j<num_splines; j++)
1701 eval_UBspline_3d_z_vgl (norm_splines[j], x, y, z, &(norm_vals[j]),
1702 &(norm_grads[3*j]), &(norm_lapl[j]));
1703 for (int j=0; j<num_splines; j++) {
1704 // Check value
1705 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
1706 return -4;
1707
1708 // Check gradients
1709 for (int n=0; n<3; n++)
1710 if (zdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-10))
1711 return -5;
1712
1713 // Check laplacian
1714 if (zdiff (norm_lapl[j], multi_lapl[j], 1.0e-10))
1715 return -6;
1716 }
1717
1718
1719 ///////////////////////
1720 // Check VGH routine //
1721 ///////////////////////
1722 eval_multi_UBspline_3d_z_vgh (multi_spline, x, y, z,
1723 multi_vals, multi_grads, multi_hess);
1724 for (int j=0; j<num_splines; j++)
1725 eval_UBspline_3d_z_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
1726 &(norm_grads[3*j]), &(norm_hess[9*j]));
1727 for (int j=0; j<num_splines; j++) {
1728 // Check value
1729 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
1730 return -7;
1731
1732 // Check gradients
1733 for (int n=0; n<3; n++)
1734 if (zdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-12))
1735 return -8;
1736
1737 // Check hessian
1738 for (int n=0; n<9; n++)
1739 if (zdiff (norm_hess[9*j+n], multi_hess[9*j+n], 1.0e-10)) {
1740 fprintf (stderr, "\nj = %d n = %d \n", j, n);
1741 fprintf (stderr, "norm_hess[j] = %1.14e + %1.14ei\n",
1742 real(norm_hess[9*j+n]), imag(norm_hess[9*j+n]));
1743 fprintf (stderr, "multi_hess[j] = %1.14e + %1.15ei\n",
1744 real(multi_hess[9*j+n]), imag(multi_hess[9*j+n]));
1745 return -9;
1746 }
1747 }
1748 }
1749 return 0;
1750 }
1751
1752
test_complex_double_vgh()1753 void test_complex_double_vgh()
1754 {
1755 int Nx=73; int Ny=91; int Nz = 29;
1756 int num_splines = 128;
1757
1758 Ugrid x_grid, y_grid, z_grid;
1759 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
1760 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
1761 z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
1762
1763 BCtype_z xBC, yBC, zBC;
1764 xBC.lCode = xBC.rCode = PERIODIC;
1765 yBC.lCode = yBC.rCode = PERIODIC;
1766 zBC.lCode = zBC.rCode = PERIODIC;
1767
1768 // First, create splines the normal way
1769 UBspline_3d_z* norm_splines[num_splines];
1770 multi_UBspline_3d_z *multi_spline;
1771
1772 // First, create multispline
1773 multi_spline = create_multi_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC,
1774 num_splines);
1775
1776 complex_double data[Nx*Ny*Nz];
1777 // Now, create normal splines and set multispline data
1778 for (int i=0; i<num_splines; i++) {
1779 for (int j=0; j<Nx*Ny*Nz; j++)
1780 data[j] = complex<double>((drand48()-0.5), (drand48()-0.5));
1781 norm_splines[i] = create_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
1782 set_multi_UBspline_3d_z (multi_spline, i, data);
1783 }
1784
1785 fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
1786 real(norm_splines[19]->coefs[227]),
1787 imag(norm_splines[19]->coefs[227]));
1788 fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
1789 real(multi_spline->coefs[19+227*multi_spline->z_stride]),
1790 imag(multi_spline->coefs[19+227*multi_spline->z_stride]));
1791
1792 // Now, test random values
1793 int num_vals = 100;
1794 complex_double multi_vals[num_splines], norm_vals[num_splines];
1795 complex_double multi_grads[3*num_splines], norm_grads[3*num_splines];
1796 complex_double multi_lapl[num_splines], norm_lapl[num_splines];
1797 complex_double multi_hess[9*num_splines], norm_hess[9*num_splines];
1798 for (int i=0; i<num_vals; i++) {
1799 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1800 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1801 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1802
1803 ///////////////////////
1804 // Check VGH routine //
1805 ///////////////////////
1806 eval_multi_UBspline_3d_z_vgh (multi_spline, x, y, z,
1807 multi_vals, multi_grads, multi_hess);
1808 for (int j=0; j<num_splines; j++)
1809 eval_UBspline_3d_z_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
1810 &(norm_grads[3*j]), &(norm_hess[9*j]));
1811 for (int j=0; j<num_splines; j++) {
1812 // Check value
1813 if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12)) {
1814 fprintf (stderr, "Error! norm_vals[j] = %1.14e + %1.14ei\n",
1815 real(norm_vals[j]), imag(norm_vals[j]));
1816 fprintf (stderr, " multi_vals[j] = %1.14e + %1.14ei\n",
1817 real(multi_vals[j]), imag(multi_vals[j]));
1818 }
1819 // Check gradients
1820 for (int n=0; n<3; n++) {
1821 if (zdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-12)) {
1822 fprintf (stderr, "n=%d\n", n);
1823 fprintf (stderr, "Error! norm_grads[j] = %1.14e + %1.14ei\n",
1824 real(norm_grads[3*j+n]), imag(norm_grads[3*j+n]));
1825 fprintf (stderr, " multi_grads[j] = %1.14e + %1.14ei\n",
1826 real(multi_grads[3*j+n]), imag(multi_grads[3*j+n]));
1827 }
1828 }
1829 // Check hessian
1830 for (int n=0; n<9; n++) {
1831 if (zdiff (norm_hess[9*j+n], multi_hess[9*j+n], 1.0e-10)) {
1832 fprintf (stderr, "Error! norm_hess[j] = %1.14e + %1.14ei\n",
1833 real(norm_hess[9*j+n]), imag(norm_hess[9*j+n]));
1834 fprintf (stderr, " multi_hess[j] = %1.14e + %1.14ei\n",
1835 real(multi_hess[9*j+n]), imag(multi_hess[9*j+n]));
1836 }
1837 }
1838 }
1839 }
1840
1841 num_vals = 100000;
1842
1843 // Now do timing
1844 clock_t norm_start, norm_end, multi_start, multi_end, rand_start, rand_end;
1845 rand_start = clock();
1846 for (int i=0; i<num_vals; i++) {
1847 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1848 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1849 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1850 }
1851 rand_end = clock();
1852
1853 norm_start = clock();
1854 for (int i=0; i<num_vals; i++) {
1855 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1856 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1857 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1858
1859 for (int j=0; j<num_splines; j++)
1860 eval_UBspline_3d_z_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
1861 &(norm_grads[3*j]), &(norm_hess[9*j]));
1862 }
1863 norm_end = clock();
1864
1865 multi_start = clock();
1866 for (int i=0; i<num_vals; i++) {
1867 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1868 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1869 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1870 eval_multi_UBspline_3d_z_vgh (multi_spline, x, y, z, multi_vals, multi_grads, multi_hess);
1871 }
1872 multi_end = clock();
1873
1874 fprintf (stderr, "Normal spline time = %1.5f\n",
1875 (double)(norm_end-norm_start+rand_start-rand_end)/CLOCKS_PER_SEC);
1876 fprintf (stderr, "Multi spline time = %1.5f\n",
1877 (double)(multi_end-multi_start+rand_start-rand_end)/CLOCKS_PER_SEC);
1878
1879 }
1880
1881
test_double()1882 void test_double()
1883 {
1884 int Nx=73; int Ny=91; int Nz = 29;
1885 int num_splines = 201;
1886
1887 Ugrid x_grid, y_grid, z_grid;
1888 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
1889 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
1890 z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
1891
1892 BCtype_d xBC, yBC, zBC;
1893 xBC.lCode = xBC.rCode = PERIODIC;
1894 yBC.lCode = yBC.rCode = PERIODIC;
1895 zBC.lCode = zBC.rCode = PERIODIC;
1896
1897 // First, create splines the normal way
1898 UBspline_3d_d* norm_splines[num_splines];
1899 multi_UBspline_3d_d *multi_spline;
1900
1901 // First, create multispline
1902 multi_spline = create_multi_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC,
1903 num_splines);
1904
1905 double data[Nx*Ny*Nz];
1906 // Now, create normal splines and set multispline data
1907 for (int i=0; i<num_splines; i++) {
1908 for (int j=0; j<Nx*Ny*Nz; j++)
1909 data[j] = (drand48()-0.5);
1910 norm_splines[i] = create_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
1911 set_multi_UBspline_3d_d (multi_spline, i, data);
1912 }
1913
1914 fprintf (stderr, "norm coef = %1.14e\n",
1915 norm_splines[19]->coefs[227]);
1916 fprintf (stderr, "multi coef = %1.14e\n",
1917 multi_spline->coefs[19+227*multi_spline->z_stride]);
1918
1919 // Now, test random values
1920 int num_vals = 100;
1921 double multi_vals[num_splines], norm_vals[num_splines];
1922 for (int i=0; i<num_vals; i++) {
1923 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1924 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1925 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1926
1927 eval_multi_UBspline_3d_d (multi_spline, x, y, z,
1928 multi_vals);
1929 for (int j=0; j<num_splines; j++)
1930 eval_UBspline_3d_d (norm_splines[j], x, y, z, &(norm_vals[j]));
1931 for (int j=0; j<num_splines; j++) {
1932 // Check value
1933 double diff = norm_vals[j] - multi_vals[j];
1934 if (fabs(diff) > 1.0e-12) {
1935 fprintf (stderr, "Error! norm_vals[j] = %1.14e\n",
1936 norm_vals[j]);
1937 fprintf (stderr, " multi_vals[j] = %1.14e\n",
1938 multi_vals[j]);
1939 }
1940 }
1941 }
1942
1943 num_vals = 100000;
1944
1945 // Now do timing
1946 clock_t norm_start, norm_end, multi_start, multi_end, rand_start, rand_end;
1947 rand_start = clock();
1948 for (int i=0; i<num_vals; i++) {
1949 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1950 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1951 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1952 }
1953 rand_end = clock();
1954
1955 norm_start = clock();
1956 for (int i=0; i<num_vals; i++) {
1957 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1958 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1959 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1960
1961 for (int j=0; j<num_splines; j++)
1962 eval_UBspline_3d_d (norm_splines[j], x, y, z, &(norm_vals[j]));
1963 }
1964 norm_end = clock();
1965
1966 multi_start = clock();
1967 for (int i=0; i<num_vals; i++) {
1968 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
1969 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
1970 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
1971 eval_multi_UBspline_3d_d (multi_spline, x, y, z, multi_vals);
1972 }
1973 multi_end = clock();
1974
1975 fprintf (stderr, "Normal spline time = %1.5f\n",
1976 (double)(norm_end-norm_start+rand_start-rand_end)/CLOCKS_PER_SEC);
1977 fprintf (stderr, "Multi spline time = %1.5f\n",
1978 (double)(multi_end-multi_start+rand_start-rand_end)/CLOCKS_PER_SEC);
1979
1980 }
1981
1982
1983
test_double_vgh()1984 void test_double_vgh()
1985 {
1986 int Nx=73; int Ny=91; int Nz = 29;
1987 int num_splines = 128;
1988
1989 Ugrid x_grid, y_grid, z_grid;
1990 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
1991 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
1992 z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
1993
1994 BCtype_d xBC, yBC, zBC;
1995 xBC.lCode = xBC.rCode = PERIODIC;
1996 yBC.lCode = yBC.rCode = PERIODIC;
1997 zBC.lCode = zBC.rCode = PERIODIC;
1998
1999 // First, create splines the normal way
2000 UBspline_3d_d* norm_splines[num_splines];
2001 multi_UBspline_3d_d *multi_spline;
2002
2003 // First, create multispline
2004 multi_spline = create_multi_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC,
2005 num_splines);
2006
2007 double data[Nx*Ny*Nz];
2008 // Now, create normal splines and set multispline data
2009 for (int i=0; i<num_splines; i++) {
2010 for (int j=0; j<Nx*Ny*Nz; j++)
2011 data[j] = (drand48()-0.5);
2012 norm_splines[i] = create_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
2013 set_multi_UBspline_3d_d (multi_spline, i, data);
2014 }
2015
2016 fprintf (stderr, "norm coef = %1.14e\n",
2017 norm_splines[19]->coefs[227]);
2018 fprintf (stderr, "multi coef = %1.14e\n",
2019 multi_spline->coefs[19+227*multi_spline->z_stride]);
2020
2021 // Now, test random values
2022 int num_vals = 100;
2023 double multi_vals[num_splines], norm_vals[num_splines];
2024 double multi_grads[3*num_splines], norm_grads[3*num_splines];
2025 double multi_hess[9*num_splines], norm_hess[9*num_splines];
2026 for (int i=0; i<num_vals; i++) {
2027 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
2028 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
2029 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
2030
2031 eval_multi_UBspline_3d_d_vgh (multi_spline, x, y, z,
2032 multi_vals, multi_grads, multi_hess);
2033 for (int j=0; j<num_splines; j++)
2034 eval_UBspline_3d_d_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
2035 &(norm_grads[3*j]), &(norm_hess[9*j]));
2036 for (int j=0; j<num_splines; j++) {
2037 // Check value
2038 double diff = norm_vals[j] - multi_vals[j];
2039 if (fabs(diff) > 1.0e-12) {
2040 fprintf (stderr, "j = %d\n", j);
2041 fprintf (stderr, "Error! norm_vals[j] = %1.14e\n",
2042 norm_vals[j]);
2043 fprintf (stderr, " multi_vals[j] = %1.14e\n",
2044 multi_vals[j]);
2045 }
2046 // Check gradients
2047 for (int n=0; n<3; n++) {
2048 diff = norm_grads[3*j+n] - multi_grads[3*j+n];
2049 if (fabs(diff) > 1.0e-12) {
2050 fprintf (stderr, "n=%d\n", n);
2051 fprintf (stderr, "Error! norm_grads[j] = %1.14e\n",
2052 norm_grads[3*j+n]);
2053 fprintf (stderr, " multi_grads[j] = %1.14e\n",
2054 multi_grads[3*j+n]);
2055 }
2056 }
2057 // Check hessian
2058 for (int n=0; n<9; n++) {
2059 diff = norm_hess[9*j+n] - multi_hess[9*j+n];
2060 if (fabs(diff) > 1.0e-10) {
2061 fprintf (stderr, "Error! norm_hess[j] = %1.14e\n",
2062 norm_hess[9*j+n]);
2063 fprintf (stderr, " multi_hess[j] = %1.14e\n",
2064 multi_hess[9*j+n]);
2065 }
2066 }
2067 }
2068 }
2069
2070 num_vals = 100000;
2071
2072 // Now do timing
2073 clock_t norm_start, norm_end, multi_start, multi_end, rand_start, rand_end;
2074 rand_start = clock();
2075 for (int i=0; i<num_vals; i++) {
2076 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
2077 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
2078 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
2079 }
2080 rand_end = clock();
2081
2082 norm_start = clock();
2083 for (int i=0; i<num_vals; i++) {
2084 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
2085 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
2086 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
2087
2088 for (int j=0; j<num_splines; j++)
2089 eval_UBspline_3d_d_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
2090 &(norm_grads[3*j]), &(norm_hess[9*j]));
2091 }
2092 norm_end = clock();
2093
2094 multi_start = clock();
2095 for (int i=0; i<num_vals; i++) {
2096 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
2097 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
2098 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
2099 eval_multi_UBspline_3d_d_vgh (multi_spline, x, y, z, multi_vals, multi_grads, multi_hess);
2100 }
2101 multi_end = clock();
2102
2103 fprintf (stderr, "Normal spline time = %1.5f\n",
2104 (double)(norm_end-norm_start+rand_start-rand_end)/CLOCKS_PER_SEC);
2105 fprintf (stderr, "Multi spline time = %1.5f\n",
2106 (double)(multi_end-multi_start+rand_start-rand_end)/CLOCKS_PER_SEC);
2107
2108 }
2109
PrintPassFail(int code)2110 void PrintPassFail (int code)
2111 {
2112 char green[100], normal[100], red[100];
2113 snprintf (green, 100, "%c[0;32;47m", 0x1b);
2114 snprintf (normal, 100, "%c[0;30;47m", 0x1b);
2115 snprintf (red, 100, "%c[0;31;47m", 0x1b);
2116
2117 if (code == 0)
2118 fprintf (stderr, "PASSED\n");
2119 else
2120 fprintf (stderr, "FAILED: code = %d\n", code);
2121 }
2122
2123
main()2124 main()
2125 {
2126 int code;
2127 //test_complex_double();
2128 //test_complex_double_vgh();
2129
2130 fprintf (stderr, "Testing 1D real single-precision multiple cubic B-spline routines: ");
2131 code = test_1d_float_all(); PrintPassFail (code);
2132 fprintf (stderr, "Testing 2D real single-precision multiple cubic B-spline routines: ");
2133 code = test_2d_float_all(); PrintPassFail (code);
2134 fprintf (stderr, "Testing 3D real single-precision multiple cubic B-spline routines: ");
2135 code = test_3d_float_all(); PrintPassFail (code);
2136
2137 fprintf (stderr, "Testing 1D real double-precision multiple cubic B-spline routines: ");
2138 code = test_1d_double_all(); PrintPassFail (code);
2139 fprintf (stderr, "Testing 2D real double-precision multiple cubic B-spline routines: ");
2140 code = test_2d_double_all(); PrintPassFail (code);
2141 fprintf (stderr, "Testing 3D real double-precision multiple cubic B-spline routines: ");
2142 code = test_3d_double_all(); PrintPassFail (code);
2143
2144 fprintf (stderr, "Testing 1D complex single-precision multiple cubic B-spline routines: ");
2145 code = test_1d_complex_float_all(); PrintPassFail (code);
2146 fprintf (stderr, "Testing 2D complex single-precision multiple cubic B-spline routines: ");
2147 code = test_2d_complex_float_all(); PrintPassFail (code);
2148 fprintf (stderr, "Testing 3D complex single-precision multiple cubic B-spline routines: ");
2149 code = test_3d_complex_float_all(); PrintPassFail (code);
2150
2151 fprintf (stderr, "Testing 1D complex double-precision multiple cubic B-spline routines: ");
2152 code = test_1d_complex_double_all(); PrintPassFail (code);
2153 fprintf (stderr, "Testing 2D complex double-precision multiple cubic B-spline routines: ");
2154 code = test_2d_complex_double_all(); PrintPassFail (code);
2155 fprintf (stderr, "Testing 3D complex double-precision multiple cubic B-spline routines: ");
2156 code = test_3d_complex_double_all(); PrintPassFail (code);
2157 //test_double();
2158 //test_double_vgh();
2159 }
2160