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 "multi_nubspline.h"
9 #include "bspline.h"
10 #include "nubspline.h"
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <math.h>
14 #include <time.h>
15 #ifdef _OPENMP
16 #include <omp.h>
17 #endif
18
19 double drand48();
20
get_time()21 inline double get_time()
22 {
23 #ifdef _OPENMP
24 fprintf(stderr, "Using omp_get_wtime().\n");
25 return omp_get_wtime();
26 #else
27 return (double)clock() / (double)CLOCKS_PER_SEC;
28 #endif
29 }
30
31
diff(double a,double b,double tol)32 inline double diff (double a, double b, double tol)
33 {
34 if (fabs(a-b) > tol)
35 return 1;
36 else
37 return 0;
38 }
39
40
41 int
test_3d_double_all()42 test_3d_double_all()
43 {
44 int Nx=73; int Ny=91; int Nz = 29;
45 int num_splines = 128;
46
47 Ugrid x_grid, y_grid, z_grid;
48 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
49 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
50 z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
51
52 BCtype_d xBC, yBC, zBC;
53 xBC.lCode = xBC.rCode = PERIODIC;
54 yBC.lCode = yBC.rCode = PERIODIC;
55 zBC.lCode = zBC.rCode = PERIODIC;
56
57 // First, create splines the normal way
58 UBspline_3d_d* norm_splines[num_splines];
59 multi_UBspline_3d_d *multi_spline;
60
61 // First, create multispline
62 multi_spline = create_multi_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC,
63 num_splines);
64
65 double data[Nx*Ny*Nz];
66 // Now, create normal splines and set multispline data
67 for (int i=0; i<num_splines; i++) {
68 for (int j=0; j<Nx*Ny*Nz; j++)
69 data[j] = (drand48()-0.5);// + (drand48()-0.5)*1.0i;
70 norm_splines[i] = create_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
71 set_multi_UBspline_3d_d (multi_spline, i, data);
72 }
73
74 // Now, test random values
75 int num_vals = 100;
76 double multi_vals[num_splines], norm_vals[num_splines];
77 double multi_grads[3*num_splines], norm_grads[3*num_splines];
78 double multi_lapl[num_splines], norm_lapl[num_splines];
79 double multi_hess[9*num_splines], norm_hess[9*num_splines];
80 for (int i=0; i<num_vals; i++) {
81 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
82 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
83 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
84
85 //////////////////////
86 // Check V routine //
87 //////////////////////
88 eval_multi_UBspline_3d_d (multi_spline, x, y, z, multi_vals);
89 for (int j=0; j<num_splines; j++)
90 eval_UBspline_3d_d (norm_splines[j], x, y, z, &(norm_vals[j]));
91 for (int j=0; j<num_splines; j++) {
92 // Check value
93 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
94 return -1;
95 }
96
97
98
99 ///////////////////////
100 // Check VG routine //
101 ///////////////////////
102 eval_multi_UBspline_3d_d_vg (multi_spline, x, y, z,
103 multi_vals, multi_grads);
104 for (int j=0; j<num_splines; j++)
105 eval_UBspline_3d_d_vg (norm_splines[j], x, y, z, &(norm_vals[j]),
106 &(norm_grads[3*j]));
107 for (int j=0; j<num_splines; j++) {
108 // Check value
109 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
110 return -1;
111
112 // Check gradients
113 for (int n=0; n<3; n++)
114 if (diff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-12))
115 return -2;
116 }
117
118
119 ///////////////////////
120 // Check VGL routine //
121 ///////////////////////
122 eval_multi_UBspline_3d_d_vgl (multi_spline, x, y, z,
123 multi_vals, multi_grads, multi_lapl);
124 for (int j=0; j<num_splines; j++)
125 eval_UBspline_3d_d_vgl (norm_splines[j], x, y, z, &(norm_vals[j]),
126 &(norm_grads[3*j]), &(norm_lapl[j]));
127 for (int j=0; j<num_splines; j++) {
128 // Check value
129 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
130 return -3;
131
132 // Check gradients
133 for (int n=0; n<3; n++)
134 if (diff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-10))
135 return -4;
136
137 // Check laplacian
138 if (diff (norm_lapl[j], multi_lapl[j], 1.0e-10))
139 return -5;
140 }
141
142
143 ///////////////////////
144 // Check VGH routine //
145 ///////////////////////
146 eval_multi_UBspline_3d_d_vgh (multi_spline, x, y, z,
147 multi_vals, multi_grads, multi_hess);
148 for (int j=0; j<num_splines; j++)
149 eval_UBspline_3d_d_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
150 &(norm_grads[3*j]), &(norm_hess[9*j]));
151 for (int j=0; j<num_splines; j++) {
152 // Check value
153 if (diff(norm_vals[j], multi_vals[j], 1.0e-12))
154 return -6;
155
156 // Check gradients
157 for (int n=0; n<3; n++)
158 if (diff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-12))
159 return -7;
160
161 // Check hessian
162 for (int n=0; n<9; n++)
163 if (diff (norm_hess[9*j+n], multi_hess[9*j+n], 1.0e-10))
164 return -8;
165 }
166 }
167 return 0;
168 }
169
170
171
172
173 /////////////////////////////////////////////
174 // Single-precision complex test functions //
175 /////////////////////////////////////////////
176 inline int
cdiff(complex_float a,complex_float b,double tol)177 cdiff (complex_float a, complex_float b, double tol)
178 {
179 double rdiff = fabs(creal(a) - creal(b));
180 double idiff = fabs(cimag(a) - cimag(b));
181 if (rdiff > tol || idiff > tol)
182 return 1;
183 else
184 return 0;
185 }
186
187
188
189 /////////////////////////////////////////////
190 // Double-precision complex test functions //
191 /////////////////////////////////////////////
test_complex_double()192 void test_complex_double()
193 {
194 int Nx=73; int Ny=91; int Nz = 29;
195 int num_splines = 128;
196
197 Ugrid x_grid, y_grid, z_grid;
198 x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
199 y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
200 z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
201
202 BCtype_z xBC, yBC, zBC;
203 xBC.lCode = xBC.rCode = PERIODIC;
204 yBC.lCode = yBC.rCode = PERIODIC;
205 zBC.lCode = zBC.rCode = PERIODIC;
206
207 // First, create splines the normal way
208 UBspline_3d_z* norm_splines[num_splines];
209 multi_UBspline_3d_z *multi_spline;
210
211 // First, create multispline
212 multi_spline = create_multi_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC,
213 num_splines);
214
215 complex_double data[Nx*Ny*Nz];
216 // Now, create normal splines and set multispline data
217 for (int i=0; i<num_splines; i++) {
218 for (int j=0; j<Nx*Ny*Nz; j++)
219 data[j] = (drand48()-0.5);// + (drand48()-0.5)*1.0i;
220 norm_splines[i] = create_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
221 set_multi_UBspline_3d_z (multi_spline, i, data);
222 }
223
224 fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
225 creal(norm_splines[19]->coefs[227]),
226 cimag(norm_splines[19]->coefs[227]));
227 fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
228 creal(multi_spline->coefs[19+227*multi_spline->z_stride]),
229 cimag(multi_spline->coefs[19+227*multi_spline->z_stride]));
230 //return;
231
232 // Now, test random values
233 int num_vals = 100;
234 complex_double multi_vals[num_splines], norm_vals[num_splines];
235 for (int i=0; i<num_vals; i++) {
236 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
237 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
238 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
239
240 eval_multi_UBspline_3d_z (multi_spline, x, y, z, multi_vals);
241 for (int j=0; j<num_splines; j++)
242 eval_UBspline_3d_z (norm_splines[j], x, y, z, &(norm_vals[j]));
243 for (int j=0; j<num_splines; j++) {
244 double rdiff = creal(norm_vals[j]) - creal(multi_vals[j]);
245 double idiff = cimag(norm_vals[j]) - cimag(multi_vals[j]);
246 if (fabs(rdiff) > 1.0e-12 || fabs(idiff) > 1.0e-12) {
247 fprintf (stderr, "Error! norm_vals[j] = %1.14e + %1.14ei\n",
248 creal(norm_vals[j]), cimag(norm_vals[j]));
249 fprintf (stderr, " multi_vals[j] = %1.14e + %1.14ei\n",
250 creal(multi_vals[j]), cimag(multi_vals[j]));
251 }
252 }
253 }
254
255 num_vals = 100000;
256
257 // Now do timing
258 double norm_start, norm_end, multi_start, multi_end, rand_start, rand_end;
259 rand_start = get_time();
260 for (int i=0; i<num_vals; i++) {
261 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
262 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
263 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
264 }
265 rand_end = get_time();
266
267 norm_start = get_time();
268 for (int i=0; i<num_vals; i++) {
269 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
270 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
271 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
272
273 for (int j=0; j<num_splines; j++)
274 eval_UBspline_3d_z (norm_splines[j], x, y, z, &(norm_vals[j]));
275 }
276 norm_end = get_time();
277
278 multi_start = get_time();
279 for (int i=0; i<num_vals; i++) {
280 double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
281 double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
282 double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
283 eval_multi_UBspline_3d_z (multi_spline, x, y, z, multi_vals);
284 }
285 multi_end = get_time();
286
287 fprintf (stderr, "Normal spline time = %1.5f\n",
288 (double)(norm_end-norm_start+rand_start-rand_end)/CLOCKS_PER_SEC);
289 fprintf (stderr, "Multi spline time = %1.5f\n",
290 (double)(multi_end-multi_start+rand_start-rand_end)/CLOCKS_PER_SEC);
291
292 }
293
294
295 inline int
zdiff(complex_double a,complex_double b,double tol)296 zdiff (complex_double a, complex_double b, double tol)
297 {
298 double rdiff = fabs(creal(a) - creal(b));
299 double idiff = fabs(cimag(a) - cimag(b));
300 if (rdiff > tol || idiff > tol)
301 return 1;
302 else
303 return 0;
304 }
305
306
307 // int
308 // test_3d_complex_double_all()
309 // {
310 // int Nx=73; int Ny=91; int Nz = 29;
311 // int num_splines = 21;
312
313 // Ugrid x_grid, y_grid, z_grid;
314 // x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
315 // y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
316 // z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
317
318 // BCtype_z xBC, yBC, zBC;
319 // xBC.lCode = xBC.rCode = PERIODIC;
320 // yBC.lCode = yBC.rCode = PERIODIC;
321 // zBC.lCode = zBC.rCode = PERIODIC;
322
323 // // First, create splines the normal way
324 // UBspline_3d_z* norm_splines[num_splines];
325 // multi_UBspline_3d_z *multi_spline;
326
327 // // First, create multispline
328 // multi_spline = create_multi_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC,
329 // num_splines);
330
331 // complex_double data[Nx*Ny*Nz];
332 // // Now, create normal splines and set multispline data
333 // for (int i=0; i<num_splines; i++) {
334 // for (int j=0; j<Nx*Ny*Nz; j++)
335 // data[j] = (drand48()-0.5) + (drand48()-0.5)*1.0i;
336 // norm_splines[i] = create_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
337 // set_multi_UBspline_3d_z (multi_spline, i, data);
338 // }
339
340 // // fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
341 // // creal(norm_splines[19]->coefs[227]),
342 // // cimag(norm_splines[19]->coefs[227]));
343 // // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
344 // // creal(multi_spline->coefs[19+227*multi_spline->z_stride]),
345 // // cimag(multi_spline->coefs[19+227*multi_spline->z_stride]));
346
347 // // Now, test random values
348 // int num_vals = 100;
349 // complex_double multi_vals[num_splines], norm_vals[num_splines];
350 // complex_double multi_grads[3*num_splines], norm_grads[3*num_splines];
351 // complex_double multi_lapl[num_splines], norm_lapl[num_splines];
352 // complex_double multi_hess[9*num_splines], norm_hess[9*num_splines];
353 // for (int i=0; i<num_vals; i++) {
354 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
355 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
356 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
357
358 // ///////////////////////
359 // // Check value only //
360 // ///////////////////////
361 // eval_multi_UBspline_3d_z (multi_spline, x, y, z, multi_vals);
362 // for (int j=0; j<num_splines; j++)
363 // eval_UBspline_3d_z (norm_splines[j], x, y, z, &(norm_vals[j]));
364 // for (int j=0; j<num_splines; j++)
365 // // Check value
366 // if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
367 // return -1;
368
369 // ///////////////////////
370 // // Check VG routine //
371 // ///////////////////////
372 // eval_multi_UBspline_3d_z_vg (multi_spline, x, y, z,
373 // multi_vals, multi_grads);
374 // for (int j=0; j<num_splines; j++)
375 // eval_UBspline_3d_z_vg (norm_splines[j], x, y, z, &(norm_vals[j]),
376 // &(norm_grads[3*j]));
377 // for (int j=0; j<num_splines; j++) {
378 // // Check value
379 // if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
380 // return -2;
381
382 // // Check gradients
383 // for (int n=0; n<3; n++)
384 // if (zdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-12))
385 // return -3;
386 // }
387
388
389 // ///////////////////////
390 // // Check VGL routine //
391 // ///////////////////////
392 // eval_multi_UBspline_3d_z_vgl (multi_spline, x, y, z,
393 // multi_vals, multi_grads, multi_lapl);
394 // for (int j=0; j<num_splines; j++)
395 // eval_UBspline_3d_z_vgl (norm_splines[j], x, y, z, &(norm_vals[j]),
396 // &(norm_grads[3*j]), &(norm_lapl[j]));
397 // for (int j=0; j<num_splines; j++) {
398 // // Check value
399 // if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
400 // return -4;
401
402 // // Check gradients
403 // for (int n=0; n<3; n++)
404 // if (zdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-10))
405 // return -5;
406
407 // // Check laplacian
408 // if (zdiff (norm_lapl[j], multi_lapl[j], 1.0e-10))
409 // return -6;
410 // }
411
412
413 // ///////////////////////
414 // // Check VGH routine //
415 // ///////////////////////
416 // eval_multi_UBspline_3d_z_vgh (multi_spline, x, y, z,
417 // multi_vals, multi_grads, multi_hess);
418 // for (int j=0; j<num_splines; j++)
419 // eval_UBspline_3d_z_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
420 // &(norm_grads[3*j]), &(norm_hess[9*j]));
421 // for (int j=0; j<num_splines; j++) {
422 // // Check value
423 // if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12))
424 // return -7;
425
426 // // Check gradients
427 // for (int n=0; n<3; n++)
428 // if (zdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-12))
429 // return -8;
430
431 // // Check hessian
432 // for (int n=0; n<9; n++)
433 // if (zdiff (norm_hess[9*j+n], multi_hess[9*j+n], 1.0e-10)) {
434 // fprintf (stderr, "\nj = %d n = %d \n", j, n);
435 // fprintf (stderr, "norm_hess[j] = %1.14e + %1.14ei\n",
436 // creal(norm_hess[9*j+n]), cimag(norm_hess[9*j+n]));
437 // fprintf (stderr, "multi_hess[j] = %1.14e + %1.15ei\n",
438 // creal(multi_hess[9*j+n]), cimag(multi_hess[9*j+n]));
439 // return -9;
440 // }
441 // }
442 // }
443 // return 0;
444 // }
445
446
447 // void test_complex_double_vgh()
448 // {
449 // int Nx=73; int Ny=91; int Nz = 29;
450 // int num_splines = 128;
451
452 // Ugrid x_grid, y_grid, z_grid;
453 // x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
454 // y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
455 // z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
456
457 // BCtype_z xBC, yBC, zBC;
458 // xBC.lCode = xBC.rCode = PERIODIC;
459 // yBC.lCode = yBC.rCode = PERIODIC;
460 // zBC.lCode = zBC.rCode = PERIODIC;
461
462 // // First, create splines the normal way
463 // UBspline_3d_z* norm_splines[num_splines];
464 // multi_UBspline_3d_z *multi_spline;
465
466 // // First, create multispline
467 // multi_spline = create_multi_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC,
468 // num_splines);
469
470 // complex_double data[Nx*Ny*Nz];
471 // // Now, create normal splines and set multispline data
472 // for (int i=0; i<num_splines; i++) {
473 // for (int j=0; j<Nx*Ny*Nz; j++)
474 // data[j] = (drand48()-0.5) + (drand48()-0.5)*1.0i;
475 // norm_splines[i] = create_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
476 // set_multi_UBspline_3d_z (multi_spline, i, data);
477 // }
478
479 // fprintf (stderr, "norm coef = %1.14e + %1.14ei\n",
480 // creal(norm_splines[19]->coefs[227]),
481 // cimag(norm_splines[19]->coefs[227]));
482 // fprintf (stderr, "multi coef = %1.14e + %1.14ei\n",
483 // creal(multi_spline->coefs[19+227*multi_spline->z_stride]),
484 // cimag(multi_spline->coefs[19+227*multi_spline->z_stride]));
485
486 // // Now, test random values
487 // int num_vals = 100;
488 // complex_double multi_vals[num_splines], norm_vals[num_splines];
489 // complex_double multi_grads[3*num_splines], norm_grads[3*num_splines];
490 // complex_double multi_lapl[num_splines], norm_lapl[num_splines];
491 // complex_double multi_hess[9*num_splines], norm_hess[9*num_splines];
492 // for (int i=0; i<num_vals; i++) {
493 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
494 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
495 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
496
497 // ///////////////////////
498 // // Check VGH routine //
499 // ///////////////////////
500 // eval_multi_UBspline_3d_z_vgh (multi_spline, x, y, z,
501 // multi_vals, multi_grads, multi_hess);
502 // for (int j=0; j<num_splines; j++)
503 // eval_UBspline_3d_z_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
504 // &(norm_grads[3*j]), &(norm_hess[9*j]));
505 // for (int j=0; j<num_splines; j++) {
506 // // Check value
507 // if (zdiff(norm_vals[j], multi_vals[j], 1.0e-12)) {
508 // fprintf (stderr, "Error! norm_vals[j] = %1.14e + %1.14ei\n",
509 // creal(norm_vals[j]), cimag(norm_vals[j]));
510 // fprintf (stderr, " multi_vals[j] = %1.14e + %1.14ei\n",
511 // creal(multi_vals[j]), cimag(multi_vals[j]));
512 // }
513 // // Check gradients
514 // for (int n=0; n<3; n++) {
515 // if (zdiff (norm_grads[3*j+n], multi_grads[3*j+n], 1.0e-12)) {
516 // fprintf (stderr, "n=%d\n", n);
517 // fprintf (stderr, "Error! norm_grads[j] = %1.14e + %1.14ei\n",
518 // creal(norm_grads[3*j+n]), cimag(norm_grads[3*j+n]));
519 // fprintf (stderr, " multi_grads[j] = %1.14e + %1.14ei\n",
520 // creal(multi_grads[3*j+n]), cimag(multi_grads[3*j+n]));
521 // }
522 // }
523 // // Check hessian
524 // for (int n=0; n<9; n++) {
525 // if (zdiff (norm_hess[9*j+n], multi_hess[9*j+n], 1.0e-10)) {
526 // fprintf (stderr, "Error! norm_hess[j] = %1.14e + %1.14ei\n",
527 // creal(norm_hess[9*j+n]), cimag(norm_hess[9*j+n]));
528 // fprintf (stderr, " multi_hess[j] = %1.14e + %1.14ei\n",
529 // creal(multi_hess[9*j+n]), cimag(multi_hess[9*j+n]));
530 // }
531 // }
532 // }
533 // }
534
535 // num_vals = 100000;
536
537 // // Now do timing
538 // clock_t norm_start, norm_end, multi_start, multi_end, rand_start, rand_end;
539 // rand_start = get_time();
540 // for (int i=0; i<num_vals; i++) {
541 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
542 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
543 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
544 // }
545 // rand_end = get_time();
546
547 // norm_start = get_time();
548 // for (int i=0; i<num_vals; i++) {
549 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
550 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
551 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
552
553 // for (int j=0; j<num_splines; j++)
554 // eval_UBspline_3d_z_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
555 // &(norm_grads[3*j]), &(norm_hess[9*j]));
556 // }
557 // norm_end = get_time();
558
559 // multi_start = get_time();
560 // for (int i=0; i<num_vals; i++) {
561 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
562 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
563 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
564 // eval_multi_UBspline_3d_z_vgh (multi_spline, x, y, z, multi_vals, multi_grads, multi_hess);
565 // }
566 // multi_end = get_time();
567
568 // fprintf (stderr, "Normal spline time = %1.5f\n",
569 // (double)(norm_end-norm_start+rand_start-rand_end)/CLOCKS_PER_SEC);
570 // fprintf (stderr, "Multi spline time = %1.5f\n",
571 // (double)(multi_end-multi_start+rand_start-rand_end)/CLOCKS_PER_SEC);
572
573 // }
574
575
576 // void test_double()
577 // {
578 // int Nx=73; int Ny=91; int Nz = 29;
579 // int num_splines = 201;
580
581 // Ugrid x_grid, y_grid, z_grid;
582 // x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
583 // y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
584 // z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
585
586 // BCtype_d xBC, yBC, zBC;
587 // xBC.lCode = xBC.rCode = PERIODIC;
588 // yBC.lCode = yBC.rCode = PERIODIC;
589 // zBC.lCode = zBC.rCode = PERIODIC;
590
591 // // First, create splines the normal way
592 // UBspline_3d_d* norm_splines[num_splines];
593 // multi_UBspline_3d_d *multi_spline;
594
595 // // First, create multispline
596 // multi_spline = create_multi_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC,
597 // num_splines);
598
599 // double data[Nx*Ny*Nz];
600 // // Now, create normal splines and set multispline data
601 // for (int i=0; i<num_splines; i++) {
602 // for (int j=0; j<Nx*Ny*Nz; j++)
603 // data[j] = (drand48()-0.5) + (drand48()-0.5)*1.0i;
604 // norm_splines[i] = create_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
605 // set_multi_UBspline_3d_d (multi_spline, i, data);
606 // }
607
608 // fprintf (stderr, "norm coef = %1.14e\n",
609 // norm_splines[19]->coefs[227]);
610 // fprintf (stderr, "multi coef = %1.14e\n",
611 // multi_spline->coefs[19+227*multi_spline->z_stride]);
612
613 // // Now, test random values
614 // int num_vals = 100;
615 // double multi_vals[num_splines], norm_vals[num_splines];
616 // for (int i=0; i<num_vals; i++) {
617 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
618 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
619 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
620
621 // eval_multi_UBspline_3d_d (multi_spline, x, y, z,
622 // multi_vals);
623 // for (int j=0; j<num_splines; j++)
624 // eval_UBspline_3d_d (norm_splines[j], x, y, z, &(norm_vals[j]));
625 // for (int j=0; j<num_splines; j++) {
626 // // Check value
627 // double diff = norm_vals[j] - multi_vals[j];
628 // if (fabs(diff) > 1.0e-12) {
629 // fprintf (stderr, "Error! norm_vals[j] = %1.14e\n",
630 // norm_vals[j]);
631 // fprintf (stderr, " multi_vals[j] = %1.14e\n",
632 // multi_vals[j]);
633 // }
634 // }
635 // }
636
637 // num_vals = 100000;
638
639 // // Now do timing
640 // clock_t norm_start, norm_end, multi_start, multi_end, rand_start, rand_end;
641 // rand_start = get_time();
642 // for (int i=0; i<num_vals; i++) {
643 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
644 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
645 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
646 // }
647 // rand_end = get_time();
648
649 // norm_start = get_time();
650 // for (int i=0; i<num_vals; i++) {
651 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
652 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
653 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
654
655 // for (int j=0; j<num_splines; j++)
656 // eval_UBspline_3d_d (norm_splines[j], x, y, z, &(norm_vals[j]));
657 // }
658 // norm_end = get_time();
659
660 // multi_start = get_time();
661 // for (int i=0; i<num_vals; i++) {
662 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
663 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
664 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
665 // eval_multi_UBspline_3d_d (multi_spline, x, y, z, multi_vals);
666 // }
667 // multi_end = get_time();
668
669 // fprintf (stderr, "Normal spline time = %1.5f\n",
670 // (double)(norm_end-norm_start+rand_start-rand_end)/CLOCKS_PER_SEC);
671 // fprintf (stderr, "Multi spline time = %1.5f\n",
672 // (double)(multi_end-multi_start+rand_start-rand_end)/CLOCKS_PER_SEC);
673
674 // }
675
676
677
678 // void test_double_vgh()
679 // {
680 // int Nx=73; int Ny=91; int Nz = 29;
681 // int num_splines = 128;
682
683 // Ugrid x_grid, y_grid, z_grid;
684 // x_grid.start = 3.1; x_grid.end = 9.1; x_grid.num = Nx;
685 // y_grid.start = 8.7; y_grid.end = 12.7; y_grid.num = Ny;
686 // z_grid.start = 4.5; z_grid.end = 9.3; z_grid.num = Nz;
687
688 // BCtype_d xBC, yBC, zBC;
689 // xBC.lCode = xBC.rCode = PERIODIC;
690 // yBC.lCode = yBC.rCode = PERIODIC;
691 // zBC.lCode = zBC.rCode = PERIODIC;
692
693 // // First, create splines the normal way
694 // UBspline_3d_d* norm_splines[num_splines];
695 // multi_UBspline_3d_d *multi_spline;
696
697 // // First, create multispline
698 // multi_spline = create_multi_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC,
699 // num_splines);
700
701 // double data[Nx*Ny*Nz];
702 // // Now, create normal splines and set multispline data
703 // for (int i=0; i<num_splines; i++) {
704 // for (int j=0; j<Nx*Ny*Nz; j++)
705 // data[j] = (drand48()-0.5) + (drand48()-0.5)*1.0i;
706 // norm_splines[i] = create_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC, data);
707 // set_multi_UBspline_3d_d (multi_spline, i, data);
708 // }
709
710 // fprintf (stderr, "norm coef = %1.14e\n",
711 // norm_splines[19]->coefs[227]);
712 // fprintf (stderr, "multi coef = %1.14e\n",
713 // multi_spline->coefs[19+227*multi_spline->z_stride]);
714
715 // // Now, test random values
716 // int num_vals = 100;
717 // double multi_vals[num_splines], norm_vals[num_splines];
718 // double multi_grads[3*num_splines], norm_grads[3*num_splines];
719 // double multi_hess[9*num_splines], norm_hess[9*num_splines];
720 // for (int i=0; i<num_vals; i++) {
721 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
722 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
723 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
724
725 // eval_multi_UBspline_3d_d_vgh (multi_spline, x, y, z,
726 // multi_vals, multi_grads, multi_hess);
727 // for (int j=0; j<num_splines; j++)
728 // eval_UBspline_3d_d_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
729 // &(norm_grads[3*j]), &(norm_hess[9*j]));
730 // for (int j=0; j<num_splines; j++) {
731 // // Check value
732 // double diff = norm_vals[j] - multi_vals[j];
733 // if (fabs(diff) > 1.0e-12) {
734 // fprintf (stderr, "j = %d\n", j);
735 // fprintf (stderr, "Error! norm_vals[j] = %1.14e\n",
736 // norm_vals[j]);
737 // fprintf (stderr, " multi_vals[j] = %1.14e\n",
738 // multi_vals[j]);
739 // }
740 // // Check gradients
741 // for (int n=0; n<3; n++) {
742 // diff = norm_grads[3*j+n] - multi_grads[3*j+n];
743 // if (fabs(diff) > 1.0e-12) {
744 // fprintf (stderr, "n=%d\n", n);
745 // fprintf (stderr, "Error! norm_grads[j] = %1.14e\n",
746 // norm_grads[3*j+n]);
747 // fprintf (stderr, " multi_grads[j] = %1.14e\n",
748 // multi_grads[3*j+n]);
749 // }
750 // }
751 // // Check hessian
752 // for (int n=0; n<9; n++) {
753 // diff = norm_hess[9*j+n] - multi_hess[9*j+n];
754 // if (fabs(diff) > 1.0e-10) {
755 // fprintf (stderr, "Error! norm_hess[j] = %1.14e\n",
756 // norm_hess[9*j+n]);
757 // fprintf (stderr, " multi_hess[j] = %1.14e\n",
758 // multi_hess[9*j+n]);
759 // }
760 // }
761 // }
762 // }
763
764 // num_vals = 100000;
765
766 // // Now do timing
767 // clock_t norm_start, norm_end, multi_start, multi_end, rand_start, rand_end;
768 // rand_start = get_time();
769 // for (int i=0; i<num_vals; i++) {
770 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
771 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
772 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
773 // }
774 // rand_end = get_time();
775
776 // norm_start = get_time();
777 // for (int i=0; i<num_vals; i++) {
778 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
779 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
780 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
781
782 // for (int j=0; j<num_splines; j++)
783 // eval_UBspline_3d_d_vgh (norm_splines[j], x, y, z, &(norm_vals[j]),
784 // &(norm_grads[3*j]), &(norm_hess[9*j]));
785 // }
786 // norm_end = get_time();
787
788 // multi_start = get_time();
789 // for (int i=0; i<num_vals; i++) {
790 // double rx = drand48(); double x = rx*x_grid.start + (1.0-rx)*x_grid.end;
791 // double ry = drand48(); double y = ry*y_grid.start + (1.0-ry)*y_grid.end;
792 // double rz = drand48(); double z = rz*z_grid.start + (1.0-rz)*z_grid.end;
793 // eval_multi_UBspline_3d_d_vgh (multi_spline, x, y, z, multi_vals, multi_grads, multi_hess);
794 // }
795 // multi_end = get_time();
796
797 // fprintf (stderr, "Normal spline time = %1.5f\n",
798 // (double)(norm_end-norm_start+rand_start-rand_end)/CLOCKS_PER_SEC);
799 // fprintf (stderr, "Multi spline time = %1.5f\n",
800 // (double)(multi_end-multi_start+rand_start-rand_end)/CLOCKS_PER_SEC);
801
802 // }
803
PrintPassFail(int code)804 void PrintPassFail (int code)
805 {
806 char green[100], normal[100], red[100];
807 snprintf (green, 100, "%c[0;32;47m", 0x1b);
808 snprintf (normal, 100, "%c[0;30;47m", 0x1b);
809 snprintf (red, 100, "%c[0;31;47m", 0x1b);
810
811 if (code == 0)
812 fprintf (stderr, "PASSED\n");
813 else
814 fprintf (stderr, "FAILED: code = %d\n", code);
815 }
816
817
main()818 main()
819 {
820 int code;
821 //test_complex_double();
822 //test_complex_double_vgh();
823
824 // fprintf (stderr, "Testing 1D complex double-precision multiple nonuniform cubic B-spline routines: ");
825 // code = test_1d_NUB_complex_double_all(); PrintPassFail (code);
826
827 // fprintf (stderr, "Testing 1D real single-precision multiple cubic B-spline routines: ");
828 // code = test_1d_float_all(); PrintPassFail (code);
829 // fprintf (stderr, "Testing 2D real single-precision multiple cubic B-spline routines: ");
830 // code = test_2d_float_all(); PrintPassFail (code);
831 // fprintf (stderr, "Testing 3D real single-precision multiple cubic B-spline routines: ");
832 // code = test_3d_float_all(); PrintPassFail (code);
833
834 // fprintf (stderr, "Testing 1D real double-precision multiple cubic B-spline routines: ");
835 // code = test_1d_double_all(); PrintPassFail (code);
836 // fprintf (stderr, "Testing 2D real double-precision multiple cubic B-spline routines: ");
837 // code = test_2d_double_all(); PrintPassFail (code);
838 fprintf (stderr, "Testing 3D real double-precision multiple cubic B-spline routines: ");
839 code = test_3d_double_all(); PrintPassFail (code);
840
841 // fprintf (stderr, "Testing 1D complex single-precision multiple cubic B-spline routines: ");
842 // code = test_1d_complex_float_all(); PrintPassFail (code);
843 // fprintf (stderr, "Testing 2D complex single-precision multiple cubic B-spline routines: ");
844 // code = test_2d_complex_float_all(); PrintPassFail (code);
845 // fprintf (stderr, "Testing 3D complex single-precision multiple cubic B-spline routines: ");
846 // code = test_3d_complex_float_all(); PrintPassFail (code);
847
848 // fprintf (stderr, "Testing 1D complex double-precision multiple cubic B-spline routines: ");
849 // code = test_1d_complex_double_all(); PrintPassFail (code);
850 // fprintf (stderr, "Testing 2D complex double-precision multiple cubic B-spline routines: ");
851 // code = test_2d_complex_double_all(); PrintPassFail (code);
852 // fprintf (stderr, "Testing 3D complex double-precision multiple cubic B-spline routines: ");
853 // code = test_3d_complex_double_all(); PrintPassFail (code);
854
855
856 //test_double();
857 //test_double_vgh();
858 }
859