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 "nubspline.h"
8 #include <stdio.h>
9 #include <assert.h>
10 #include <stdlib.h>
11 #include <time.h>
12 #include <math.h>
13 #include <string.h>
14
15 #ifndef M_PI
16 #define M_PI 3.1415926535897932384626433
17 #endif
18
19 double drand48();
20
21 void
PrintPassFail(bool pass)22 PrintPassFail(bool pass)
23 {
24 if (pass)
25 // Print green "Passed"
26 fprintf (stderr, "%c[32mPassed%c[0m\n", 0x1B, 0x1B);
27 else
28 // Print red "Failed"
29 fprintf (stderr, "%c[31mFailed%c[0m\n", 0x1B, 0x1B);
30 }
31
PrintTest(char * name,bool pass)32 void PrintTest (char *name, bool pass)
33 {
34 int n = strlen (name);
35 fprintf (stderr, "%s:", name);
36 for (int i=n; i<57; i++)
37 fprintf (stderr, " ");
38 PrintPassFail (pass);
39 }
40
41
42 bool
TestCenterGrid()43 TestCenterGrid()
44 {
45 fprintf (stderr, "Testing CenterGrid: ");
46 bool passed = true;
47 NUgrid* grid = create_center_grid (-5.0, 7.0, 6.0, 200);
48
49 for (int i=0; i<10000; i++) {
50 double x = -5.0+12.0*drand48();
51 int lo = (*grid->reverse_map)(grid, x);
52 assert (x >= grid->points[lo]);
53 assert (x <= grid->points[lo+1]);
54 }
55 PrintPassFail (passed);
56 return passed;
57 }
58
59
60 bool
TestGeneralGrid()61 TestGeneralGrid()
62 {
63 fprintf (stderr, "Testing GeneralGrid: ");
64 bool passed = true;
65 NUgrid* centgrid = create_center_grid (-5.0, 7.0, 6.0, 200);
66 NUgrid* grid = create_general_grid (centgrid->points, 200);
67 for (int i=0; i<10000; i++) {
68 double x = -5.0+12.0*drand48();
69 int lo = (*grid->reverse_map)(grid, x);
70 passed = passed && (x >= grid->points[lo]);
71 passed = passed && (x <= grid->points[lo+1]);
72 }
73 PrintPassFail (passed);
74 return passed;
75 }
76
77 bool
close_float(float x,float y)78 close_float (float x, float y)
79 {
80 float max = fmaxf (x, y);
81 return (fabs(x-y)/max < 1.0e-5);
82 }
83
84 bool
TestNUB_1d_s()85 TestNUB_1d_s()
86 {
87 double start = -5.0;
88 double end = 7.0;
89 int N = 200;
90 NUgrid* grid = create_center_grid (start, end, 6.0, N);
91 bool passed = true;
92 float data[N];
93 for (int i=0; i<N; i++)
94 data[i] = -1.0 + 2.0*drand48();
95 BCtype_s bc;
96
97 // Create spline with PBC
98 fprintf (stderr, "Testing 1D single-precision periodic boundary conditions:\n");
99 bc.lCode = PERIODIC; bc.rCode = PERIODIC;
100 NUBspline_1d_s *periodic = create_NUBspline_1d_s (grid, bc, data);
101 float sval, sgrad, slapl, eval, egrad, elapl;
102 eval_NUBspline_1d_s_vgl (periodic, start, &sval, &sgrad, &slapl);
103 eval_NUBspline_1d_s_vgl (periodic, end , &eval, &egrad, &elapl);
104 bool v_passed, grad_passed, lapl_passed;
105 v_passed = close_float (sval, eval);
106 grad_passed = close_float (sgrad, egrad);
107 lapl_passed = close_float (slapl, elapl);
108 PrintTest ("Value", v_passed);
109 PrintTest ("First derivative", grad_passed);
110 PrintTest ("Second derivative", lapl_passed);
111 passed = passed && v_passed && grad_passed && lapl_passed;
112
113 double x = grid->points[26];
114 float val;
115 eval_NUBspline_1d_s (periodic, x, &val);
116 bool interp_passed = close_float (val, data[26]);
117 PrintTest ("Interpolation", interp_passed);
118 passed = passed && interp_passed;
119
120 // Create spline with fixed first derivative:
121 bc.lCode = DERIV1; bc.lVal = 1.5;
122 bc.rCode = DERIV1; bc.rVal = -0.3;
123 NUBspline_1d_s *fixed_first = create_NUBspline_1d_s (grid, bc, data);
124 fprintf (stderr, "Testing 1D single-precsion fixed first derivative boundary conditions: \n");
125 eval_NUBspline_1d_s_vg (fixed_first, start, &sval, &sgrad);
126 eval_NUBspline_1d_s_vg (fixed_first, end, &eval, &egrad);
127 bool bc_passed = close_float (sgrad, 1.5) && close_float (egrad, -0.3);
128 PrintTest ("Boundary conditions", bc_passed);
129 x = grid->points[26];
130 eval_NUBspline_1d_s (periodic, x, &val);
131 interp_passed = close_float (val, data[26]);
132 PrintTest ("Interpolation", interp_passed);
133 passed = passed && interp_passed && bc_passed;
134
135 // Create spline with fixed second derivative:
136 bc.lCode = DERIV2; bc.lVal = 1.5;
137 bc.rCode = DERIV2; bc.rVal = -0.3;
138 NUBspline_1d_s *fixed_second = create_NUBspline_1d_s (grid, bc, data);
139 fprintf (stderr, "Testing 1d_s fixed second derivative boundary conditions: \n");
140 eval_NUBspline_1d_s_vgl (fixed_second, start, &sval, &sgrad, &slapl);
141 eval_NUBspline_1d_s_vgl (fixed_second, end, &eval, &egrad, &elapl);
142 bc_passed = close_float (slapl, 1.5) && close_float (elapl, -0.3);
143 fprintf (stderr, "slapl = %1.8f elapl = %1.8f\n", slapl, elapl);
144 PrintTest ("Boundary conditions", bc_passed);
145 x = grid->points[26];
146 eval_NUBspline_1d_s (periodic, x, &val);
147 interp_passed = close_float (val, data[26]);
148 PrintTest ("Interpolation", interp_passed);
149 passed = passed && interp_passed && bc_passed;
150
151 return passed;
152 }
153
154 void
GridSpeedTest()155 GridSpeedTest()
156 {
157 NUgrid* centgrid = create_center_grid (-5.0, 7.0, 6.0, 2000);
158 NUgrid* gengrid = create_general_grid (centgrid->points, 2000);
159 int centsum=0, gensum=0;
160
161 clock_t rstart, rend, cstart, cend, gstart, gend;
162
163 rstart = clock();
164 for (int i=0; i<100000000; i++) {
165 double x = -5.0 + 12.0*drand48();
166 }
167 rend = clock();
168
169 cstart = clock();
170 for (int i=0; i<100000000; i++) {
171 double x = -5.0 + 12.0*drand48();
172 centsum += (*centgrid->reverse_map)(centgrid, x);
173 }
174 cend = clock();
175
176 gstart = clock();
177 for (int i=0; i<100000000; i++) {
178 double x = -5.0 + 12.0*drand48();
179 gensum += (*gengrid->reverse_map)(gengrid, x);
180 }
181 gend = clock();
182
183 double cent_time = (double)(cend-cstart+rstart-rend)/(double)CLOCKS_PER_SEC;
184 double gen_time = (double)(gend-gstart+rstart-rend)/(double)CLOCKS_PER_SEC;
185 fprintf (stderr, "%d %d\n", centsum, gensum);
186 fprintf (stderr, "center_grid time = %1.3f s.\n", cent_time);
187 fprintf (stderr, "general_grid time = %1.3f s.\n", gen_time);
188 }
189
190 void
TestNUBasis()191 TestNUBasis()
192 {
193 NUgrid* centgrid = create_center_grid (-5.0, 7.0, 10.0, 20);
194 NUBasis* basis = create_NUBasis (centgrid, true);
195
196 double bfuncs[4];
197 for (double x=-5.0; x<=7.0; x+=0.001) {
198 get_NUBasis_funcs_d (basis, x, bfuncs);
199 fprintf (stderr, "%1.12f %1.12f %1.12f %1.12f %1.12f\n",
200 x, bfuncs[0], bfuncs[1], bfuncs[2], bfuncs[3]);
201 }
202 }
203
204 void
TestNUBspline()205 TestNUBspline()
206 {
207 NUgrid* centgrid = create_center_grid (-5.0, 7.0, 10.0, 20);
208 NUBasis* basis = create_NUBasis (centgrid, true);
209 float data[20];
210 for (int i=0; i<20; i++) {
211 double x = centgrid->points[i];
212 double angle = (x+5.0)/12.0 * 2.0*M_PI;
213 data[i] = sin(angle);
214 }
215 BCtype_s bc;
216 // bc.lCode = PERIODIC; bc.rCode = PERIODIC;
217 bc.lCode = DERIV1; bc.lVal = 2.0*M_PI/12.0;
218 bc.rCode = DERIV1; bc.rVal = 2.0*M_PI/12.0;
219 //bc.lCode = NATURAL; bc.rCode = FLAT;
220 NUBspline_1d_s *spline = create_NUBspline_1d_s (centgrid, bc, data);
221 for (double x=-5.0; x<=7.0; x+=0.001) {
222 float val, deriv;
223 eval_NUBspline_1d_s_vg (spline, x, &val, &deriv);
224 double angle = (x+5.0)/12.0 * 2.0*M_PI;
225 fprintf (stderr, "%1.16e %1.16e %1.16e %1.16e\n", x, val,
226 sin(angle), deriv);
227 }
228 }
229
230
231 void
TestNUBspline_d()232 TestNUBspline_d()
233 {
234 NUgrid* centgrid = create_center_grid (-5.0, 7.0, 10.0, 20);
235 NUBasis* basis = create_NUBasis (centgrid, true);
236 double data[20];
237 for (int i=0; i<20; i++) {
238 double x = centgrid->points[i];
239 double angle = (x+5.0)/12.0 * 2.0*M_PI;
240 data[i] = sin(angle);
241 }
242 BCtype_d bc;
243 // bc.lCode = PERIODIC; bc.rCode = PERIODIC;
244 bc.lCode = DERIV1; bc.lVal = 2.0*M_PI/12.0;
245 bc.rCode = DERIV1; bc.rVal = 2.0*M_PI/12.0;
246 //bc.lCode = NATURAL; bc.rCode = FLAT;
247 NUBspline_1d_d *spline = create_NUBspline_1d_d (centgrid, bc, data);
248 for (double x=-5.0; x<=7.0; x+=0.001) {
249 double val, deriv;
250 eval_NUBspline_1d_d_vg (spline, x, &val, &deriv);
251 double angle = (x+5.0)/12.0 * 2.0*M_PI;
252 fprintf (stderr, "%1.16e %1.16e %1.16e %1.16e\n", x, val,
253 sin(angle), deriv);
254 }
255 }
256
257
258 void
TestNUB_2d_s()259 TestNUB_2d_s()
260 {
261 int Mx=30, My=35;
262 NUgrid *x_grid = create_center_grid (-3.0, 4.0, 7.5, Mx);
263 NUgrid *y_grid = create_center_grid (-1.0, 9.0, 3.5, My);
264 float data[Mx*My];
265 for (int ix=0; ix<Mx; ix++)
266 for (int iy=0; iy<My; iy++)
267 data[ix*My+iy] = -1.0+2.0*drand48();
268
269 BCtype_s xBC, yBC;
270 xBC.lCode = PERIODIC;
271 yBC.lCode = PERIODIC;
272 // xBC.lCode = FLAT; xBC.rCode = FLAT;
273 // yBC.lCode = FLAT; yBC.rCode = FLAT;
274
275 NUBspline_2d_s *spline = create_NUBspline_2d_s (x_grid, y_grid, xBC, yBC, data);
276
277 int xFine = 400;
278 int yFine = 400;
279 FILE *fout = fopen ("2d_s.dat", "w");
280 double xi = x_grid->start;
281 double xf = x_grid->end;// + x_grid->points[1] - x_grid->points[0];
282 double yi = y_grid->start;
283 double yf = y_grid->end;// + y_grid->points[1] - y_grid->points[0];
284 for (int ix=0; ix<xFine; ix++) {
285 double x = xi+ (double)ix/(double)(xFine)*(xf-xi);
286 for (int iy=0; iy<yFine; iy++) {
287 double y = yi + (double)iy/(double)(yFine)*(yf-yi);
288 float val;
289 eval_NUBspline_2d_s (spline, x, y, &val);
290 fprintf (fout, "%1.16e ", val);
291 }
292 fprintf (fout, "\n");
293 }
294 fclose (fout);
295 }
296
297
298 void
TestNUB_2d_c()299 TestNUB_2d_c()
300 {
301 int Mx=30, My=35;
302 NUgrid *x_grid = create_center_grid (-3.0, 4.0, 7.5, Mx);
303 NUgrid *y_grid = create_center_grid (-1.0, 9.0, 3.5, My);
304 complex_float data[Mx*My];
305 for (int ix=0; ix<Mx; ix++)
306 for (int iy=0; iy<My; iy++)
307 data[ix*My+iy] = -1.0+2.0*drand48() + 1.0fi*(-1.0+2.0*drand48());
308
309 BCtype_c xBC, yBC;
310 xBC.lCode = PERIODIC;
311 yBC.lCode = PERIODIC;
312 // xBC.lCode = FLAT; xBC.rCode = FLAT;
313 // yBC.lCode = FLAT; yBC.rCode = FLAT;
314
315 NUBspline_2d_c *spline = create_NUBspline_2d_c (x_grid, y_grid, xBC, yBC, data);
316
317 int xFine = 400;
318 int yFine = 400;
319 FILE *rout = fopen ("2d_r.dat", "w");
320 FILE *iout = fopen ("2d_i.dat", "w");
321 double xi = x_grid->start;
322 double xf = x_grid->end;// + x_grid->points[1] - x_grid->points[0];
323 double yi = y_grid->start;
324 double yf = y_grid->end;// + y_grid->points[1] - y_grid->points[0];
325 for (int ix=0; ix<xFine; ix++) {
326 double x = xi+ (double)ix/(double)(xFine)*(xf-xi);
327 for (int iy=0; iy<yFine; iy++) {
328 double y = yi + (double)iy/(double)(yFine)*(yf-yi);
329 complex_float val, grad[2], hess[4];
330 eval_NUBspline_2d_c_vgh (spline, x, y, &val, grad, hess);
331 fprintf (rout, "%1.16e ", crealf(val));
332 fprintf (iout, "%1.16e ", cimagf(val));
333 }
334 fprintf (rout, "\n");
335 fprintf (iout, "\n");
336 }
337 fclose (rout);
338 fclose (iout);
339 }
340
341 void
TestNUB_3d_s()342 TestNUB_3d_s()
343 {
344 int Mx=20, My=27, Mz=23;
345 NUgrid *x_grid = create_center_grid (-3.0, 4.0, 7.5, Mx);
346 NUgrid *y_grid = create_center_grid (-1.0, 9.0, 3.5, My);
347 NUgrid *z_grid = create_center_grid (-1.8, 2.0, 2.8, Mz);
348 float data[Mx*My*Mz];
349 for (int ix=0; ix<Mx; ix++)
350 for (int iy=0; iy<My; iy++)
351 for (int iz=0; iz<Mz; iz++)
352 data[(ix*My+iy)*Mz+iz] = -1.0+2.0*drand48();
353
354 BCtype_s xBC, yBC, zBC;
355 // xBC.lCode = PERIODIC;
356 // yBC.lCode = PERIODIC;
357 xBC.lCode = PERIODIC; xBC.rCode = PERIODIC;
358 yBC.lCode = PERIODIC; yBC.rCode = PERIODIC;
359 zBC.lCode = PERIODIC; zBC.rCode = PERIODIC;
360
361 NUBspline_3d_s *spline = create_NUBspline_3d_s (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
362
363 int xFine = 200, yFine = 200, zFine=200;
364 FILE *fout = fopen ("3d_s.dat", "w");
365 double xi = x_grid->start; double xf = x_grid->end;
366 double yi = y_grid->start; double yf = y_grid->end;
367 double zi = z_grid->start; double zf = z_grid->end;
368 for (int ix=0; ix<xFine; ix++) {
369 double x = xi+ (double)ix/(double)(xFine)*(xf-xi);
370 for (int iy=0; iy<yFine; iy++) {
371 double y = yi + (double)iy/(double)(yFine)*(yf-yi);
372 for (int iz=0; iz<zFine; iz++) {
373 double z = zi + (double)iz/(double)(zFine)*(zf-zi);
374 float val, grad[3], hess[9];
375 eval_NUBspline_3d_s_vgh (spline, x, y, z, &val, grad, hess);
376 fprintf (fout, "%1.16e ", val);
377 }
378 }
379 fprintf (fout, "\n");
380 }
381 fclose (fout);
382 fprintf (stderr, "spline->sp_code = %d\n", spline->sp_code);
383 destroy_Bspline (spline);
384 }
385
386
387 void
TestNUB_3d_d()388 TestNUB_3d_d()
389 {
390 int Mx=20, My=27, Mz=23;
391 NUgrid *x_grid = create_center_grid (-3.0, 4.0, 7.5, Mx);
392 NUgrid *y_grid = create_center_grid (-1.0, 9.0, 3.5, My);
393 NUgrid *z_grid = create_center_grid (-1.8, 2.0, 2.8, Mz);
394 double data[Mx*My*Mz];
395 for (int ix=0; ix<Mx; ix++)
396 for (int iy=0; iy<My; iy++)
397 for (int iz=0; iz<Mz; iz++)
398 data[(ix*My+iy)*Mz+iz] = -1.0+2.0*drand48();
399
400 BCtype_d xBC, yBC, zBC;
401 // xBC.lCode = PERIODIC;
402 // yBC.lCode = PERIODIC;
403 xBC.lCode = PERIODIC; xBC.rCode = PERIODIC;
404 yBC.lCode = PERIODIC; yBC.rCode = PERIODIC;
405 zBC.lCode = PERIODIC; zBC.rCode = PERIODIC;
406
407 NUBspline_3d_d *spline = create_NUBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
408
409 int xFine = 200, yFine = 200, zFine=200;
410 FILE *fout = fopen ("3d_d.dat", "w");
411 double xi = x_grid->start; double xf = x_grid->end;
412 double yi = y_grid->start; double yf = y_grid->end;
413 double zi = z_grid->start; double zf = z_grid->end;
414 for (int ix=0; ix<xFine; ix++) {
415 double x = xi+ (double)ix/(double)(xFine)*(xf-xi);
416 for (int iy=0; iy<yFine; iy++) {
417 double y = yi + (double)iy/(double)(yFine)*(yf-yi);
418 for (int iz=0; iz<zFine; iz++) {
419 double z = zi + (double)iz/(double)(zFine)*(zf-zi);
420 double val, grad[3], hess[9];
421 eval_NUBspline_3d_d_vgh (spline, x, y, z, &val, grad, hess);
422 fprintf (fout, "%1.16e ", val);
423 }
424 }
425 fprintf (fout, "\n");
426 }
427 fclose (fout);
428 fprintf (stderr, "spline->sp_code = %d\n", spline->sp_code);
429 destroy_Bspline (spline);
430 }
431
432 void
TestNUB_3d_c()433 TestNUB_3d_c()
434 {
435 int Mx=20, My=27, Mz=23;
436 NUgrid *x_grid = create_center_grid (-3.0, 4.0, 7.5, Mx);
437 NUgrid *y_grid = create_center_grid (-1.0, 9.0, 3.5, My);
438 NUgrid *z_grid = create_center_grid (-1.8, 2.0, 2.8, Mz);
439 complex_float data[Mx*My*Mz];
440 for (int ix=0; ix<Mx; ix++)
441 for (int iy=0; iy<My; iy++)
442 for (int iz=0; iz<Mz; iz++)
443 data[(ix*My+iy)*Mz+iz] = -1.0+2.0*drand48() + 1.0if*(-1.0+2.0*drand48());
444
445 BCtype_c xBC, yBC, zBC;
446 // xBC.lCode = PERIODIC;
447 // yBC.lCode = PERIODIC;
448 xBC.lCode = PERIODIC; xBC.rCode = PERIODIC;
449 yBC.lCode = PERIODIC; yBC.rCode = PERIODIC;
450 zBC.lCode = PERIODIC; zBC.rCode = PERIODIC;
451
452 NUBspline_3d_c *spline = create_NUBspline_3d_c (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
453
454 int xFine = 200, yFine = 200, zFine=200;
455 FILE *rout = fopen ("3d_r.dat", "w");
456 FILE *iout = fopen ("3d_i.dat", "w");
457 double xi = x_grid->start; double xf = x_grid->end;
458 double yi = y_grid->start; double yf = y_grid->end;
459 double zi = z_grid->start; double zf = z_grid->end;
460 for (int ix=0; ix<xFine; ix++) {
461 double x = xi+ (double)ix/(double)(xFine)*(xf-xi);
462 for (int iy=0; iy<yFine; iy++) {
463 double y = yi + (double)iy/(double)(yFine)*(yf-yi);
464 for (int iz=0; iz<zFine; iz++) {
465 double z = zi + (double)iz/(double)(zFine)*(zf-zi);
466 complex_float val, grad[3], hess[9];
467 eval_NUBspline_3d_c_vgh (spline, x, y, z, &val, grad, hess);
468 fprintf (rout, "%1.16e ", crealf(val));
469 fprintf (iout, "%1.16e ", cimagf(val));
470 }
471 }
472 fprintf (rout, "\n");
473 fprintf (iout, "\n");
474 }
475 fclose (rout);
476 fclose (iout);
477 }
478
479
480 void
TestNUB_3d_z()481 TestNUB_3d_z()
482 {
483 int Mx=20, My=27, Mz=23;
484 NUgrid *x_grid = create_center_grid (-3.0, 4.0, 7.5, Mx);
485 NUgrid *y_grid = create_center_grid (-1.0, 9.0, 3.5, My);
486 NUgrid *z_grid = create_center_grid (-1.8, 2.0, 2.8, Mz);
487 complex_double data[Mx*My*Mz];
488 for (int ix=0; ix<Mx; ix++)
489 for (int iy=0; iy<My; iy++)
490 for (int iz=0; iz<Mz; iz++)
491 data[(ix*My+iy)*Mz+iz] = -1.0+2.0*drand48() + 1.0if*(-1.0+2.0*drand48());
492
493 BCtype_z xBC, yBC, zBC;
494 // xBC.lCode = PERIODIC;
495 // yBC.lCode = PERIODIC;
496 xBC.lCode = PERIODIC; xBC.rCode = PERIODIC;
497 yBC.lCode = PERIODIC; yBC.rCode = PERIODIC;
498 zBC.lCode = PERIODIC; zBC.rCode = PERIODIC;
499
500 NUBspline_3d_z *spline = create_NUBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
501
502 int xFine = 200, yFine = 200, zFine=200;
503 FILE *rout = fopen ("3d_r.dat", "w");
504 FILE *iout = fopen ("3d_i.dat", "w");
505 double xi = x_grid->start; double xf = x_grid->end;
506 double yi = y_grid->start; double yf = y_grid->end;
507 double zi = z_grid->start; double zf = z_grid->end;
508 for (int ix=0; ix<xFine; ix++) {
509 double x = xi+ (double)ix/(double)(xFine)*(xf-xi);
510 for (int iy=0; iy<yFine; iy++) {
511 double y = yi + (double)iy/(double)(yFine)*(yf-yi);
512 for (int iz=0; iz<zFine; iz++) {
513 double z = zi + (double)iz/(double)(zFine)*(zf-zi);
514 complex_double val, grad[3], hess[9];
515 eval_NUBspline_3d_z_vgh (spline, x, y, z, &val, grad, hess);
516 fprintf (rout, "%1.16e ", crealf(val));
517 fprintf (iout, "%1.16e ", cimagf(val));
518 }
519 }
520 fprintf (rout, "\n");
521 fprintf (iout, "\n");
522 }
523 fclose (rout);
524 fclose (iout);
525 }
526
527 void
SpeedNUB_3d_s()528 SpeedNUB_3d_s()
529 {
530 int Mx=200, My=200, Mz=200;
531 NUgrid *x_grid = create_center_grid (-3.0, 4.0, 1.0001, Mx);
532 NUgrid *y_grid = create_center_grid (-1.0, 9.0, 1.0001, My);
533 NUgrid *z_grid = create_center_grid (-1.8, 2.0, 1.0001, Mz);
534 float *data;
535 data = malloc (sizeof(float)*Mx*My*Mz);
536 for (int ix=0; ix<Mx; ix++)
537 for (int iy=0; iy<My; iy++)
538 for (int iz=0; iz<Mz; iz++)
539 data[(ix*My+iy)*Mz+iz] = -1.0+2.0*drand48();
540
541 BCtype_s xBC, yBC, zBC;
542 // xBC.lCode = PERIODIC;
543 // yBC.lCode = PERIODIC;
544 xBC.lCode = PERIODIC; xBC.rCode = PERIODIC;
545 yBC.lCode = PERIODIC; yBC.rCode = PERIODIC;
546 zBC.lCode = PERIODIC; zBC.rCode = PERIODIC;
547
548 NUBspline_3d_s *spline = create_NUBspline_3d_s (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
549
550 float val, grad[3], hess[9];
551 clock_t start, end, rstart, rend;
552 rstart = clock();
553 for (int i=0; i<10000000; i++) {
554 double x = x_grid->start+ 0.9999*drand48()*(x_grid->end - x_grid->start);
555 double y = y_grid->start+ 0.9999*drand48()*(y_grid->end - y_grid->start);
556 double z = z_grid->start+ 0.9999*drand48()*(z_grid->end - z_grid->start);
557 }
558 rend = clock();
559 start = clock();
560 for (int i=0; i<10000000; i++) {
561 double x = x_grid->start+ 0.9999*drand48()*(x_grid->end - x_grid->start);
562 double y = y_grid->start+ 0.9999*drand48()*(y_grid->end - y_grid->start);
563 double z = z_grid->start+ 0.9999*drand48()*(z_grid->end - z_grid->start);
564 eval_NUBspline_3d_s_vgh (spline, x, y, z, &val, grad, hess);
565 }
566 end = clock();
567 fprintf (stderr, "10,000,000 evalations in %f seconds.\n",
568 (double)(end-start-(rend-rstart))/(double)CLOCKS_PER_SEC);
569 }
570
571
572 void
SpeedNUB_3d_z()573 SpeedNUB_3d_z()
574 {
575 int Mx=200, My=200, Mz=200;
576 NUgrid *x_grid = create_center_grid (-3.0, 4.0, 7.5, Mx);
577 NUgrid *y_grid = create_center_grid (-1.0, 9.0, 3.5, My);
578 NUgrid *z_grid = create_center_grid (-1.8, 2.0, 2.8, Mz);
579 complex_double *data = malloc (sizeof(complex_double)*Mx*My*Mz);
580 for (int ix=0; ix<Mx; ix++)
581 for (int iy=0; iy<My; iy++)
582 for (int iz=0; iz<Mz; iz++)
583 data[(ix*My+iy)*Mz+iz] = -1.0+2.0*drand48() + 1.0if*(-1.0+2.0*drand48());
584
585 BCtype_z xBC, yBC, zBC;
586 xBC.lCode = PERIODIC; xBC.rCode = PERIODIC;
587 yBC.lCode = PERIODIC; yBC.rCode = PERIODIC;
588 zBC.lCode = PERIODIC; zBC.rCode = PERIODIC;
589
590 NUBspline_3d_z *spline = create_NUBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
591 complex_double val, grad[3], hess[9];
592 clock_t start, end, rstart, rend;
593 rstart = clock();
594 for (int i=0; i<10000000; i++) {
595 double x = x_grid->start+ 0.9999*drand48()*(x_grid->end - x_grid->start);
596 double y = y_grid->start+ 0.9999*drand48()*(y_grid->end - y_grid->start);
597 double z = z_grid->start+ 0.9999*drand48()*(z_grid->end - z_grid->start);
598 }
599 rend = clock();
600 start = clock();
601 for (int i=0; i<10000000; i++) {
602 double x = x_grid->start+ 0.9999*drand48()*(x_grid->end - x_grid->start);
603 double y = y_grid->start+ 0.9999*drand48()*(y_grid->end - y_grid->start);
604 double z = z_grid->start+ 0.9999*drand48()*(z_grid->end - z_grid->start);
605 eval_NUBspline_3d_z_vgh (spline, x, y, z, &val, grad, hess);
606 }
607 end = clock();
608 fprintf (stderr, "10,000,000 evalations in %f seconds.\n",
609 (double)(end-start-(rend-rstart))/(double)CLOCKS_PER_SEC);
610 }
611
612
613 void
TestNUB_2d_d()614 TestNUB_2d_d()
615 {
616 int Mx=30, My=35;
617 NUgrid *x_grid = create_center_grid (-3.0, 4.0, 7.5, Mx);
618 NUgrid *y_grid = create_center_grid (-1.0, 9.0, 3.5, My);
619 double data[Mx*My];
620 for (int ix=0; ix<Mx; ix++)
621 for (int iy=0; iy<My; iy++)
622 data[ix*My+iy] = -1.0+2.0*drand48();
623
624 BCtype_d xBC, yBC;
625 xBC.lCode = PERIODIC;
626 yBC.lCode = PERIODIC;
627 // xBC.lCode = FLAT; xBC.rCode = FLAT;
628 // yBC.lCode = FLAT; yBC.rCode = FLAT;
629
630
631
632 NUBspline_2d_d *spline = create_NUBspline_2d_d (x_grid, y_grid, xBC, yBC, data);
633
634 int xFine = 400;
635 int yFine = 400;
636 FILE *fout = fopen ("2d_d.dat", "w");
637 double xi = x_grid->start;
638 double xf = x_grid->end;// + x_grid->points[1] - x_grid->points[0];
639 double yi = y_grid->start;
640 double yf = y_grid->end;// + y_grid->points[1] - y_grid->points[0];
641 for (int ix=0; ix<xFine; ix++) {
642 double x = xi+ (double)ix/(double)(xFine)*(xf-xi);
643 for (int iy=0; iy<yFine; iy++) {
644 double y = yi + (double)iy/(double)(yFine)*(yf-yi);
645 double val;
646 eval_NUBspline_2d_d (spline, x, y, &val);
647 fprintf (fout, "%1.16e ", val);
648 }
649 fprintf (fout, "\n");
650 }
651 fclose (fout);
652 }
653
main()654 int main()
655 {
656 // TestCenterGrid();
657 // TestGeneralGrid();
658 // GridSpeedTest();
659 // TestNUBasis();
660 // TestNUBasis();
661 TestNUBspline_d();
662 // TestNUB_2d_s();
663 // TestNUB_2d_c();
664 // TestNUB_3d_c();
665 // SpeedNUB_3d_s();
666 // TestNUB_2d_d();
667 // TestNUB_3d_d();
668 // TestNUB_3d_z();
669 //SpeedNUB_3d_z();
670 // bool passed = TestNUB_1d_s();
671 }
672
673