Krita Source Code Documentation
Loading...
Searching...
No Matches
multi_bspline_create.cpp
Go to the documentation of this file.
1
2// einspline: a library for creating and evaluating B-splines //
3// Copyright (C) 2007 Kenneth P. Esler, Jr. //
4// //
5// This program is free software; you can redistribute it and/or modify //
6// it under the terms of the GNU General Public License as published by //
7// the Free Software Foundation; either version 2 of the License, or //
8// (at your option) any later version. //
9// //
10// This program is distributed in the hope that it will be useful, //
11// but WITHOUT ANY WARRANTY; without even the implied warranty of //
12// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
13// GNU General Public License for more details. //
14// //
15// You should have received a copy of the GNU General Public License //
16// along with this program; if not, write to the Free Software //
17// Foundation, Inc., 51 Franklin Street, Fifth Floor, //
18// Boston, MA 02110-1301 USA //
20
22#ifndef _XOPEN_SOURCE
23#define _XOPEN_SOURCE 600
24#endif
25#ifndef __USE_XOPEN2K
26 #define __USE_XOPEN2K
27#endif
28#include <stdlib.h>
29#include <stdio.h>
30
31int posix_memalign(void **memptr, size_t alignment, size_t size);
32
38void init_sse_data();
39
40void
42 double *data, intptr_t dstride,
43 double *coefs, intptr_t cstride);
44
45void
46solve_deriv_interp_1d_s (float bands[], float coefs[],
47 int M, int cstride);
48
49// On input, bands should be filled with:
50// row 0 : abcdInitial from boundary conditions
51// rows 1:M: basis functions in first 3 cols, data in last
52// row M+1 : abcdFinal from boundary conditions
53// cstride gives the stride between values in coefs.
54// On exit, coefs with contain interpolating B-spline coefs
55void
56solve_periodic_interp_1d_s (float bands[], float coefs[],
57 int M, int cstride);
58
59// On input, bands should be filled with:
60// row 0 : abcdInitial from boundary conditions
61// rows 1:M: basis functions in first 3 cols, data in last
62// row M+1 : abcdFinal from boundary conditions
63// cstride gives the stride between values in coefs.
64// On exit, coefs with contain interpolating B-spline coefs
65void
66solve_antiperiodic_interp_1d_s (float bands[], float coefs[],
67 int M, int cstride);
68
69void
71 float *data, intptr_t dstride,
72 float *coefs, intptr_t cstride);
73
79
80// On input, bands should be filled with:
81// row 0 : abcdInitial from boundary conditions
82// rows 1:M: basis functions in first 3 cols, data in last
83// row M+1 : abcdFinal from boundary conditions
84// cstride gives the stride between values in coefs.
85// On exit, coefs with contain interpolating B-spline coefs
87create_multi_UBspline_1d_s (Ugrid x_grid, BCtype_s xBC, int num_splines)
88{
89 // Create new spline
90 multi_UBspline_1d_s* restrict spline = static_cast<multi_UBspline_1d_s*>(malloc(sizeof(multi_UBspline_1d_s)));
91 if (!spline) {
92 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_1d_s.\n");
93 abort();
94 }
95 spline->spcode = MULTI_U1D;
96 spline->tcode = SINGLE_REAL;
97 spline->xBC = xBC; spline->x_grid = x_grid;
98 spline->num_splines = num_splines;
99
100 // Setup internal variables
101 int Mx = x_grid.num;
102 int Nx;
103
104 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC) {
105 x_grid.delta = (x_grid.end-x_grid.start)/(double)(x_grid.num);
106 Nx = Mx+3;
107 }
108 else {
109 x_grid.delta = (x_grid.end-x_grid.start)/(double)(x_grid.num-1);
110 Nx = Mx+2;
111 }
112
113 int N = num_splines;
114#ifdef HAVE_SSE
115 if (N % 4)
116 N += 4 - (N % 4);
117#endif
118
119 spline->x_stride = N;
120 x_grid.delta_inv = 1.0/x_grid.delta;
121 spline->x_grid = x_grid;
122#ifndef HAVE_POSIX_MEMALIGN
123 spline->coefs = (float*)malloc (sizeof(float)*Nx*N);
124#else
125 posix_memalign ((void**)&spline->coefs, 64, (sizeof(float)*Nx*N));
126#endif
127#ifdef HAVE_SSE
128 init_sse_data();
129#endif
130 if (!spline->coefs) {
131 fprintf (stderr, "Out of memory allocating spline coefficient in create_multi_UBspline_1d_s.\n");
132 abort();
133 }
134
135
136 return spline;
137}
138
139void
141 float *data)
142{
143 float *coefs = spline->coefs + num;
144 int xs = spline->x_stride;
145 find_coefs_1d_s (spline->x_grid, spline->xBC, data, 1,
146 coefs, xs);
147}
148
149
152 BCtype_s xBC, BCtype_s yBC, int num_splines)
153{
154 // Create new spline
155 multi_UBspline_2d_s* restrict spline = static_cast<multi_UBspline_2d_s*>(malloc(sizeof(multi_UBspline_2d_s)));
156 if (!spline) {
157 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_2d_s.\n");
158 abort();
159 }
160 spline->spcode = MULTI_U2D;
161 spline->tcode = SINGLE_REAL;
162 spline->xBC = xBC;
163 spline->yBC = yBC;
164 spline->num_splines = num_splines;
165 // Setup internal variables
166 int Mx = x_grid.num;
167 int My = y_grid.num;
168 int Nx, Ny;
169
170 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)
171 Nx = Mx+3;
172 else
173 Nx = Mx+2;
174 x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
175
176 x_grid.delta_inv = 1.0/x_grid.delta;
177 spline->x_grid = x_grid;
178
179 if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)
180 Ny = My+3;
181 else
182 Ny = My+2;
183 y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
184 y_grid.delta_inv = 1.0/y_grid.delta;
185 spline->y_grid = y_grid;
186
187 int N = num_splines;
188#ifdef HAVE_SSE
189 if (N % 4)
190 N += 4 - (N % 4);
191#endif
192
193
194 spline->x_stride = Ny*N;
195 spline->y_stride = N;
196
197#ifndef HAVE_POSIX_MEMALIGN
198 spline->coefs = (float*)malloc (sizeof(float)*Nx*Ny*N);
199#else
200 posix_memalign ((void**)&spline->coefs, 64,
201 sizeof(float)*Nx*Ny*N);
202#endif
203#ifdef HAVE_SSE
205#endif
206 if (!spline->coefs) {
207 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_2d_s.\n");
208 abort();
209 }
210
211 return spline;
212}
213
214void
215set_multi_UBspline_2d_s (multi_UBspline_2d_s* spline, int num, float *data)
216{
217 int Mx = spline->x_grid.num;
218 int My = spline->y_grid.num;
219 int Nx, Ny;
220
221 if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)
222 Nx = Mx+3;
223 else
224 Nx = Mx+2;
225 if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)
226 Ny = My+3;
227 else
228 Ny = My+2;
229
230 float *coefs = spline->coefs + num;
231 int ys = spline->y_stride;
232 // First, solve in the X-direction
233 for (int iy=0; iy<My; iy++) {
234 intptr_t doffset = iy;
235 intptr_t coffset = iy*ys;
236 find_coefs_1d_s (spline->x_grid, spline->xBC, data+doffset, (intptr_t)My,
237 coefs+coffset, (intptr_t)Ny*ys);
238 }
239
240 // Now, solve in the Y-direction
241 for (int ix=0; ix<Nx; ix++) {
242 intptr_t doffset = ix*Ny*ys;
243 intptr_t coffset = ix*Ny*ys;
244 find_coefs_1d_s (spline->y_grid, spline->yBC, coefs+doffset, (intptr_t)ys,
245 coefs+coffset, (intptr_t)ys);
246 }
247}
248
249
252 BCtype_s xBC, BCtype_s yBC, BCtype_s zBC,
253 int num_splines)
254{
255 // Create new spline
256 multi_UBspline_3d_s* restrict spline = static_cast<multi_UBspline_3d_s*>(malloc(sizeof(multi_UBspline_3d_s)));
257 if (!spline) {
258 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_3d_s.\n");
259 abort();
260 }
261 spline->spcode = MULTI_U3D;
262 spline->tcode = SINGLE_REAL;
263 spline->xBC = xBC;
264 spline->yBC = yBC;
265 spline->zBC = zBC;
266 spline->num_splines = num_splines;
267 // Setup internal variables
268 int Mx = x_grid.num; int My = y_grid.num; int Mz = z_grid.num;
269 int Nx, Ny, Nz;
270
271 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)
272 Nx = Mx+3;
273 else
274 Nx = Mx+2;
275 x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
276 x_grid.delta_inv = 1.0/x_grid.delta;
277 spline->x_grid = x_grid;
278
279 if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)
280 Ny = My+3;
281 else
282 Ny = My+2;
283 y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
284 y_grid.delta_inv = 1.0/y_grid.delta;
285 spline->y_grid = y_grid;
286
287 if (zBC.lCode == PERIODIC || zBC.lCode == ANTIPERIODIC)
288 Nz = Mz+3;
289 else
290 Nz = Mz+2;
291 z_grid.delta = (z_grid.end - z_grid.start)/(double)(Nz-3);
292 z_grid.delta_inv = 1.0/z_grid.delta;
293 spline->z_grid = z_grid;
294
295 int N = num_splines;
296#ifdef HAVE_SSE
297 if (N % 4)
298 N += 4 - (N % 4);
299#endif
300
301 spline->x_stride = Ny*Nz*N;
302 spline->y_stride = Nz*N;
303 spline->z_stride = N;
304
305#ifndef HAVE_POSIX_MEMALIGN
306 spline->coefs = new float[Nx*Ny*Nz*N];
307#else
308 posix_memalign ((void**)&spline->coefs, 64,
309 (sizeof(float)*Nx*Ny*Nz*N));
310#endif
311#ifdef HAVE_SSE
313#endif
314 if (!spline->coefs) {
315 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_3d_s.\n");
316 abort();
317 }
318
319 return spline;
320}
321
322void
323set_multi_UBspline_3d_s (multi_UBspline_3d_s* spline, int num, float *data)
324{
325 int Mx = spline->x_grid.num;
326 int My = spline->y_grid.num;
327 int Mz = spline->z_grid.num;
328 int Nx, Ny, Nz;
329
330 if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)
331 Nx = Mx+3;
332 else
333 Nx = Mx+2;
334 if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)
335 Ny = My+3;
336 else
337 Ny = My+2;
338 if (spline->zBC.lCode == PERIODIC || spline->zBC.lCode == ANTIPERIODIC)
339 Nz = Mz+3;
340 else
341 Nz = Mz+2;
342
343 float *coefs = spline->coefs + num;
344
345 int zs = spline->z_stride;
346 // First, solve in the X-direction
347 for (int iy=0; iy<My; iy++)
348 for (int iz=0; iz<Mz; iz++) {
349 intptr_t doffset = iy*Mz+iz;
350 intptr_t coffset = (iy*Nz+iz)*zs;
351 find_coefs_1d_s (spline->x_grid, spline->xBC, data+doffset, (intptr_t)My*Mz,
352 coefs+coffset, (intptr_t)Ny*Nz*zs);
353 }
354
355 // Now, solve in the Y-direction
356 for (int ix=0; ix<Nx; ix++)
357 for (int iz=0; iz<Nz; iz++) {
358 intptr_t doffset = (ix*Ny*Nz + iz)*zs;
359 intptr_t coffset = (ix*Ny*Nz + iz)*zs;
360 find_coefs_1d_s (spline->y_grid, spline->yBC, coefs+doffset, (intptr_t)Nz*zs,
361 coefs+coffset, (intptr_t)Nz*zs);
362 }
363
364 // Now, solve in the Z-direction
365 for (int ix=0; ix<Nx; ix++)
366 for (int iy=0; iy<Ny; iy++) {
367 intptr_t doffset = ((ix*Ny+iy)*Nz)*zs;
368 intptr_t coffset = ((ix*Ny+iy)*Nz)*zs;
369 find_coefs_1d_s (spline->z_grid, spline->zBC, coefs+doffset,
370 (intptr_t)zs, coefs+coffset, (intptr_t)zs);
371 }
372}
373
374
380
381// On input, bands should be filled with:
382// row 0 : abcdInitial from boundary conditions
383// rows 1:M: basis functions in first 3 cols, data in last
384// row M+1 : abcdFinal from boundary conditions
385// cstride gives the stride between values in coefs.
386// On exit, coefs with contain interpolating B-spline coefs
388create_multi_UBspline_1d_c (Ugrid x_grid, BCtype_c xBC, int num_splines)
389{
390 // Create new spline
391 multi_UBspline_1d_c* restrict spline = static_cast<multi_UBspline_1d_c*>(malloc(sizeof(multi_UBspline_1d_c)));
392 if (!spline) {
393 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_1d_c.\n");
394 abort();
395 }
396 spline->spcode = MULTI_U1D;
397 spline->tcode = SINGLE_COMPLEX;
398 spline->xBC = xBC;
399 spline->num_splines = num_splines;
400 // Setup internal variables
401 int M = x_grid.num;
402 int N;
403
404 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC) {
405 x_grid.delta = (x_grid.end-x_grid.start)/(double)(x_grid.num);
406 N = M+3;
407 }
408 else {
409 x_grid.delta = (x_grid.end-x_grid.start)/(double)(x_grid.num-1);
410 N = M+2;
411 }
412
413 x_grid.delta_inv = 1.0/x_grid.delta;
414 spline->x_grid = x_grid;
415 spline->x_stride = num_splines;
416#ifndef HAVE_POSIX_MEMALIGN
417 spline->coefs = (complex_float*)malloc (2*sizeof(float)*N*num_splines);
418#else
419 posix_memalign ((void**)&spline->coefs, 64, 2*sizeof(float)*N*num_splines);
420#endif
421#ifdef HAVE_SSE
422 init_sse_data();
423#endif
424 if (!spline->coefs) {
425 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_1d_c.\n");
426 abort();
427 }
428
429 return spline;
430}
431
432void
434{
435 complex_float *coefs = spline->coefs + num;
436
437 BCtype_s xBC_r, xBC_i;
438 xBC_r.lCode = spline->xBC.lCode; xBC_r.rCode = spline->xBC.rCode;
439 xBC_r.lVal = spline->xBC.lVal_r; xBC_r.rVal = spline->xBC.rVal_r;
440 xBC_i.lCode = spline->xBC.lCode; xBC_i.rCode = spline->xBC.rCode;
441 xBC_i.lVal = spline->xBC.lVal_i; xBC_i.rVal = spline->xBC.rVal_i;
442
443 int xs = spline->x_stride;
444 // Real part
445 find_coefs_1d_s (spline->x_grid, xBC_r,
446 (float*)data, (intptr_t)2, (float*)coefs, (intptr_t)2*xs);
447 // Imaginarty part
448 find_coefs_1d_s (spline->x_grid, xBC_i,
449 ((float*)data)+1, (intptr_t)2, ((float*)coefs+1), (intptr_t)2*xs);
450}
451
452
453
456 BCtype_c xBC, BCtype_c yBC, int num_splines)
457{
458 // Create new spline
459 multi_UBspline_2d_c* restrict spline = static_cast<multi_UBspline_2d_c*>(malloc(sizeof(multi_UBspline_2d_c)));
460 if (!spline) {
461 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_2d_c.\n");
462 abort();
463 }
464 spline->spcode = MULTI_U2D;
465 spline->tcode = SINGLE_COMPLEX;
466 spline->xBC = xBC;
467 spline->yBC = yBC;
468 spline->num_splines = num_splines;
469
470 // Setup internal variables
471 int Mx = x_grid.num;
472 int My = y_grid.num;
473 int Nx, Ny;
474
475 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)
476 Nx = Mx+3;
477 else
478 Nx = Mx+2;
479 x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
480 x_grid.delta_inv = 1.0/x_grid.delta;
481 spline->x_grid = x_grid;
482
483 if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)
484 Ny = My+3;
485 else
486 Ny = My+2;
487 y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
488 y_grid.delta_inv = 1.0/y_grid.delta;
489 spline->y_grid = y_grid;
490
491 int N = num_splines;
492#ifdef HAVE_SSE
493 if (N % 2)
494 N++;
495#endif
496
497 spline->x_stride = Ny*N;
498 spline->y_stride = N;
499
500#ifndef HAVE_POSIX_MEMALIGN
501 spline->coefs = (complex_float*)malloc (2*sizeof(float)*Nx*Ny*N);
502 spline->lapl2 = (complex_float*)malloc (4*sizeof(float)*N);
503#else
504 posix_memalign ((void**)&spline->coefs, 64,
505 2*sizeof(float)*Nx*Ny*N);
506 posix_memalign ((void**)&spline->lapl2, 64,
507 4*sizeof(float)*N);
508#endif
509#ifdef HAVE_SSE
511#endif
512 if (!spline->coefs || !spline->lapl2) {
513 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_2d_c.\n");
514 abort();
515 }
516 return spline;
517}
518
519
520void
522{
523 // Setup internal variables
524 int Mx = spline->x_grid.num;
525 int My = spline->y_grid.num;
526 int Nx, Ny;
527
528 complex_float* coefs = spline->coefs + num;
529
530 if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)
531 Nx = Mx+3;
532 else
533 Nx = Mx+2;
534 if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)
535 Ny = My+3;
536 else
537 Ny = My+2;
538
539 BCtype_s xBC_r, xBC_i, yBC_r, yBC_i;
540 xBC_r.lCode = spline->xBC.lCode; xBC_r.rCode = spline->xBC.rCode;
541 xBC_r.lVal = spline->xBC.lVal_r; xBC_r.rVal = spline->xBC.rVal_r;
542 xBC_i.lCode = spline->xBC.lCode; xBC_i.rCode = spline->xBC.rCode;
543 xBC_i.lVal = spline->xBC.lVal_i; xBC_i.rVal = spline->xBC.rVal_i;
544 yBC_r.lCode = spline->yBC.lCode; yBC_r.rCode = spline->yBC.rCode;
545 yBC_r.lVal = spline->yBC.lVal_r; yBC_r.rVal = spline->yBC.rVal_r;
546 yBC_i.lCode = spline->yBC.lCode; yBC_i.rCode = spline->yBC.rCode;
547 yBC_i.lVal = spline->yBC.lVal_i; yBC_i.rVal = spline->yBC.rVal_i;
548
549 int ys = spline->y_stride;
550 // First, solve in the X-direction
551 for (int iy=0; iy<My; iy++) {
552 intptr_t doffset = (2*iy);
553 intptr_t coffset = (2*iy)*ys;
554 // Real part
555 find_coefs_1d_s (spline->x_grid, xBC_r, ((float*)data)+doffset, (intptr_t)2*My,
556 (float*)coefs+coffset, (intptr_t)2*Ny*ys);
557 // Imag part
558 find_coefs_1d_s (spline->x_grid, xBC_i, ((float*)data)+doffset+1, (intptr_t)2*My,
559 ((float*)coefs)+coffset+1, (intptr_t)2*Ny*ys);
560 }
561
562 // Now, solve in the Y-direction
563 for (int ix=0; ix<Nx; ix++) {
564 intptr_t doffset = (2*ix*Ny)*ys;
565 intptr_t coffset = (2*ix*Ny)*ys;
566 // Real part
567 find_coefs_1d_s (spline->y_grid, yBC_r, ((float*)coefs)+doffset,
568 (intptr_t)2*ys, ((float*)coefs)+coffset, (intptr_t)2*ys);
569 // Imag part
570 find_coefs_1d_s (spline->y_grid, yBC_i, ((float*)coefs)+doffset+1,
571 (intptr_t)2*ys, ((float*)coefs)+coffset+1, (intptr_t)2*ys);
572 }
573}
574
577 BCtype_c xBC, BCtype_c yBC, BCtype_c zBC,
578 int num_splines)
579{
580 // Create new spline
581 multi_UBspline_3d_c* restrict spline = static_cast<multi_UBspline_3d_c*>(malloc(sizeof(multi_UBspline_3d_c)));
582 if (!spline) {
583 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_3d_c.\n");
584 abort();
585 }
586 spline->spcode = MULTI_U3D;
587 spline->tcode = SINGLE_COMPLEX;
588 spline->xBC = xBC;
589 spline->yBC = yBC;
590 spline->zBC = zBC;
591 spline->num_splines = num_splines;
592
593 // Setup internal variables
594 int Mx = x_grid.num; int My = y_grid.num; int Mz = z_grid.num;
595 int Nx, Ny, Nz;
596
597 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)
598 Nx = Mx+3;
599 else
600 Nx = Mx+2;
601 x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
602 x_grid.delta_inv = 1.0/x_grid.delta;
603 spline->x_grid = x_grid;
604
605 if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)
606 Ny = My+3;
607 else
608 Ny = My+2;
609 y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
610 y_grid.delta_inv = 1.0/y_grid.delta;
611 spline->y_grid = y_grid;
612
613 if (zBC.lCode == PERIODIC || zBC.lCode == ANTIPERIODIC)
614 Nz = Mz+3;
615 else
616 Nz = Mz+2;
617 z_grid.delta = (z_grid.end - z_grid.start)/(double)(Nz-3);
618 z_grid.delta_inv = 1.0/z_grid.delta;
619 spline->z_grid = z_grid;
620
621 int N = spline->num_splines;
622#ifdef HAVE_SSE
623 if (N % 2)
624 N++;
625#endif
626
627 spline->x_stride = Ny*Nz*N;
628 spline->y_stride = Nz*N;
629 spline->z_stride = N;
630
631#ifndef HAVE_POSIX_MEMALIGN
632 spline->coefs = (complex_float*)malloc (2*sizeof(float)*Nx*Ny*Nz*N);
633 spline->lapl3 = (complex_float*)malloc (6*sizeof(float)*N);
634#else
635 posix_memalign ((void**)&spline->coefs, 64, 2*sizeof(float)*Nx*Ny*Nz*N);
636 posix_memalign ((void**)&spline->lapl3, 64, 6*sizeof(float)*N);
637#endif
638#ifdef HAVE_SSE
640#endif
641 if (!spline->coefs || !spline->lapl3) {
642 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_3d_c.\n");
643 abort();
644 }
645
646 return spline;
647}
648
649void
651{
652 // Setup internal variables
653 int Mx = spline->x_grid.num;
654 int My = spline->y_grid.num;
655 int Mz = spline->z_grid.num;
656 int Nx, Ny, Nz;
657
658 if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)
659 Nx = Mx+3;
660 else
661 Nx = Mx+2;
662 if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)
663 Ny = My+3;
664 else
665 Ny = My+2;
666 if (spline->zBC.lCode == PERIODIC || spline->zBC.lCode == ANTIPERIODIC)
667 Nz = Mz+3;
668 else
669 Nz = Mz+2;
670
671 BCtype_s xBC_r, xBC_i, yBC_r, yBC_i, zBC_r, zBC_i;
672 xBC_r.lCode = spline->xBC.lCode; xBC_r.rCode = spline->xBC.rCode;
673 xBC_r.lVal = spline->xBC.lVal_r; xBC_r.rVal = spline->xBC.rVal_r;
674 xBC_i.lCode = spline->xBC.lCode; xBC_i.rCode = spline->xBC.rCode;
675 xBC_i.lVal = spline->xBC.lVal_i; xBC_i.rVal = spline->xBC.rVal_i;
676 yBC_r.lCode = spline->yBC.lCode; yBC_r.rCode = spline->yBC.rCode;
677 yBC_r.lVal = spline->yBC.lVal_r; yBC_r.rVal = spline->yBC.rVal_r;
678 yBC_i.lCode = spline->yBC.lCode; yBC_i.rCode = spline->yBC.rCode;
679 yBC_i.lVal = spline->yBC.lVal_i; yBC_i.rVal = spline->yBC.rVal_i;
680 zBC_r.lCode = spline->zBC.lCode; zBC_r.rCode = spline->zBC.rCode;
681 zBC_r.lVal = spline->zBC.lVal_r; zBC_r.rVal = spline->zBC.rVal_r;
682 zBC_i.lCode = spline->zBC.lCode; zBC_i.rCode = spline->zBC.rCode;
683 zBC_i.lVal = spline->zBC.lVal_i; zBC_i.rVal = spline->zBC.rVal_i;
684
685 complex_float *coefs = spline->coefs + num;
686 int zs = spline->z_stride;
687 // First, solve in the X-direction
688 for (int iy=0; iy<My; iy++)
689 for (int iz=0; iz<Mz; iz++) {
690 intptr_t doffset = 2*(iy*Mz+iz);
691 intptr_t coffset = 2*(iy*Nz+iz)*zs;
692 // Real part
693 find_coefs_1d_s (spline->x_grid, xBC_r,
694 ((float*)data)+doffset, (intptr_t)2*My*Mz,
695 ((float*)coefs)+coffset, (intptr_t)2*Ny*Nz*zs);
696 // Imag part
697 find_coefs_1d_s (spline->x_grid, xBC_i,
698 ((float*)data)+doffset+1, (intptr_t)2*My*Mz,
699 ((float*)coefs)+coffset+1, (intptr_t)2*Ny*Nz*zs);
700 }
701
702 // Now, solve in the Y-direction
703 for (int ix=0; ix<Nx; ix++)
704 for (int iz=0; iz<Nz; iz++) {
705 intptr_t doffset = 2*(ix*Ny*Nz + iz)*zs;
706 intptr_t coffset = 2*(ix*Ny*Nz + iz)*zs;
707 // Real part
708 find_coefs_1d_s (spline->y_grid, yBC_r,
709 ((float*)coefs)+doffset, (intptr_t)2*Nz*zs,
710 ((float*)coefs)+coffset, (intptr_t)2*Nz*zs);
711 // Imag part
712 find_coefs_1d_s (spline->y_grid, yBC_i,
713 ((float*)coefs)+doffset+1, (intptr_t)2*Nz*zs,
714 ((float*)coefs)+coffset+1, (intptr_t)2*Nz*zs);
715 }
716
717 // Now, solve in the Z-direction
718 for (int ix=0; ix<Nx; ix++)
719 for (int iy=0; iy<Ny; iy++) {
720 intptr_t doffset = 2*((ix*Ny+iy)*Nz)*zs;
721 intptr_t coffset = 2*((ix*Ny+iy)*Nz)*zs;
722 // Real part
723 find_coefs_1d_s (spline->z_grid, zBC_r,
724 ((float*)coefs)+doffset, (intptr_t)2*zs,
725 ((float*)coefs)+coffset, (intptr_t)2*zs);
726 // Imag part
727 find_coefs_1d_s (spline->z_grid, zBC_i,
728 ((float*)coefs)+doffset+1, (intptr_t)2*zs,
729 ((float*)coefs)+coffset+1, (intptr_t)2*zs);
730 }
731}
732
733
739
740// On input, bands should be filled with:
741// row 0 : abcdInitial from boundary conditions
742// rows 1:M: basis functions in first 3 cols, data in last
743// row M+1 : abcdFinal from boundary conditions
744// cstride gives the stride between values in coefs.
745// On exit, coefs with contain interpolating B-spline coefs
746void
747solve_deriv_interp_1d_d (double bands[], double coefs[],
748 int M, int cstride);
749
750// On input, bands should be filled with:
751// row 0 : abcdInitial from boundary conditions
752// rows 1:M: basis functions in first 3 cols, data in last
753// row M+1 : abcdFinal from boundary conditions
754// cstride gives the stride between values in coefs.
755// On exit, coefs with contain interpolating B-spline coefs
756void
757solve_periodic_interp_1d_d (double bands[], double coefs[],
758 int M, int cstride);
759
760void
762 double *data, intptr_t dstride,
763 double *coefs, intptr_t cstride);
764
766create_multi_UBspline_1d_d (Ugrid x_grid, BCtype_d xBC, int num_splines)
767{
768 // Create new spline
769 multi_UBspline_1d_d* restrict spline = static_cast<multi_UBspline_1d_d*>(malloc(sizeof(multi_UBspline_1d_d)));
770 if (!spline) {
771 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_1d_d.\n");
772 abort();
773 }
774 spline->spcode = MULTI_U1D;
775 spline->tcode = DOUBLE_REAL;
776 spline->xBC = xBC;
777 spline->num_splines = num_splines;
778
779 // Setup internal variables
780 int Mx = x_grid.num;
781 int Nx;
782
783 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC) {
784 x_grid.delta = (x_grid.end-x_grid.start)/(double)(x_grid.num);
785 Nx = Mx+3;
786 }
787 else {
788 x_grid.delta = (x_grid.end-x_grid.start)/(double)(x_grid.num-1);
789 Nx = Mx+2;
790 }
791
792 x_grid.delta_inv = 1.0/x_grid.delta;
793 spline->x_grid = x_grid;
794
795 int N = num_splines;
796#ifdef HAVE_SSE2
797 // We must pad to keep data aligned for SSE operations
798 if (N & 1)
799 N++;
800#endif
801 spline->x_stride = N;
802
803#ifndef HAVE_POSIX_MEMALIGN
804 spline->coefs = (double*)malloc (sizeof(double)*Nx*N);
805
806#else
807 posix_memalign ((void**)&spline->coefs, 64, sizeof(double)*Nx*N);
808#endif
809#ifdef HAVE_SSE2
811#endif
812 if (!spline->coefs) {
813 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_1d_d.\n");
814 abort();
815 }
816
817 return spline;
818}
819
820void
821set_multi_UBspline_1d_d (multi_UBspline_1d_d* spline, int num, double *data)
822{
823 double *coefs = spline->coefs + num;
824 int xs = spline->x_stride;
825 find_coefs_1d_d (spline->x_grid, spline->xBC, data, 1, coefs, xs);
826}
827
828void
829set_multi_UBspline_1d_d_BC (multi_UBspline_1d_d* spline, int num, double *data,
830 BCtype_d xBC)
831{
832 double *coefs = spline->coefs + num;
833 int xs = spline->x_stride;
834 find_coefs_1d_d (spline->x_grid, xBC, data, 1, coefs, xs);
835}
836
837
840 BCtype_d xBC, BCtype_d yBC, int num_splines)
841{
842 // Create new spline
843 multi_UBspline_2d_d* restrict spline = static_cast<multi_UBspline_2d_d*>(malloc(sizeof(multi_UBspline_2d_d)));
844 if (!spline) {
845 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_2d_d.\n");
846 abort();
847 }
848 spline->spcode = MULTI_U2D;
849 spline->tcode = DOUBLE_REAL;
850 spline->xBC = xBC;
851 spline->yBC = yBC;
852 spline->num_splines = num_splines;
853
854 // Setup internal variables
855 int Mx = x_grid.num;
856 int My = y_grid.num;
857 int Nx, Ny;
858
859 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)
860 Nx = Mx+3;
861 else
862 Nx = Mx+2;
863 x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
864 x_grid.delta_inv = 1.0/x_grid.delta;
865 spline->x_grid = x_grid;
866
867 if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)
868 Ny = My+3;
869 else
870 Ny = My+2;
871 y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
872 y_grid.delta_inv = 1.0/y_grid.delta;
873 spline->y_grid = y_grid;
874
875 int N = num_splines;
876
877#ifdef HAVE_SSE2
878 // We must pad to keep data align for SSE operations
879 if (num_splines & 1)
880 N++;
881#endif
882 spline->x_stride = Ny*N;
883 spline->y_stride = N;
884
885#ifndef HAVE_POSIX_MEMALIGN
886 spline->coefs = (double*)malloc (sizeof(double)*Nx*Ny*N);
887#else
888 posix_memalign ((void**)&spline->coefs, 64, (sizeof(double)*Nx*Ny*N));
889#endif
890#ifdef HAVE_SSE2
892#endif
893 if (!spline->coefs) {
894 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_2d_d.\n");
895 abort();
896 }
897
898 return spline;
899}
900
901void
902set_multi_UBspline_2d_d (multi_UBspline_2d_d* spline, int num, double *data)
903{
904 int Mx = spline->x_grid.num;
905 int My = spline->y_grid.num;
906 int Nx, Ny;
907 double *coefs = spline->coefs + num;
908
909 if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)
910 Nx = Mx+3;
911 else
912 Nx = Mx+2;
913
914 if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)
915 Ny = My+3;
916 else
917 Ny = My+2;
918
919 int ys = spline->y_stride;
920 // First, solve in the X-direction
921 for (int iy=0; iy<My; iy++) {
922 intptr_t doffset = iy;
923 intptr_t coffset = iy*ys;
924 find_coefs_1d_d (spline->x_grid, spline->xBC,
925 data+doffset, (intptr_t)My,
926 coefs+coffset, (intptr_t)Ny*ys);
927 }
928
929 // Now, solve in the Y-direction
930 for (int ix=0; ix<Nx; ix++) {
931 intptr_t doffset = ix*Ny*ys;
932 intptr_t coffset = ix*Ny*ys;
933 find_coefs_1d_d (spline->y_grid, spline->yBC,
934 coefs+doffset, (intptr_t)ys,
935 coefs+coffset, (intptr_t)ys);
936 }
937}
938
939
942 BCtype_d xBC, BCtype_d yBC, BCtype_d zBC,
943 int num_splines)
944{
945 // Create new spline
947#ifdef HAVE_POSIX_MEMALIGN
948 posix_memalign ((void**)&spline, 64, (size_t)sizeof(multi_UBspline_3d_d));
949#else
950 spline = static_cast<multi_UBspline_3d_d*>(malloc(sizeof(multi_UBspline_3d_d)));
951#endif
952 if (!spline) {
953 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_3d_d.\n");
954 abort();
955 }
956 spline->spcode = MULTI_U3D;
957 spline->tcode = DOUBLE_REAL;
958 spline->xBC = xBC;
959 spline->yBC = yBC;
960 spline->zBC = zBC;
961 spline->num_splines = num_splines;
962
963 // Setup internal variables
964 int Mx = x_grid.num; int My = y_grid.num; int Mz = z_grid.num;
965 int Nx, Ny, Nz;
966
967 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)
968 Nx = Mx+3;
969 else
970 Nx = Mx+2;
971 x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
972 x_grid.delta_inv = 1.0/x_grid.delta;
973 spline->x_grid = x_grid;
974
975 if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)
976 Ny = My+3;
977 else
978 Ny = My+2;
979 y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
980 y_grid.delta_inv = 1.0/y_grid.delta;
981 spline->y_grid = y_grid;
982
983 if (zBC.lCode == PERIODIC || zBC.lCode == ANTIPERIODIC)
984 Nz = Mz+3;
985 else
986 Nz = Mz+2;
987 z_grid.delta = (z_grid.end - z_grid.start)/(double)(Nz-3);
988 z_grid.delta_inv = 1.0/z_grid.delta;
989 spline->z_grid = z_grid;
990
991
992 int N = num_splines;
993#if defined HAVE_SSE2 || defined HAVE_VSX
994 // We must pad to keep data align for SSE operations
995 if (N & 1)
996 N++;
997#endif
998
999 spline->x_stride = Ny*Nz*N;
1000 spline->y_stride = Nz*N;
1001 spline->z_stride = N;
1002
1003#ifdef HAVE_POSIX_MEMALIGN
1004 posix_memalign ((void**)&spline->coefs, 64, ((size_t)sizeof(double)*Nx*Ny*Nz*N));
1005#else
1006 spline->coefs = new double[Nx*Ny*Nz*N];
1007#endif
1008#ifdef HAVE_SSE2
1009 init_sse_data();
1010#endif
1011 if (!spline->coefs) {
1012 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_3d_d.\n");
1013 abort();
1014 }
1015
1016 return spline;
1017}
1018
1019void
1020set_multi_UBspline_3d_d (multi_UBspline_3d_d* spline, int num, double *data)
1021{
1022 int Mx = spline->x_grid.num;
1023 int My = spline->y_grid.num;
1024 int Mz = spline->z_grid.num;
1025 int Nx, Ny, Nz;
1026
1027 if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)
1028 Nx = Mx+3;
1029 else
1030 Nx = Mx+2;
1031 if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)
1032 Ny = My+3;
1033 else
1034 Ny = My+2;
1035 if (spline->zBC.lCode == PERIODIC || spline->zBC.lCode == ANTIPERIODIC)
1036 Nz = Mz+3;
1037 else
1038 Nz = Mz+2;
1039
1040 double *coefs = spline->coefs + num;
1041 intptr_t zs = spline->z_stride;
1042
1043 // First, solve in the X-direction
1044 for (int iy=0; iy<My; iy++)
1045 for (int iz=0; iz<Mz; iz++) {
1046 intptr_t doffset = iy*Mz+iz;
1047 intptr_t coffset = (iy*Nz+iz)*zs;
1048 find_coefs_1d_d (spline->x_grid, spline->xBC,
1049 data+doffset, (intptr_t)My*Mz,
1050 coefs+coffset, (intptr_t)Ny*Nz*zs);
1051 }
1052
1053 // Now, solve in the Y-direction
1054 for (int ix=0; ix<Nx; ix++)
1055 for (int iz=0; iz<Nz; iz++) {
1056 intptr_t doffset = (ix*Ny*Nz + iz)*zs;
1057 intptr_t coffset = (ix*Ny*Nz + iz)*zs;
1058 find_coefs_1d_d (spline->y_grid, spline->yBC,
1059 coefs+doffset, (intptr_t)Nz*zs,
1060 coefs+coffset, (intptr_t)Nz*zs);
1061 }
1062
1063 // Now, solve in the Z-direction
1064 for (int ix=0; ix<Nx; ix++)
1065 for (int iy=0; iy<Ny; iy++) {
1066 intptr_t doffset = (ix*Ny+iy)*Nz*zs;
1067 intptr_t coffset = (ix*Ny+iy)*Nz*zs;
1068 find_coefs_1d_d (spline->z_grid, spline->zBC,
1069 coefs+doffset, (intptr_t)zs,
1070 coefs+coffset, (intptr_t)zs);
1071 }
1072}
1073
1074
1080
1081// On input, bands should be filled with:
1082// row 0 : abcdInitial from boundary conditions
1083// rows 1:M: basis functions in first 3 cols, data in last
1084// row M+1 : abcdFinal from boundary conditions
1085// cstride gives the stride between values in coefs.
1086// On exit, coefs with contain interpolating B-spline coefs
1087
1088
1090create_multi_UBspline_1d_z (Ugrid x_grid, BCtype_z xBC, int num_splines)
1091{
1092 // Create new spline
1093 multi_UBspline_1d_z* restrict spline = static_cast<multi_UBspline_1d_z*>(malloc(sizeof(multi_UBspline_1d_z)));
1094 if (!spline) {
1095 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_1d_z.\n");
1096 abort();
1097 }
1098 spline->spcode = MULTI_U1D;
1099 spline->tcode = DOUBLE_COMPLEX;
1100 spline->xBC = xBC;
1101 spline->num_splines = num_splines;
1102
1103 // Setup internal variables
1104 int Mx = x_grid.num;
1105 int Nx;
1106
1107 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC) {
1108 x_grid.delta = (x_grid.end-x_grid.start)/(double)(x_grid.num);
1109 Nx = Mx+3;
1110 }
1111 else {
1112 x_grid.delta = (x_grid.end-x_grid.start)/(double)(x_grid.num-1);
1113 Nx = Mx+2;
1114 }
1115
1116 x_grid.delta_inv = 1.0/x_grid.delta;
1117 spline->x_grid = x_grid;
1118 spline->x_stride = num_splines;
1119
1120#ifndef HAVE_POSIX_MEMALIGN
1121 spline->coefs = (complex_double*)malloc (2*sizeof(double)*Nx*num_splines);
1122#else
1123 posix_memalign ((void**)&spline->coefs, 64, 2*sizeof(double)*Nx*num_splines);
1124#endif
1125#ifdef HAVE_SSE2
1126 init_sse_data();
1127#endif
1128 if (!spline->coefs) {
1129 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_1d_z.\n");
1130 abort();
1131 }
1132
1133 return spline;
1134}
1135
1136void
1138{
1139// int Mx = spline->x_grid.num;
1140// Set but not used
1141// int Nx;
1142
1143 complex_double *coefs = spline->coefs + num;
1144
1145// if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)
1146// Nx = Mx+3;
1147// else
1148// Nx = Mx+2;
1149
1150 BCtype_d xBC_r, xBC_i;
1151 xBC_r.lCode = spline->xBC.lCode; xBC_r.rCode = spline->xBC.rCode;
1152 xBC_r.lVal = spline->xBC.lVal_r; xBC_r.rVal = spline->xBC.rVal_r;
1153 xBC_i.lCode = spline->xBC.lCode; xBC_i.rCode = spline->xBC.rCode;
1154 xBC_i.lVal = spline->xBC.lVal_i; xBC_i.rVal = spline->xBC.rVal_i;
1155 int xs = spline->x_stride;
1156 // Real part
1157 find_coefs_1d_d (spline->x_grid, xBC_r,
1158 (double*)data, (intptr_t)2,
1159 ((double*)coefs), (intptr_t)2*xs);
1160 // Imaginary part
1161 find_coefs_1d_d (spline->x_grid, xBC_i,
1162 ((double*)data)+1, (intptr_t)2,
1163 ((double*)coefs)+1, (intptr_t)2*xs);
1164
1165}
1166
1167void
1169 complex_double *data, BCtype_z xBC)
1170{
1171// int Mx = spline->x_grid.num;
1172// int Nx;
1173
1174 complex_double *coefs = spline->coefs + num;
1175
1176// if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)
1177// Nx = Mx+3;
1178// else
1179// Nx = Mx+2;
1180
1181 BCtype_d xBC_r, xBC_i;
1182 xBC_r.lCode = xBC.lCode; xBC_r.rCode = xBC.rCode;
1183 xBC_r.lVal = xBC.lVal_r; xBC_r.rVal = xBC.rVal_r;
1184 xBC_i.lCode = xBC.lCode; xBC_i.rCode = xBC.rCode;
1185 xBC_i.lVal = xBC.lVal_i; xBC_i.rVal = xBC.rVal_i;
1186 int xs = spline->x_stride;
1187 // Real part
1188 find_coefs_1d_d (spline->x_grid, xBC_r,
1189 (double*)data, (intptr_t)2,
1190 ((double*)coefs), (intptr_t)2*xs);
1191 // Imaginary part
1192 find_coefs_1d_d (spline->x_grid, xBC_i,
1193 ((double*)data)+1, (intptr_t)2,
1194 ((double*)coefs)+1, (intptr_t)2*xs);
1195}
1196
1197
1200 BCtype_z xBC, BCtype_z yBC, int num_splines)
1201{
1202 // Create new spline
1203 multi_UBspline_2d_z* restrict spline = static_cast<multi_UBspline_2d_z*>(malloc(sizeof(multi_UBspline_2d_z)));
1204 if (!spline) {
1205 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_2d_z.\n");
1206 abort();
1207 }
1208 spline->spcode = MULTI_U2D;
1209 spline->tcode = DOUBLE_COMPLEX;
1210 spline->xBC = xBC;
1211 spline->yBC = yBC;
1212 spline->num_splines = num_splines;
1213
1214 // Setup internal variables
1215 int Mx = x_grid.num;
1216 int My = y_grid.num;
1217 int Nx, Ny;
1218
1219 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)
1220 Nx = Mx+3;
1221 else
1222 Nx = Mx+2;
1223 x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
1224 x_grid.delta_inv = 1.0/x_grid.delta;
1225 spline->x_grid = x_grid;
1226
1227 if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)
1228 Ny = My+3;
1229 else
1230 Ny = My+2;
1231 y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
1232 y_grid.delta_inv = 1.0/y_grid.delta;
1233 spline->y_grid = y_grid;
1234 spline->x_stride = Ny*num_splines;
1235 spline->y_stride = num_splines;
1236
1237#ifndef HAVE_POSIX_MEMALIGN
1238 spline->coefs = (complex_double*)malloc (2*sizeof(double)*Nx*Ny*num_splines);
1239 spline->lapl2 = (complex_double*)malloc (4*sizeof(double)*num_splines);
1240#else
1241 posix_memalign ((void**)&spline->coefs, 64, 2*sizeof(double)*Nx*Ny*num_splines);
1242 posix_memalign ((void**)&spline->lapl2, 64, 4*sizeof(double)*num_splines);
1243#endif
1244#ifdef HAVE_SSE2
1245 init_sse_data();
1246#endif
1247 if (!spline->coefs || !spline->lapl2) {
1248 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_2d_z.\n");
1249 abort();
1250 }
1251
1252 return spline;
1253}
1254
1255
1256void
1258 complex_double *data)
1259{
1260 int Mx = spline->x_grid.num;
1261 int My = spline->y_grid.num;
1262 int Nx, Ny;
1263
1264 if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)
1265 Nx = Mx+3;
1266 else
1267 Nx = Mx+2;
1268 if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)
1269 Ny = My+3;
1270 else
1271 Ny = My+2;
1272
1273 BCtype_d xBC_r, xBC_i, yBC_r, yBC_i;
1274 xBC_r.lCode = spline->xBC.lCode; xBC_r.rCode = spline->xBC.rCode;
1275 xBC_r.lVal = spline->xBC.lVal_r; xBC_r.rVal = spline->xBC.rVal_r;
1276 xBC_i.lCode = spline->xBC.lCode; xBC_i.rCode = spline->xBC.rCode;
1277 xBC_i.lVal = spline->xBC.lVal_i; xBC_i.rVal = spline->xBC.rVal_i;
1278 yBC_r.lCode = spline->yBC.lCode; yBC_r.rCode = spline->yBC.rCode;
1279 yBC_r.lVal = spline->yBC.lVal_r; yBC_r.rVal = spline->yBC.rVal_r;
1280 yBC_i.lCode = spline->yBC.lCode; yBC_i.rCode = spline->yBC.rCode;
1281 yBC_i.lVal = spline->yBC.lVal_i; yBC_i.rVal = spline->yBC.rVal_i;
1282
1283 complex_double *coefs = spline->coefs + num;
1284 int ys = spline->y_stride;
1285
1286 // First, solve in the X-direction
1287 for (int iy=0; iy<My; iy++) {
1288 intptr_t doffset = 2*iy;
1289 intptr_t coffset = 2*iy*ys;
1290 // Real part
1291 find_coefs_1d_d (spline->x_grid, xBC_r,
1292 ((double*)data+doffset), (intptr_t)2*My,
1293 (double*)coefs+coffset, (intptr_t)2*Ny*ys);
1294 // Imag part
1295 find_coefs_1d_d (spline->x_grid, xBC_i,
1296 ((double*)data)+doffset+1, (intptr_t)2*My,
1297 ((double*)coefs)+coffset+1, (intptr_t)2*Ny*ys);
1298 }
1299
1300 // Now, solve in the Y-direction
1301 for (int ix=0; ix<Nx; ix++) {
1302 intptr_t doffset = 2*ix*Ny*ys;
1303 intptr_t coffset = 2*ix*Ny*ys;
1304 // Real part
1305 find_coefs_1d_d (spline->y_grid, yBC_r,
1306 ((double*)coefs)+doffset, (intptr_t)2*ys,
1307 (double*)coefs+coffset, (intptr_t)2*ys);
1308 // Imag part
1309 find_coefs_1d_d (spline->y_grid, yBC_i,
1310 (double*)coefs+doffset+1, (intptr_t)2*ys,
1311 ((double*)coefs)+coffset+1, (intptr_t)2*ys);
1312 }
1313}
1314
1315
1316
1319 BCtype_z xBC, BCtype_z yBC, BCtype_z zBC,
1320 int num_splines)
1321{
1322 // Create new spline
1323 multi_UBspline_3d_z* restrict spline = static_cast<multi_UBspline_3d_z*>(malloc(sizeof(multi_UBspline_3d_z)));
1324 if (!spline) {
1325 fprintf (stderr, "Out of memory allocating spline in create_multi_UBspline_3d_z.\n");
1326 abort();
1327 }
1328 spline->spcode = MULTI_U3D;
1329 spline->tcode = DOUBLE_COMPLEX;
1330 spline->xBC = xBC;
1331 spline->yBC = yBC;
1332 spline->zBC = zBC;
1333 spline->num_splines = num_splines;
1334
1335 // Setup internal variables
1336 int Mx = x_grid.num; int My = y_grid.num; int Mz = z_grid.num;
1337 int Nx, Ny, Nz;
1338
1339 if (xBC.lCode == PERIODIC || xBC.lCode == ANTIPERIODIC)
1340 Nx = Mx+3;
1341 else
1342 Nx = Mx+2;
1343 x_grid.delta = (x_grid.end - x_grid.start)/(double)(Nx-3);
1344 x_grid.delta_inv = 1.0/x_grid.delta;
1345 spline->x_grid = x_grid;
1346
1347 if (yBC.lCode == PERIODIC || yBC.lCode == ANTIPERIODIC)
1348 Ny = My+3;
1349 else
1350 Ny = My+2;
1351 y_grid.delta = (y_grid.end - y_grid.start)/(double)(Ny-3);
1352 y_grid.delta_inv = 1.0/y_grid.delta;
1353 spline->y_grid = y_grid;
1354
1355 if (zBC.lCode == PERIODIC || zBC.lCode == ANTIPERIODIC)
1356 Nz = Mz+3;
1357 else
1358 Nz = Mz+2;
1359 z_grid.delta = (z_grid.end - z_grid.start)/(double)(Nz-3);
1360 z_grid.delta_inv = 1.0/z_grid.delta;
1361 spline->z_grid = z_grid;
1362
1363 int N = num_splines;
1364#ifdef HAVE_SSE2
1365 if (N & 3)
1366 N += 4-(N & 3);
1367#endif
1368
1369 spline->x_stride = (intptr_t)Ny*(intptr_t)Nz*N;
1370 spline->y_stride = Nz*N;
1371 spline->z_stride = N;
1372
1373#ifndef HAVE_POSIX_MEMALIGN
1374 spline->coefs = new complex_double[Nx*Ny*Nz*N];
1375 spline->lapl3 = new complex_double[3*N];
1376#else
1377 posix_memalign ((void**)&spline->coefs, 64, (size_t)2*sizeof(double)*Nx*Ny*Nz*N);
1378 posix_memalign ((void**)&spline->lapl3, 64, 6*sizeof(double)*N);
1379#endif
1380
1381#ifdef HAVE_SSE2
1382 init_sse_data();
1383#endif
1384 if (!spline->coefs || !spline->lapl3) {
1385 fprintf (stderr, "Out of memory allocating spline coefficients in create_multi_UBspline_3d_z.\n");
1386 abort();
1387 }
1388
1389 return spline;
1390}
1391
1392void
1394{
1395 // Setup internal variables
1396 int Mx = spline->x_grid.num;
1397 int My = spline->y_grid.num;
1398 int Mz = spline->z_grid.num;
1399 int Nx, Ny, Nz;
1400
1401 if (spline->xBC.lCode == PERIODIC || spline->xBC.lCode == ANTIPERIODIC)
1402 Nx = Mx+3;
1403 else
1404 Nx = Mx+2;
1405 if (spline->yBC.lCode == PERIODIC || spline->yBC.lCode == ANTIPERIODIC)
1406 Ny = My+3;
1407 else
1408 Ny = My+2;
1409 if (spline->zBC.lCode == PERIODIC || spline->zBC.lCode == ANTIPERIODIC)
1410 Nz = Mz+3;
1411 else
1412 Nz = Mz+2;
1413
1414 BCtype_d xBC_r, xBC_i, yBC_r, yBC_i, zBC_r, zBC_i;
1415 xBC_r.lCode = spline->xBC.lCode; xBC_r.rCode = spline->xBC.rCode;
1416 xBC_r.lVal = spline->xBC.lVal_r; xBC_r.rVal = spline->xBC.rVal_r;
1417 xBC_i.lCode = spline->xBC.lCode; xBC_i.rCode = spline->xBC.rCode;
1418 xBC_i.lVal = spline->xBC.lVal_i; xBC_i.rVal = spline->xBC.rVal_i;
1419 yBC_r.lCode = spline->yBC.lCode; yBC_r.rCode = spline->yBC.rCode;
1420 yBC_r.lVal = spline->yBC.lVal_r; yBC_r.rVal = spline->yBC.rVal_r;
1421 yBC_i.lCode = spline->yBC.lCode; yBC_i.rCode = spline->yBC.rCode;
1422 yBC_i.lVal = spline->yBC.lVal_i; yBC_i.rVal = spline->yBC.rVal_i;
1423 zBC_r.lCode = spline->zBC.lCode; zBC_r.rCode = spline->zBC.rCode;
1424 zBC_r.lVal = spline->zBC.lVal_r; zBC_r.rVal = spline->zBC.rVal_r;
1425 zBC_i.lCode = spline->zBC.lCode; zBC_i.rCode = spline->zBC.rCode;
1426 zBC_i.lVal = spline->zBC.lVal_i; zBC_i.rVal = spline->zBC.rVal_i;
1427
1428 complex_double *coefs = spline->coefs + num;
1429
1430// unused variable
1431// int N = spline->num_splines;
1432 int zs = spline->z_stride;
1433 // First, solve in the X-direction
1434 for (int iy=0; iy<My; iy++)
1435 for (int iz=0; iz<Mz; iz++) {
1436 intptr_t doffset = 2*(iy*Mz+iz);
1437 intptr_t coffset = 2*(iy*Nz+iz)*zs;
1438 // Real part
1439 find_coefs_1d_d (spline->x_grid, xBC_r,
1440 ((double*)data)+doffset, (intptr_t)2*My*Mz,
1441 ((double*)coefs)+coffset, (intptr_t)2*Ny*Nz*zs);
1442 // Imag part
1443 find_coefs_1d_d (spline->x_grid, xBC_i,
1444 ((double*)data)+doffset+1, (intptr_t)2*My*Mz,
1445 ((double*)coefs)+coffset+1, (intptr_t)2*Ny*Nz*zs);
1446 }
1447
1448 // Now, solve in the Y-direction
1449 for (int ix=0; ix<Nx; ix++)
1450 for (int iz=0; iz<Nz; iz++) {
1451 intptr_t doffset = 2*(ix*Ny*Nz + iz)*zs;
1452 intptr_t coffset = 2*(ix*Ny*Nz + iz)*zs;
1453 // Real part
1454 find_coefs_1d_d (spline->y_grid, yBC_r,
1455 ((double*)coefs)+doffset, (intptr_t)2*Nz*zs,
1456 ((double*)coefs)+coffset, (intptr_t)2*Nz*zs);
1457 // Imag part
1458 find_coefs_1d_d (spline->y_grid, yBC_i,
1459 ((double*)coefs)+doffset+1, (intptr_t)2*Nz*zs,
1460 ((double*)coefs)+coffset+1, (intptr_t)2*Nz*zs);
1461 }
1462
1463 // Now, solve in the Z-direction
1464 for (int ix=0; ix<Nx; ix++)
1465 for (int iy=0; iy<Ny; iy++) {
1466 intptr_t doffset = 2*((ix*Ny+iy)*Nz)*zs;
1467 intptr_t coffset = 2*((ix*Ny+iy)*Nz)*zs;
1468 // Real part
1469 find_coefs_1d_d (spline->z_grid, zBC_r,
1470 ((double*)coefs)+doffset, (intptr_t)2*zs,
1471 ((double*)coefs)+coffset, (intptr_t)2*zs);
1472 // Imag part
1473 find_coefs_1d_d (spline->z_grid, zBC_i,
1474 ((double*)coefs)+doffset+1, (intptr_t)2*zs,
1475 ((double*)coefs)+coffset+1, (intptr_t)2*zs);
1476 }
1477}
1478
1479
1480void
1482{
1483 free(spline->coefs);
1484 free(spline);
1485}
complex float complex_float
@ DOUBLE_REAL
@ SINGLE_REAL
@ DOUBLE_COMPLEX
@ SINGLE_COMPLEX
complex double complex_double
@ ANTIPERIODIC
@ PERIODIC
@ MULTI_U1D
@ MULTI_U2D
@ MULTI_U3D
#define restrict
void set_multi_UBspline_2d_z(multi_UBspline_2d_z *spline, int num, complex_double *data)
void find_coefs_1d_s(Ugrid grid, BCtype_s bc, float *data, intptr_t dstride, float *coefs, intptr_t cstride)
void set_multi_UBspline_3d_d(multi_UBspline_3d_d *spline, int num, double *data)
void set_multi_UBspline_1d_d(multi_UBspline_1d_d *spline, int num, double *data)
multi_UBspline_1d_c * create_multi_UBspline_1d_c(Ugrid x_grid, BCtype_c xBC, int num_splines)
multi_UBspline_2d_d * create_multi_UBspline_2d_d(Ugrid x_grid, Ugrid y_grid, BCtype_d xBC, BCtype_d yBC, int num_splines)
void set_multi_UBspline_2d_d(multi_UBspline_2d_d *spline, int num, double *data)
void init_sse_data()
void set_multi_UBspline_1d_z_BC(multi_UBspline_1d_z *spline, int num, complex_double *data, BCtype_z xBC)
void set_multi_UBspline_2d_c(multi_UBspline_2d_c *spline, int num, complex_float *data)
void set_multi_UBspline_3d_s(multi_UBspline_3d_s *spline, int num, float *data)
multi_UBspline_3d_z * create_multi_UBspline_3d_z(Ugrid x_grid, Ugrid y_grid, Ugrid z_grid, BCtype_z xBC, BCtype_z yBC, BCtype_z zBC, int num_splines)
void solve_periodic_interp_1d_d(double bands[], double coefs[], int M, int cstride)
multi_UBspline_3d_d * create_multi_UBspline_3d_d(Ugrid x_grid, Ugrid y_grid, Ugrid z_grid, BCtype_d xBC, BCtype_d yBC, BCtype_d zBC, int num_splines)
void set_multi_UBspline_1d_c(multi_UBspline_1d_c *spline, int num, complex_float *data)
void solve_deriv_interp_1d_d(double bands[], double coefs[], int M, int cstride)
void solve_deriv_interp_1d_s(float bands[], float coefs[], int M, int cstride)
void set_multi_UBspline_3d_c(multi_UBspline_3d_c *spline, int num, complex_float *data)
multi_UBspline_2d_s * create_multi_UBspline_2d_s(Ugrid x_grid, Ugrid y_grid, BCtype_s xBC, BCtype_s yBC, int num_splines)
void set_multi_UBspline_3d_z(multi_UBspline_3d_z *spline, int num, complex_double *data)
multi_UBspline_2d_z * create_multi_UBspline_2d_z(Ugrid x_grid, Ugrid y_grid, BCtype_z xBC, BCtype_z yBC, int num_splines)
multi_UBspline_1d_z * create_multi_UBspline_1d_z(Ugrid x_grid, BCtype_z xBC, int num_splines)
void solve_periodic_interp_1d_s(float bands[], float coefs[], int M, int cstride)
void find_coefs_1d_d(Ugrid grid, BCtype_d bc, double *data, intptr_t dstride, double *coefs, intptr_t cstride)
multi_UBspline_1d_d * create_multi_UBspline_1d_d(Ugrid x_grid, BCtype_d xBC, int num_splines)
void set_multi_UBspline_2d_s(multi_UBspline_2d_s *spline, int num, float *data)
multi_UBspline_3d_c * create_multi_UBspline_3d_c(Ugrid x_grid, Ugrid y_grid, Ugrid z_grid, BCtype_c xBC, BCtype_c yBC, BCtype_c zBC, int num_splines)
multi_UBspline_2d_c * create_multi_UBspline_2d_c(Ugrid x_grid, Ugrid y_grid, BCtype_c xBC, BCtype_c yBC, int num_splines)
void set_multi_UBspline_1d_d_BC(multi_UBspline_1d_d *spline, int num, double *data, BCtype_d xBC)
multi_UBspline_1d_s * create_multi_UBspline_1d_s(Ugrid x_grid, BCtype_s xBC, int num_splines)
multi_UBspline_3d_s * create_multi_UBspline_3d_s(Ugrid x_grid, Ugrid y_grid, Ugrid z_grid, BCtype_s xBC, BCtype_s yBC, BCtype_s zBC, int num_splines)
void set_multi_UBspline_1d_z(multi_UBspline_1d_z *spline, int num, complex_double *data)
int posix_memalign(void **memptr, size_t alignment, size_t size)
void destroy_multi_UBspline(Bspline *spline)
void set_multi_UBspline_1d_s(multi_UBspline_1d_s *spline, int num, float *data)
void solve_antiperiodic_interp_1d_s(float bands[], float coefs[], int M, int cstride)
bc_code rCode
bc_code lCode
float lVal_r
float rVal_r
float rVal_i
float lVal_i
double lVal
bc_code rCode
double rVal
bc_code lCode
float lVal
float rVal
bc_code lCode
bc_code rCode
double lVal_r
double lVal_i
bc_code lCode
double rVal_i
bc_code rCode
double rVal_r
void *restrict coefs
double end
double delta
double start
double delta_inv
int num
complex_float *restrict coefs
complex_double *restrict coefs
complex_float *restrict coefs
complex_double *restrict coefs
complex_float *restrict coefs
complex_double *restrict coefs