Krita Source Code Documentation
Loading...
Searching...
No Matches
nubspline_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
21#include "nubspline_create.h"
22#include <math.h>
23#include <assert.h>
24#ifndef _XOPEN_SOURCE
25 #define _XOPEN_SOURCE 600
26#endif
27#ifndef __USE_XOPEN2K
28 #define __USE_XOPEN2K
29#endif
30#include <stdlib.h>
31#include <stdio.h>
32
34// Notes on conventions: //
35// Below, M (and Mx, My, Mz) represent the number of //
36// data points to be interpolated. With derivative //
37// boundary conditions, it is equal to the number of //
38// grid points. With periodic boundary conditions, //
39// it is one less than the number of grid points. //
40// N (and Nx, Ny, Nz) is the number of B-spline //
41// coefficients, which is #(grid points)+2 for all //
42// boundary conditions. //
44
45
51void
52solve_NUB_deriv_interp_1d_s (NUBasis* restrict basis,
53 float* restrict data, int datastride,
54 float* restrict p, int pstride,
55 float abcdInitial[4], float abcdFinal[4])
56{
57 int M = basis->grid->num_points;
58 int N = M+2;
59 // Banded matrix storage. The first three elements in the
60 // tinyvector store the tridiagonal coefficients. The last element
61 // stores the RHS data.
62#ifdef HAVE_C_VARARRAYS
63 float bands[4*N];
64#else
65 float *bands = new float[4*N];
66#endif
67
68
69 // Fill up bands
70 for (int i=0; i<4; i++) {
71 bands[i] = abcdInitial[i];
72 bands[4*(N-1)+i] = abcdFinal[i];
73 }
74 for (int i=0; i<M; i++) {
75 get_NUBasis_funcs_si (basis, i, &(bands[4*(i+1)]));
76 bands[4*(i+1)+3] = data[datastride*i];
77 }
78
79 // Now solve:
80 // First and last rows are different
81 bands[4*0+1] /= bands[4*0+0];
82 bands[4*0+2] /= bands[4*0+0];
83 bands[4*0+3] /= bands[4*0+0];
84 bands[4*0+0] = 1.0;
85 bands[4*1+1] -= bands[4*1+0]*bands[4*0+1];
86 bands[4*1+2] -= bands[4*1+0]*bands[4*0+2];
87 bands[4*1+3] -= bands[4*1+0]*bands[4*0+3];
88 bands[4*0+0] = 0.0;
89 bands[4*1+2] /= bands[4*1+1];
90 bands[4*1+3] /= bands[4*1+1];
91 bands[4*1+1] = 1.0;
92
93 // Now do rows 2 through M+1
94 for (int row=2; row < N-1; row++) {
95 bands[4*(row)+1] -= bands[4*(row)+0]*bands[4*(row-1)+2];
96 bands[4*(row)+3] -= bands[4*(row)+0]*bands[4*(row-1)+3];
97 bands[4*(row)+2] /= bands[4*(row)+1];
98 bands[4*(row)+3] /= bands[4*(row)+1];
99 bands[4*(row)+0] = 0.0;
100 bands[4*(row)+1] = 1.0;
101 }
102
103 // Do last row
104 bands[4*(M+1)+1] -= bands[4*(M+1)+0]*bands[4*(M-1)+2];
105 bands[4*(M+1)+3] -= bands[4*(M+1)+0]*bands[4*(M-1)+3];
106 bands[4*(M+1)+2] -= bands[4*(M+1)+1]*bands[4*(M)+2];
107 bands[4*(M+1)+3] -= bands[4*(M+1)+1]*bands[4*(M)+3];
108 bands[4*(M+1)+3] /= bands[4*(M+1)+2];
109 bands[4*(M+1)+2] = 1.0;
110
111 p[pstride*(M+1)] = bands[4*(M+1)+3];
112 // Now back substitute up
113 for (int row=M; row>0; row--)
114 p[pstride*(row)] = bands[4*(row)+3] - bands[4*(row)+2]*p[pstride*(row+1)];
115
116 // Finish with first row
117 p[0] = bands[4*(0)+3] - bands[4*(0)+1]*p[pstride*1] - bands[4*(0)+2]*p[pstride*2];
118
119#ifndef HAVE_C_VARARRAYS
120 delete[] bands;
121#endif
122}
123
124
125
126// The number of elements in data should be one less than the number
127// of grid points
128void
129solve_NUB_periodic_interp_1d_s (NUBasis* restrict basis,
130 float* restrict data, int datastride,
131 float* restrict p, int pstride)
132{
133 int M = basis->grid->num_points-1;
134
135 // Banded matrix storage. The first three elements in each row
136 // store the tridiagonal coefficients. The last element
137 // stores the RHS data.
138#ifdef HAVE_C_VARARRAYS
139 float bands[4*M], lastCol[M];
140#else
141 float *bands = new float[4*M];
142 float *lastCol = new float[ M];
143#endif
144
145 // Fill up bands
146 for (int i=0; i<M; i++) {
147 get_NUBasis_funcs_si (basis, i, &(bands[4*i]));
148 bands[4*(i)+3] = data[i*datastride];
149 }
150
151 // Now solve:
152 // First and last rows are different
153 bands[4*(0)+2] /= bands[4*(0)+1];
154 bands[4*(0)+0] /= bands[4*(0)+1];
155 bands[4*(0)+3] /= bands[4*(0)+1];
156 bands[4*(0)+1] = 1.0;
157 bands[4*(M-1)+1] -= bands[4*(M-1)+2]*bands[4*(0)+0];
158 bands[4*(M-1)+3] -= bands[4*(M-1)+2]*bands[4*(0)+3];
159 bands[4*(M-1)+2] = -bands[4*(M-1)+2]*bands[4*(0)+2];
160 lastCol[0] = bands[4*(0)+0];
161
162 for (int row=1; row < (M-1); row++) {
163 bands[4*(row)+1] -= bands[4*(row)+0] * bands[4*(row-1)+2];
164 bands[4*(row)+3] -= bands[4*(row)+0] * bands[4*(row-1)+3];
165 lastCol[row] = -bands[4*(row)+0] * lastCol[row-1];
166 bands[4*(row)+0] = 0.0;
167 bands[4*(row)+2] /= bands[4*(row)+1];
168 bands[4*(row)+3] /= bands[4*(row)+1];
169 lastCol[row] /= bands[4*(row)+1];
170 bands[4*(row)+1] = 1.0;
171 if (row < (M-2)) {
172 bands[4*(M-1)+3] -= bands[4*(M-1)+2]*bands[4*(row)+3];
173 bands[4*(M-1)+1] -= bands[4*(M-1)+2]*lastCol[row];
174 bands[4*(M-1)+2] = -bands[4*(M-1)+2]*bands[4*(row)+2];
175 }
176 }
177
178 // Now do last row
179 // The [2] element and [0] element are now on top of each other
180 bands[4*(M-1)+0] += bands[4*(M-1)+2];
181 bands[4*(M-1)+1] -= bands[4*(M-1)+0] * (bands[4*(M-2)+2]+lastCol[M-2]);
182 bands[4*(M-1)+3] -= bands[4*(M-1)+0] * bands[4*(M-2)+3];
183 bands[4*(M-1)+3] /= bands[4*(M-1)+1];
184 p[pstride*M] = bands[4*(M-1)+3];
185 for (int row=M-2; row>=0; row--)
186 p[pstride*(row+1)] = bands[4*(row)+3] -
187 bands[4*(row)+2]*p[pstride*(row+2)] - lastCol[row]*p[pstride*M];
188
189 p[pstride* 0 ] = p[pstride*M];
190 p[pstride*(M+1)] = p[pstride*1];
191 p[pstride*(M+2)] = p[pstride*2];
192#ifndef HAVE_C_VARARRAYS
193 delete[] bands;
194 delete[] lastCol;
195#endif
196}
197
198
199
200void
201find_NUBcoefs_1d_s (NUBasis* restrict basis, BCtype_s bc,
202 float *data, int dstride,
203 float *coefs, int cstride)
204{
205 if (bc.lCode == PERIODIC)
206 solve_NUB_periodic_interp_1d_s (basis, data, dstride, coefs, cstride);
207 else {
208 int M = basis->grid->num_points;
209 // Setup boundary conditions
210 float bfuncs[4] = {};
211 float dbfuncs[4] = {};
212 float abcd_left[4] = {};
213 float abcd_right[4] = {};
214 // Left boundary
215 if (bc.lCode == FLAT || bc.lCode == NATURAL)
216 bc.lVal = 0.0;
217 if (bc.lCode == FLAT || bc.lCode == DERIV1) {
218 get_NUBasis_dfuncs_si (basis, 0, bfuncs, abcd_left);
219 abcd_left[3] = bc.lVal;
220 }
221 if (bc.lCode == NATURAL || bc.lCode == DERIV2) {
222 get_NUBasis_d2funcs_si (basis, 0, bfuncs, dbfuncs, abcd_left);
223 abcd_left[3] = bc.lVal;
224 }
225
226 // Right boundary
227 if (bc.rCode == FLAT || bc.rCode == NATURAL)
228 bc.rVal = 0.0;
229 if (bc.rCode == FLAT || bc.rCode == DERIV1) {
230 get_NUBasis_dfuncs_si (basis, M-1, bfuncs, abcd_right);
231 abcd_right[3] = bc.rVal;
232 }
233 if (bc.rCode == NATURAL || bc.rCode == DERIV2) {
234 get_NUBasis_d2funcs_si (basis, M-1, bfuncs, dbfuncs, abcd_right);
235 abcd_right[3] = bc.rVal;
236 }
237 // Now, solve for coefficients
238 solve_NUB_deriv_interp_1d_s (basis, data, dstride, coefs, cstride,
239 abcd_left, abcd_right);
240 }
241}
242
243
244
245
246NUBspline_1d_s *
247create_NUBspline_1d_s (NUgrid* x_grid, BCtype_s xBC, float *data)
248{
249 // First, create the spline structure
250 NUBspline_1d_s* spline = static_cast<NUBspline_1d_s*>(malloc(sizeof(NUBspline_1d_s)));
251 if (spline == NULL)
252 return spline;
253 spline->sp_code = NU1D;
254 spline->t_code = SINGLE_REAL;
255
256 // Next, create the basis
257 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
258 // M is the number of data points (but is unused)
259// int M;
260// if (xBC.lCode == PERIODIC) M = x_grid->num_points - 1;
261// else M = x_grid->num_points;
262 int N = x_grid->num_points + 2;
263
264 // Allocate coefficients and solve
265 spline->coefs = (float*)malloc (sizeof(float)*N);
266 find_NUBcoefs_1d_s (spline->x_basis, xBC, data, 1, spline->coefs, 1);
267
268 return spline;
269}
270
271NUBspline_2d_s *
272create_NUBspline_2d_s (NUgrid* x_grid, NUgrid* y_grid,
273 BCtype_s xBC, BCtype_s yBC, float *data)
274{
275 // First, create the spline structure
276 NUBspline_2d_s* spline = static_cast<NUBspline_2d_s*>(malloc(sizeof(NUBspline_2d_s)));
277 if (spline == NULL)
278 return spline;
279 spline->sp_code = NU2D;
280 spline->t_code = SINGLE_REAL;
281
282 // Next, create the bases
283 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
284 spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
285 // set but unused
286 //int Mx,
287 int My, Nx, Ny;
288// if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
289// else Mx = x_grid->num_points;
290 if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
291 else My = y_grid->num_points;
292
293 Nx = x_grid->num_points + 2;
294 Ny = y_grid->num_points + 2;
295
296 spline->x_stride = Ny;
297#ifndef HAVE_SSE2
298 spline->coefs = (float*)malloc (sizeof(float)*Nx*Ny);
299#else
300 posix_memalign ((void**)&spline->coefs, 16, sizeof(float)*Nx*Ny);
301#endif
302
303 // First, solve in the X-direction
304 for (int iy=0; iy<My; iy++) {
305 int doffset = iy;
306 int coffset = iy;
307 find_NUBcoefs_1d_s (spline->x_basis, xBC, data+doffset, My,
308 spline->coefs+coffset, Ny);
309 }
310
311 // Now, solve in the Y-direction
312 for (int ix=0; ix<Nx; ix++) {
313 int doffset = ix*Ny;
314 int coffset = ix*Ny;
315 find_NUBcoefs_1d_s (spline->y_basis, yBC, spline->coefs+doffset, 1,
316 spline->coefs+coffset, 1);
317 }
318
319 return spline;
320}
321
322
323NUBspline_3d_s *
324create_NUBspline_3d_s (NUgrid* x_grid, NUgrid* y_grid, NUgrid* z_grid,
325 BCtype_s xBC, BCtype_s yBC, BCtype_s zBC, float *data)
326{
327 // First, create the spline structure
328 NUBspline_3d_s* spline = static_cast<NUBspline_3d_s*>(malloc(sizeof(NUBspline_3d_s)));
329 if (spline == NULL)
330 return spline;
331 spline->sp_code = NU3D;
332 spline->t_code = SINGLE_REAL;
333
334 // Next, create the bases
335 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
336 spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
337 spline->z_basis = create_NUBasis (z_grid, zBC.lCode==PERIODIC);
338 // set but unused
339 // int Mx,
340 int My, Mz, Nx, Ny, Nz;
341// if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
342// else Mx = x_grid->num_points;
343 if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
344 else My = y_grid->num_points;
345 if (zBC.lCode == PERIODIC) Mz = z_grid->num_points - 1;
346 else Mz = z_grid->num_points;
347
348 Nx = x_grid->num_points + 2;
349 Ny = y_grid->num_points + 2;
350 Nz = z_grid->num_points + 2;
351
352 // Allocate coefficients and solve
353 spline->x_stride = Ny*Nz;
354 spline->y_stride = Nz;
355#ifndef HAVE_SSE2
356 spline->coefs = (float*)malloc (sizeof(float)*Nx*Ny*Nz);
357#else
358 posix_memalign ((void**)&spline->coefs, 16, sizeof(float)*Nx*Ny*Nz);
359#endif
360
361 // First, solve in the X-direction
362 for (int iy=0; iy<My; iy++)
363 for (int iz=0; iz<Mz; iz++) {
364 int doffset = iy*Mz+iz;
365 int coffset = iy*Nz+iz;
366 find_NUBcoefs_1d_s (spline->x_basis, xBC, data+doffset, My*Mz,
367 spline->coefs+coffset, Ny*Nz);
368 }
369
370 // Now, solve in the Y-direction
371 for (int ix=0; ix<Nx; ix++)
372 for (int iz=0; iz<Nz; iz++) {
373 int doffset = ix*Ny*Nz + iz;
374 int coffset = ix*Ny*Nz + iz;
375 find_NUBcoefs_1d_s (spline->y_basis, yBC, spline->coefs+doffset, Nz,
376 spline->coefs+coffset, Nz);
377 }
378
379 // Now, solve in the Z-direction
380 for (int ix=0; ix<Nx; ix++)
381 for (int iy=0; iy<Ny; iy++) {
382 int doffset = (ix*Ny+iy)*Nz;
383 int coffset = (ix*Ny+iy)*Nz;
384 find_NUBcoefs_1d_s (spline->z_basis, zBC, spline->coefs+doffset, 1,
385 spline->coefs+coffset, 1);
386 }
387 return spline;
388}
389
395void
396solve_NUB_deriv_interp_1d_d (NUBasis* restrict basis,
397 double* restrict data, int datastride,
398 double* restrict p, int pstride,
399 double abcdInitial[4], double abcdFinal[4])
400{
401 int M = basis->grid->num_points;
402 int N = M+2;
403 // Banded matrix storage. The first three elements in the
404 // tinyvector store the tridiagonal coefficients. The last element
405 // stores the RHS data.
406#ifdef HAVE_C_VARARRAYS
407 double bands[4*N];
408#else
409 double *bands = new double[4*N];
410#endif
411
412 // Fill up bands
413 for (int i=0; i<4; i++) {
414 bands[i] = abcdInitial[i];
415 bands[4*(N-1)+i] = abcdFinal[i];
416 }
417 for (int i=0; i<M; i++) {
418 get_NUBasis_funcs_di (basis, i, &(bands[4*(i+1)]));
419 bands[4*(i+1)+3] = data[datastride*i];
420 }
421 /* for (int i=0; i<4*N; i++)
422 if (isnan(bands[i]))
423 fprintf(stderr, "NAN at bands[%d].\n", i); */
424
425 // Now solve:
426 // First and last rows are different
427 bands[4*0+1] /= bands[4*0+0];
428 bands[4*0+2] /= bands[4*0+0];
429 bands[4*0+3] /= bands[4*0+0];
430 bands[4*0+0] = 1.0;
431 bands[4*1+1] -= bands[4*1+0]*bands[4*0+1];
432 bands[4*1+2] -= bands[4*1+0]*bands[4*0+2];
433 bands[4*1+3] -= bands[4*1+0]*bands[4*0+3];
434 bands[4*0+0] = 0.0;
435 bands[4*1+2] /= bands[4*1+1];
436 bands[4*1+3] /= bands[4*1+1];
437 bands[4*1+1] = 1.0;
438
439 // Now do rows 2 through M+1
440 for (int row=2; row < N-1; row++) {
441 bands[4*(row)+1] -= bands[4*(row)+0]*bands[4*(row-1)+2];
442 bands[4*(row)+3] -= bands[4*(row)+0]*bands[4*(row-1)+3];
443 bands[4*(row)+2] /= bands[4*(row)+1];
444 bands[4*(row)+3] /= bands[4*(row)+1];
445 bands[4*(row)+0] = 0.0;
446 bands[4*(row)+1] = 1.0;
447 }
448
449 // Do last row
450 bands[4*(M+1)+1] -= bands[4*(M+1)+0]*bands[4*(M-1)+2];
451 bands[4*(M+1)+3] -= bands[4*(M+1)+0]*bands[4*(M-1)+3];
452 bands[4*(M+1)+2] -= bands[4*(M+1)+1]*bands[4*(M)+2];
453 bands[4*(M+1)+3] -= bands[4*(M+1)+1]*bands[4*(M)+3];
454 bands[4*(M+1)+3] /= bands[4*(M+1)+2];
455 bands[4*(M+1)+2] = 1.0;
456
457 p[pstride*(M+1)] = bands[4*(M+1)+3];
458 // Now back substitute up
459 for (int row=M; row>0; row--)
460 p[pstride*(row)] = bands[4*(row)+3] - bands[4*(row)+2]*p[pstride*(row+1)];
461
462 // Finish with first row
463 p[0] = bands[4*(0)+3] - bands[4*(0)+1]*p[pstride*1] - bands[4*(0)+2]*p[pstride*2];
464#ifndef HAVE_C_VARARRAYS
465 delete[] bands;
466#endif
467}
468
469
470void
471solve_NUB_periodic_interp_1d_d (NUBasis* restrict basis,
472 double* restrict data, int datastride,
473 double* restrict p, int pstride)
474{
475 int M = basis->grid->num_points-1;
476
477 // Banded matrix storage. The first three elements in the
478 // tinyvector store the tridiagonal coefficients. The last element
479 // stores the RHS data.
480#ifdef HAVE_C_VARARRAYS
481 double bands[4*M], lastCol[M];
482#else
483 double *bands = new double[4*M];
484 double *lastCol = new double[ M];
485#endif
486
487 // Fill up bands
488 for (int i=0; i<M; i++) {
489 get_NUBasis_funcs_di (basis, i, &(bands[4*i]));
490 bands[4*(i)+3] = data[i*datastride];
491 }
492
493 // Now solve:
494 // First and last rows are different
495 bands[4*(0)+2] /= bands[4*(0)+1];
496 bands[4*(0)+0] /= bands[4*(0)+1];
497 bands[4*(0)+3] /= bands[4*(0)+1];
498 bands[4*(0)+1] = 1.0;
499 bands[4*(M-1)+1] -= bands[4*(M-1)+2]*bands[4*(0)+0];
500 bands[4*(M-1)+3] -= bands[4*(M-1)+2]*bands[4*(0)+3];
501 bands[4*(M-1)+2] = -bands[4*(M-1)+2]*bands[4*(0)+2];
502 lastCol[0] = bands[4*(0)+0];
503
504 for (int row=1; row < (M-1); row++) {
505 bands[4*(row)+1] -= bands[4*(row)+0] * bands[4*(row-1)+2];
506 bands[4*(row)+3] -= bands[4*(row)+0] * bands[4*(row-1)+3];
507 lastCol[row] = -bands[4*(row)+0] * lastCol[row-1];
508 bands[4*(row)+0] = 0.0;
509 bands[4*(row)+2] /= bands[4*(row)+1];
510 bands[4*(row)+3] /= bands[4*(row)+1];
511 lastCol[row] /= bands[4*(row)+1];
512 bands[4*(row)+1] = 1.0;
513 if (row < (M-2)) {
514 bands[4*(M-1)+3] -= bands[4*(M-1)+2]*bands[4*(row)+3];
515 bands[4*(M-1)+1] -= bands[4*(M-1)+2]*lastCol[row];
516 bands[4*(M-1)+2] = -bands[4*(M-1)+2]*bands[4*(row)+2];
517 }
518 }
519
520 // Now do last row
521 // The [2] element and [0] element are now on top of each other
522 bands[4*(M-1)+0] += bands[4*(M-1)+2];
523 bands[4*(M-1)+1] -= bands[4*(M-1)+0] * (bands[4*(M-2)+2]+lastCol[M-2]);
524 bands[4*(M-1)+3] -= bands[4*(M-1)+0] * bands[4*(M-2)+3];
525 bands[4*(M-1)+3] /= bands[4*(M-1)+1];
526 p[pstride*M] = bands[4*(M-1)+3];
527 for (int row=M-2; row>=0; row--)
528 p[pstride*(row+1)] = bands[4*(row)+3] -
529 bands[4*(row)+2]*p[pstride*(row+2)] - lastCol[row]*p[pstride*M];
530
531 p[pstride* 0 ] = p[pstride*M];
532 p[pstride*(M+1)] = p[pstride*1];
533 p[pstride*(M+2)] = p[pstride*2];
534#ifndef HAVE_C_VARARRAYS
535 delete[] bands;
536 delete[] lastCol;
537#endif
538}
539
540
541
542void
543find_NUBcoefs_1d_d (NUBasis* restrict basis, BCtype_d bc,
544 double *data, int dstride,
545 double *coefs, int cstride)
546{
547 if (bc.lCode == PERIODIC)
548 solve_NUB_periodic_interp_1d_d (basis, data, dstride, coefs, cstride);
549 else {
550 int M = basis->grid->num_points;
551 // Setup boundary conditions
552 double bfuncs[4] {};
553 double dbfuncs[4] {};
554 double abcd_left[4] {};
555 double abcd_right[4] {};
556
557 // Left boundary
558 if (bc.lCode == FLAT || bc.lCode == NATURAL)
559 bc.lVal = 0.0;
560 if (bc.lCode == FLAT || bc.lCode == DERIV1) {
561 get_NUBasis_dfuncs_di (basis, 0, bfuncs, abcd_left);
562 abcd_left[3] = bc.lVal;
563 }
564 if (bc.lCode == NATURAL || bc.lCode == DERIV2) {
565 get_NUBasis_d2funcs_di (basis, 0, bfuncs, dbfuncs, abcd_left);
566 abcd_left[3] = bc.lVal;
567 }
568
569 // Right boundary
570 if (bc.rCode == FLAT || bc.rCode == NATURAL)
571 bc.rVal = 0.0;
572 if (bc.rCode == FLAT || bc.rCode == DERIV1) {
573 get_NUBasis_dfuncs_di (basis, M-1, bfuncs, abcd_right);
574 abcd_right[3] = bc.rVal;
575 }
576 if (bc.rCode == NATURAL || bc.rCode == DERIV2) {
577 get_NUBasis_d2funcs_di (basis, M-1, bfuncs, dbfuncs, abcd_right);
578 abcd_right[3] = bc.rVal;
579 }
580
581 // Now, solve for coefficients
582 solve_NUB_deriv_interp_1d_d (basis, data, dstride, coefs, cstride,
583 abcd_left, abcd_right);
584 }
585}
586
587
588
589
590NUBspline_1d_d *
591create_NUBspline_1d_d (NUgrid* x_grid, BCtype_d xBC, double *data)
592{
593 // First, create the spline structure
594 NUBspline_1d_d* spline = static_cast<NUBspline_1d_d*>(malloc(sizeof(NUBspline_1d_d)));
595 if (spline == NULL)
596 return spline;
597 spline->sp_code = NU1D;
598 spline->t_code = DOUBLE_REAL;
599
600 // Next, create the basis
601 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
602 // M is the number of data points (set but unused)
603// int M;
604// if (xBC.lCode == PERIODIC) M = x_grid->num_points - 1;
605// else M = x_grid->num_points;
606 int N = x_grid->num_points + 2;
607
608 // Allocate coefficients and solve
609 spline->coefs = (double*)malloc (sizeof(double)*N);
610 find_NUBcoefs_1d_d (spline->x_basis, xBC, data, 1, spline->coefs, 1);
611
612 return spline;
613}
614
615NUBspline_2d_d *
616create_NUBspline_2d_d (NUgrid* x_grid, NUgrid* y_grid,
617 BCtype_d xBC, BCtype_d yBC, double *data)
618{
619 // First, create the spline structure
620 NUBspline_2d_d* spline = static_cast<NUBspline_2d_d*>(malloc(sizeof(NUBspline_2d_d)));
621 if (spline == NULL)
622 return spline;
623 spline->sp_code = NU2D;
624 spline->t_code = DOUBLE_REAL;
625
626 // Next, create the bases
627 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
628 spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
629
630 // int Mx, (set but unused)
631 int My, Nx, Ny;
632// if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
633// else Mx = x_grid->num_points;
634 if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
635 else My = y_grid->num_points;
636
637 Nx = x_grid->num_points + 2;
638 Ny = y_grid->num_points + 2;
639
640 spline->x_stride = Ny;
641#ifndef HAVE_SSE2
642 spline->coefs = (double*)malloc (sizeof(double)*Nx*Ny);
643#else
644 posix_memalign ((void**)&spline->coefs, 16, sizeof(double)*Nx*Ny);
645#endif
646
647 // First, solve in the X-direction
648 for (int iy=0; iy<My; iy++) {
649 int doffset = iy;
650 int coffset = iy;
651 find_NUBcoefs_1d_d (spline->x_basis, xBC, data+doffset, My,
652 spline->coefs+coffset, Ny);
653 }
654
655 // Now, solve in the Y-direction
656 for (int ix=0; ix<Nx; ix++) {
657 int doffset = ix*Ny;
658 int coffset = ix*Ny;
659 find_NUBcoefs_1d_d (spline->y_basis, yBC, spline->coefs+doffset, 1,
660 spline->coefs+coffset, 1);
661 }
662
663 return spline;
664}
665
666
667NUBspline_3d_d *
668create_NUBspline_3d_d (NUgrid* x_grid, NUgrid* y_grid, NUgrid* z_grid,
669 BCtype_d xBC, BCtype_d yBC, BCtype_d zBC, double *data)
670{
671 // First, create the spline structure
672 NUBspline_3d_d* spline = static_cast<NUBspline_3d_d*>(malloc(sizeof(NUBspline_3d_d)));
673 if (spline == NULL)
674 return spline;
675 spline->sp_code = NU3D;
676 spline->t_code = DOUBLE_REAL;
677
678 // Next, create the bases
679 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
680 spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
681 spline->z_basis = create_NUBasis (z_grid, zBC.lCode==PERIODIC);
682
683 // set but unused
684 //int Mx,
685 int My, Mz, Nx, Ny, Nz;
686// if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
687// else Mx = x_grid->num_points;
688 if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
689 else My = y_grid->num_points;
690 if (zBC.lCode == PERIODIC) Mz = z_grid->num_points - 1;
691 else Mz = z_grid->num_points;
692
693 Nx = x_grid->num_points + 2;
694 Ny = y_grid->num_points + 2;
695 Nz = z_grid->num_points + 2;
696
697 spline->x_stride = Ny*Nz;
698 spline->y_stride = Nz;
699#ifndef HAVE_SSE2
700 spline->coefs = (double*)malloc (sizeof(double)*Nx*Ny*Nz);
701#else
702 posix_memalign ((void**)&spline->coefs, 16, sizeof(double)*Nx*Ny*Nz);
703#endif
704
705 // First, solve in the X-direction
706 for (int iy=0; iy<My; iy++)
707 for (int iz=0; iz<Mz; iz++) {
708 int doffset = iy*Mz+iz;
709 int coffset = iy*Nz+iz;
710 find_NUBcoefs_1d_d (spline->x_basis, xBC, data+doffset, My*Mz,
711 spline->coefs+coffset, Ny*Nz);
712 }
713
714 // Now, solve in the Y-direction
715 for (int ix=0; ix<Nx; ix++)
716 for (int iz=0; iz<Nz; iz++) {
717 int doffset = ix*Ny*Nz + iz;
718 int coffset = ix*Ny*Nz + iz;
719 find_NUBcoefs_1d_d (spline->y_basis, yBC, spline->coefs+doffset, Nz,
720 spline->coefs+coffset, Nz);
721 }
722
723 // Now, solve in the Z-direction
724 for (int ix=0; ix<Nx; ix++)
725 for (int iy=0; iy<Ny; iy++) {
726 int doffset = (ix*Ny+iy)*Nz;
727 int coffset = (ix*Ny+iy)*Nz;
728 find_NUBcoefs_1d_d (spline->z_basis, zBC, spline->coefs+doffset, 1,
729 spline->coefs+coffset, 1);
730 }
731 return spline;
732}
733
734
740
741void
742find_NUBcoefs_1d_c (NUBasis* restrict basis, BCtype_c bc,
743 complex_float *data, int dstride,
744 complex_float *coefs, int cstride)
745{
746 BCtype_s bc_r, bc_i;
747 bc_r.lCode = bc.lCode; bc_i.lCode = bc.lCode;
748 bc_r.rCode = bc.rCode; bc_i.rCode = bc.rCode;
749 bc_r.lVal = bc.lVal_r; bc_r.rVal = bc.rVal_r;
750 bc_i.lVal = bc.lVal_i; bc_i.rVal = bc.rVal_i;
751
752 float *data_r = ((float*)data );
753 float *data_i = ((float*)data )+1;
754 float *coefs_r = ((float*)coefs);
755 float *coefs_i = ((float*)coefs)+1;
756
757 find_NUBcoefs_1d_s (basis, bc_r, data_r, 2*dstride, coefs_r, 2*cstride);
758 find_NUBcoefs_1d_s (basis, bc_i, data_i, 2*dstride, coefs_i, 2*cstride);
759}
760
761
762NUBspline_1d_c *
763create_NUBspline_1d_c (NUgrid* x_grid, BCtype_c xBC, complex_float *data)
764{
765 // First, create the spline structure
766 NUBspline_1d_c* spline = static_cast<NUBspline_1d_c*>(malloc(sizeof(NUBspline_1d_c)));
767 if (spline == NULL)
768 return spline;
769 spline->sp_code = NU1D;
770 spline->t_code = SINGLE_COMPLEX;
771
772 // Next, create the basis
773 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
774 // M is the number of data points
775// int M;
776// if (xBC.lCode == PERIODIC) M = x_grid->num_points - 1;
777// else M = x_grid->num_points;
778 int N = x_grid->num_points + 2;
779
780 // Allocate coefficients and solve
781 spline->coefs = (complex_float*)malloc (sizeof(complex_float)*N);
782 find_NUBcoefs_1d_c (spline->x_basis, xBC, data, 1, spline->coefs, 1);
783
784 return spline;
785}
786
787NUBspline_2d_c *
788create_NUBspline_2d_c (NUgrid* x_grid, NUgrid* y_grid,
789 BCtype_c xBC, BCtype_c yBC, complex_float *data)
790{
791 // First, create the spline structure
792 NUBspline_2d_c* spline = static_cast<NUBspline_2d_c*>(malloc(sizeof(NUBspline_2d_c)));
793 if (spline == NULL)
794 return spline;
795 spline->sp_code = NU2D;
796 spline->t_code = SINGLE_COMPLEX;
797
798 // Next, create the bases
799 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
800 spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
801// int Mx,
802 int My, Nx, Ny;
803// if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
804// else Mx = x_grid->num_points;
805 if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
806 else My = y_grid->num_points;
807
808 Nx = x_grid->num_points + 2;
809 Ny = y_grid->num_points + 2;
810
811 spline->x_stride = Ny;
812#ifndef HAVE_SSE2
813 spline->coefs = (complex_float*)malloc (sizeof(complex_float)*Nx*Ny);
814#else
815 posix_memalign ((void**)&spline->coefs, 16, sizeof(complex_float)*Nx*Ny);
816#endif
817
818 // First, solve in the X-direction
819 for (int iy=0; iy<My; iy++) {
820 int doffset = iy;
821 int coffset = iy;
822 find_NUBcoefs_1d_c (spline->x_basis, xBC, data+doffset, My,
823 spline->coefs+coffset, Ny);
824 }
825
826 // Now, solve in the Y-direction
827 for (int ix=0; ix<Nx; ix++) {
828 int doffset = ix*Ny;
829 int coffset = ix*Ny;
830 find_NUBcoefs_1d_c (spline->y_basis, yBC, spline->coefs+doffset, 1,
831 spline->coefs+coffset, 1);
832 }
833
834 return spline;
835}
836
837
838NUBspline_3d_c *
839create_NUBspline_3d_c (NUgrid* x_grid, NUgrid* y_grid, NUgrid* z_grid,
840 BCtype_c xBC, BCtype_c yBC, BCtype_c zBC, complex_float *data)
841{
842 // First, create the spline structure
843 NUBspline_3d_c* spline = static_cast<NUBspline_3d_c*>(malloc(sizeof(NUBspline_3d_c)));
844 if (spline == NULL)
845 return spline;
846 spline->sp_code = NU3D;
847 spline->t_code = SINGLE_COMPLEX;
848
849 // Next, create the bases
850 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
851 spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
852 spline->z_basis = create_NUBasis (z_grid, zBC.lCode==PERIODIC);
853// int Mx,
854 int My, Mz, Nx, Ny, Nz;
855// if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
856// else Mx = x_grid->num_points;
857 if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
858 else My = y_grid->num_points;
859 if (zBC.lCode == PERIODIC) Mz = z_grid->num_points - 1;
860 else Mz = z_grid->num_points;
861
862 Nx = x_grid->num_points + 2;
863 Ny = y_grid->num_points + 2;
864 Nz = z_grid->num_points + 2;
865
866 // Allocate coefficients and solve
867 spline->x_stride = Ny*Nz;
868 spline->y_stride = Nz;
869#ifndef HAVE_SSE2
870 spline->coefs = (complex_float*)malloc (sizeof(complex_float)*Nx*Ny*Nz);
871#else
872 posix_memalign ((void**)&spline->coefs, 16, sizeof(complex_float)*Nx*Ny*Nz);
873#endif
874
875 // First, solve in the X-direction
876 for (int iy=0; iy<My; iy++)
877 for (int iz=0; iz<Mz; iz++) {
878 int doffset = iy*Mz+iz;
879 int coffset = iy*Nz+iz;
880 find_NUBcoefs_1d_c (spline->x_basis, xBC, data+doffset, My*Mz,
881 spline->coefs+coffset, Ny*Nz);
882 }
883
884 // Now, solve in the Y-direction
885 for (int ix=0; ix<Nx; ix++)
886 for (int iz=0; iz<Nz; iz++) {
887 int doffset = ix*Ny*Nz + iz;
888 int coffset = ix*Ny*Nz + iz;
889 find_NUBcoefs_1d_c (spline->y_basis, yBC, spline->coefs+doffset, Nz,
890 spline->coefs+coffset, Nz);
891 }
892
893 // Now, solve in the Z-direction
894 for (int ix=0; ix<Nx; ix++)
895 for (int iy=0; iy<Ny; iy++) {
896 int doffset = (ix*Ny+iy)*Nz;
897 int coffset = (ix*Ny+iy)*Nz;
898 find_NUBcoefs_1d_c (spline->z_basis, zBC, spline->coefs+doffset, 1,
899 spline->coefs+coffset, 1);
900 }
901 return spline;
902}
903
909
910void
911find_NUBcoefs_1d_z (NUBasis* restrict basis, BCtype_z bc,
912 complex_double *data, int dstride,
913 complex_double *coefs, int cstride)
914{
915 BCtype_d bc_r, bc_i;
916 bc_r.lCode = bc.lCode; bc_i.lCode = bc.lCode;
917 bc_r.rCode = bc.rCode; bc_i.rCode = bc.rCode;
918 bc_r.lVal = bc.lVal_r; bc_r.rVal = bc.rVal_r;
919 bc_i.lVal = bc.lVal_i; bc_i.rVal = bc.rVal_i;
920
921 double *data_r = ((double*)data );
922 double *data_i = ((double*)data )+1;
923 double *coefs_r = ((double*)coefs);
924 double *coefs_i = ((double*)coefs)+1;
925
926 find_NUBcoefs_1d_d (basis, bc_r, data_r, 2*dstride, coefs_r, 2*cstride);
927 find_NUBcoefs_1d_d (basis, bc_i, data_i, 2*dstride, coefs_i, 2*cstride);
928}
929
930
931NUBspline_1d_z *
932create_NUBspline_1d_z (NUgrid* x_grid, BCtype_z xBC, complex_double *data)
933{
934 // First, create the spline structure
935 NUBspline_1d_z* spline = static_cast<NUBspline_1d_z*>(malloc(sizeof(NUBspline_1d_z)));
936 if (spline == NULL)
937 return spline;
938 spline->sp_code = NU1D;
939 spline->t_code = DOUBLE_COMPLEX;
940
941 // Next, create the basis
942 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
943 // M is the number of data points
944// int M;
945// if (xBC.lCode == PERIODIC) M = x_grid->num_points - 1;
946// else M = x_grid->num_points;
947 int N = x_grid->num_points + 2;
948
949 // Allocate coefficients and solve
950 spline->coefs = (complex_double*)malloc (sizeof(complex_double)*N);
951 find_NUBcoefs_1d_z (spline->x_basis, xBC, data, 1, spline->coefs, 1);
952
953 return spline;
954}
955
956NUBspline_2d_z *
957create_NUBspline_2d_z (NUgrid* x_grid, NUgrid* y_grid,
958 BCtype_z xBC, BCtype_z yBC, complex_double *data)
959{
960 // First, create the spline structure
961 NUBspline_2d_z* spline = static_cast<NUBspline_2d_z*>(malloc(sizeof(NUBspline_2d_z)));
962 if (spline == NULL)
963 return spline;
964 spline->sp_code = NU2D;
965 spline->t_code = DOUBLE_COMPLEX;
966
967 // Next, create the bases
968 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
969 spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
970// int Mx,
971 int My, Nx, Ny;
972// if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
973// else Mx = x_grid->num_points;
974 if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
975 else My = y_grid->num_points;
976
977 Nx = x_grid->num_points + 2;
978 Ny = y_grid->num_points + 2;
979
980 spline->x_stride = Ny;
981#ifndef HAVE_SSE2
982 spline->coefs = (complex_double*)malloc (sizeof(complex_double)*Nx*Ny);
983#else
984 posix_memalign ((void**)&spline->coefs, 16, sizeof(complex_double)*Nx*Ny);
985#endif
986
987 // First, solve in the X-direction
988 for (int iy=0; iy<My; iy++) {
989 int doffset = iy;
990 int coffset = iy;
991 find_NUBcoefs_1d_z (spline->x_basis, xBC, data+doffset, My,
992 spline->coefs+coffset, Ny);
993 }
994
995 // Now, solve in the Y-direction
996 for (int ix=0; ix<Nx; ix++) {
997 int doffset = ix*Ny;
998 int coffset = ix*Ny;
999 find_NUBcoefs_1d_z (spline->y_basis, yBC, spline->coefs+doffset, 1,
1000 spline->coefs+coffset, 1);
1001 }
1002
1003 return spline;
1004}
1005
1006
1007NUBspline_3d_z *
1008create_NUBspline_3d_z (NUgrid* x_grid, NUgrid* y_grid, NUgrid* z_grid,
1009 BCtype_z xBC, BCtype_z yBC, BCtype_z zBC, complex_double *data)
1010{
1011 // First, create the spline structure
1012 NUBspline_3d_z* spline = static_cast<NUBspline_3d_z*>(malloc(sizeof(NUBspline_3d_z)));
1013 if (spline == NULL)
1014 return spline;
1015 spline->sp_code = NU3D;
1016 spline->t_code = DOUBLE_COMPLEX;
1017 spline->x_grid = x_grid;
1018 spline->y_grid = y_grid;
1019 spline->z_grid = z_grid;
1020
1021 // Next, create the bases
1022 spline->x_basis = create_NUBasis (x_grid, xBC.lCode==PERIODIC);
1023 spline->y_basis = create_NUBasis (y_grid, yBC.lCode==PERIODIC);
1024 spline->z_basis = create_NUBasis (z_grid, zBC.lCode==PERIODIC);
1025// int Mx,
1026 int My, Mz, Nx, Ny, Nz;
1027// if (xBC.lCode == PERIODIC) Mx = x_grid->num_points - 1;
1028// else Mx = x_grid->num_points;
1029 if (yBC.lCode == PERIODIC) My = y_grid->num_points - 1;
1030 else My = y_grid->num_points;
1031 if (zBC.lCode == PERIODIC) Mz = z_grid->num_points - 1;
1032 else Mz = z_grid->num_points;
1033
1034 Nx = x_grid->num_points + 2;
1035 Ny = y_grid->num_points + 2;
1036 Nz = z_grid->num_points + 2;
1037
1038 // Allocate coefficients and solve
1039 spline->x_stride = Ny*Nz;
1040 spline->y_stride = Nz;
1041#ifndef HAVE_SSE2
1042 spline->coefs = (complex_double*)malloc (sizeof(complex_double)*Nx*Ny*Nz);
1043#else
1044 posix_memalign ((void**)&spline->coefs, 16, sizeof(complex_double)*Nx*Ny*Nz);
1045#endif
1046
1047 // First, solve in the X-direction
1048 for (int iy=0; iy<My; iy++)
1049 for (int iz=0; iz<Mz; iz++) {
1050 int doffset = iy*Mz+iz;
1051 int coffset = iy*Nz+iz;
1052 find_NUBcoefs_1d_z (spline->x_basis, xBC, data+doffset, My*Mz,
1053 spline->coefs+coffset, Ny*Nz);
1054 /* for (int ix=0; ix<Nx; ix++) {
1055 complex_double z = spline->coefs[coffset+ix*spline->x_stride];
1056 if (isnan(creal(z)))
1057 fprintf (stderr, "NAN encountered in create_NUBspline_3d_z at real part of (%d,%d,%d)\n",
1058 ix,iy,iz);
1059 if (isnan(cimag(z)))
1060 fprintf (stderr, "NAN encountered in create_NUBspline_3d_z at imag part of (%d,%d,%d)\n",
1061 ix,iy,iz);
1062 } */
1063 }
1064
1065 // Now, solve in the Y-direction
1066 for (int ix=0; ix<Nx; ix++)
1067 for (int iz=0; iz<Nz; iz++) {
1068 int doffset = ix*Ny*Nz + iz;
1069 int coffset = ix*Ny*Nz + iz;
1070 find_NUBcoefs_1d_z (spline->y_basis, yBC, spline->coefs+doffset, Nz,
1071 spline->coefs+coffset, Nz);
1072 }
1073
1074 // Now, solve in the Z-direction
1075 for (int ix=0; ix<Nx; ix++)
1076 for (int iy=0; iy<Ny; iy++) {
1077 int doffset = (ix*Ny+iy)*Nz;
1078 int coffset = (ix*Ny+iy)*Nz;
1079 find_NUBcoefs_1d_z (spline->z_basis, zBC, spline->coefs+doffset, 1,
1080 spline->coefs+coffset, 1);
1081 }
1082 return spline;
1083}
1084
1085
1086void
1087destroy_NUBspline(Bspline *spline)
1088{
1089 free(spline->coefs);
1090 switch (spline->sp_code) {
1091 case NU1D:
1092 destroy_NUBasis (((NUBspline_1d*)spline)->x_basis);
1093 break;
1094 case NU2D:
1095 destroy_NUBasis (((NUBspline_2d*)spline)->x_basis);
1096 destroy_NUBasis (((NUBspline_2d*)spline)->y_basis);
1097 break;
1098
1099 case NU3D:
1100 destroy_NUBasis (((NUBspline_3d*)spline)->x_basis);
1101 destroy_NUBasis (((NUBspline_3d*)spline)->y_basis);
1102 destroy_NUBasis (((NUBspline_3d*)spline)->z_basis);
1103 break;
1104 case U1D:
1105 case U2D:
1106 case MULTI_U1D:
1107 case MULTI_U2D:
1108 case MULTI_U3D:
1109 case MULTI_NU1D:
1110 case MULTI_NU2D:
1111 case MULTI_NU3D:
1112 default:
1113 ;
1114 }
1115 free(spline);
1116}
1117
p
const Params2D p
Definition KisBezierUtils.cpp:703
complex_float
complex float complex_float
Definition bspline_base.h:33
DOUBLE_REAL
@ DOUBLE_REAL
Definition bspline_base.h:55
SINGLE_REAL
@ SINGLE_REAL
Definition bspline_base.h:55
DOUBLE_COMPLEX
@ DOUBLE_COMPLEX
Definition bspline_base.h:55
SINGLE_COMPLEX
@ SINGLE_COMPLEX
Definition bspline_base.h:55
complex_double
complex double complex_double
Definition bspline_base.h:34
FLAT
@ FLAT
Definition bspline_base.h:50
NATURAL
@ NATURAL
Definition bspline_base.h:50
DERIV1
@ DERIV1
Definition bspline_base.h:50
DERIV2
@ DERIV2
Definition bspline_base.h:50
PERIODIC
@ PERIODIC
Definition bspline_base.h:50
MULTI_NU1D
@ MULTI_NU1D
Definition bspline_base.h:54
NU1D
@ NU1D
Definition bspline_base.h:52
MULTI_NU2D
@ MULTI_NU2D
Definition bspline_base.h:54
U1D
@ U1D
Definition bspline_base.h:51
NU2D
@ NU2D
Definition bspline_base.h:52
MULTI_U1D
@ MULTI_U1D
Definition bspline_base.h:53
MULTI_U2D
@ MULTI_U2D
Definition bspline_base.h:53
MULTI_U3D
@ MULTI_U3D
Definition bspline_base.h:53
MULTI_NU3D
@ MULTI_NU3D
Definition bspline_base.h:54
U2D
@ U2D
Definition bspline_base.h:51
NU3D
@ NU3D
Definition bspline_base.h:52
posix_memalign
int posix_memalign(void **memptr, size_t alignment, size_t size)
restrict
#define restrict
Definition local_definitions.h:4
get_NUBasis_funcs_si
void get_NUBasis_funcs_si(NUBasis *restrict basis, int i, float bfuncs[4])
Definition nubasis.cpp:97
get_NUBasis_dfuncs_si
void get_NUBasis_dfuncs_si(NUBasis *restrict basis, int i, float bfuncs[4], float dbfuncs[4])
Definition nubasis.cpp:153
destroy_NUBasis
void destroy_NUBasis(NUBasis *basis)
Definition nubasis.cpp:62
create_NUBasis
NUBasis * create_NUBasis(NUgrid *grid, bool periodic)
Definition nubasis.cpp:27
get_NUBasis_d2funcs_di
void get_NUBasis_d2funcs_di(NUBasis *restrict basis, int i, double bfuncs[4], double dbfuncs[4], double d2bfuncs[4])
Definition nubasis.cpp:413
get_NUBasis_dfuncs_di
void get_NUBasis_dfuncs_di(NUBasis *restrict basis, int i, double bfuncs[4], double dbfuncs[4])
Definition nubasis.cpp:344
get_NUBasis_funcs_di
void get_NUBasis_funcs_di(NUBasis *restrict basis, int i, double bfuncs[4])
Definition nubasis.cpp:288
get_NUBasis_d2funcs_si
void get_NUBasis_d2funcs_si(NUBasis *restrict basis, int i, float bfuncs[4], float dbfuncs[4], float d2bfuncs[4])
Definition nubasis.cpp:222
create_NUBspline_1d_c
NUBspline_1d_c * create_NUBspline_1d_c(NUgrid *x_grid, BCtype_c xBC, complex_float *data)
Definition nubspline_create.cpp:763
create_NUBspline_3d_s
NUBspline_3d_s * create_NUBspline_3d_s(NUgrid *x_grid, NUgrid *y_grid, NUgrid *z_grid, BCtype_s xBC, BCtype_s yBC, BCtype_s zBC, float *data)
Definition nubspline_create.cpp:324
find_NUBcoefs_1d_c
void find_NUBcoefs_1d_c(NUBasis *restrict basis, BCtype_c bc, complex_float *data, int dstride, complex_float *coefs, int cstride)
Definition nubspline_create.cpp:742
create_NUBspline_2d_z
NUBspline_2d_z * create_NUBspline_2d_z(NUgrid *x_grid, NUgrid *y_grid, BCtype_z xBC, BCtype_z yBC, complex_double *data)
Definition nubspline_create.cpp:957
find_NUBcoefs_1d_s
void find_NUBcoefs_1d_s(NUBasis *restrict basis, BCtype_s bc, float *data, int dstride, float *coefs, int cstride)
Definition nubspline_create.cpp:201
find_NUBcoefs_1d_z
void find_NUBcoefs_1d_z(NUBasis *restrict basis, BCtype_z bc, complex_double *data, int dstride, complex_double *coefs, int cstride)
Definition nubspline_create.cpp:911
destroy_NUBspline
void destroy_NUBspline(Bspline *spline)
Definition nubspline_create.cpp:1087
solve_NUB_periodic_interp_1d_d
void solve_NUB_periodic_interp_1d_d(NUBasis *restrict basis, double *restrict data, int datastride, double *restrict p, int pstride)
Definition nubspline_create.cpp:471
create_NUBspline_1d_d
NUBspline_1d_d * create_NUBspline_1d_d(NUgrid *x_grid, BCtype_d xBC, double *data)
Definition nubspline_create.cpp:591
create_NUBspline_1d_s
NUBspline_1d_s * create_NUBspline_1d_s(NUgrid *x_grid, BCtype_s xBC, float *data)
Definition nubspline_create.cpp:247
solve_NUB_deriv_interp_1d_s
void solve_NUB_deriv_interp_1d_s(NUBasis *restrict basis, float *restrict data, int datastride, float *restrict p, int pstride, float abcdInitial[4], float abcdFinal[4])
Definition nubspline_create.cpp:52
create_NUBspline_2d_c
NUBspline_2d_c * create_NUBspline_2d_c(NUgrid *x_grid, NUgrid *y_grid, BCtype_c xBC, BCtype_c yBC, complex_float *data)
Definition nubspline_create.cpp:788
solve_NUB_periodic_interp_1d_s
void solve_NUB_periodic_interp_1d_s(NUBasis *restrict basis, float *restrict data, int datastride, float *restrict p, int pstride)
Definition nubspline_create.cpp:129
create_NUBspline_3d_z
NUBspline_3d_z * create_NUBspline_3d_z(NUgrid *x_grid, NUgrid *y_grid, NUgrid *z_grid, BCtype_z xBC, BCtype_z yBC, BCtype_z zBC, complex_double *data)
Definition nubspline_create.cpp:1008
create_NUBspline_2d_s
NUBspline_2d_s * create_NUBspline_2d_s(NUgrid *x_grid, NUgrid *y_grid, BCtype_s xBC, BCtype_s yBC, float *data)
Definition nubspline_create.cpp:272
solve_NUB_deriv_interp_1d_d
void solve_NUB_deriv_interp_1d_d(NUBasis *restrict basis, double *restrict data, int datastride, double *restrict p, int pstride, double abcdInitial[4], double abcdFinal[4])
Definition nubspline_create.cpp:396
create_NUBspline_3d_c
NUBspline_3d_c * create_NUBspline_3d_c(NUgrid *x_grid, NUgrid *y_grid, NUgrid *z_grid, BCtype_c xBC, BCtype_c yBC, BCtype_c zBC, complex_float *data)
Definition nubspline_create.cpp:839
find_NUBcoefs_1d_d
void find_NUBcoefs_1d_d(NUBasis *restrict basis, BCtype_d bc, double *data, int dstride, double *coefs, int cstride)
Definition nubspline_create.cpp:543
create_NUBspline_1d_z
NUBspline_1d_z * create_NUBspline_1d_z(NUgrid *x_grid, BCtype_z xBC, complex_double *data)
Definition nubspline_create.cpp:932
create_NUBspline_2d_d
NUBspline_2d_d * create_NUBspline_2d_d(NUgrid *x_grid, NUgrid *y_grid, BCtype_d xBC, BCtype_d yBC, double *data)
Definition nubspline_create.cpp:616
create_NUBspline_3d_d
NUBspline_3d_d * create_NUBspline_3d_d(NUgrid *x_grid, NUgrid *y_grid, NUgrid *z_grid, BCtype_d xBC, BCtype_d yBC, BCtype_d zBC, double *data)
Definition nubspline_create.cpp:668
nubspline_create.h
BCtype_c
Definition bspline_base.h:71
BCtype_c::rCode
bc_code rCode
Definition bspline_base.h:72
BCtype_c::lCode
bc_code lCode
Definition bspline_base.h:72
BCtype_c::lVal_r
float lVal_r
Definition bspline_base.h:73
BCtype_c::rVal_r
float rVal_r
Definition bspline_base.h:73
BCtype_c::rVal_i
float rVal_i
Definition bspline_base.h:73
BCtype_c::lVal_i
float lVal_i
Definition bspline_base.h:73
BCtype_d
Definition bspline_base.h:65
BCtype_d::lVal
double lVal
Definition bspline_base.h:67
BCtype_d::rCode
bc_code rCode
Definition bspline_base.h:66
BCtype_d::rVal
double rVal
Definition bspline_base.h:67
BCtype_d::lCode
bc_code lCode
Definition bspline_base.h:66
BCtype_s
Definition bspline_base.h:59
BCtype_s::lVal
float lVal
Definition bspline_base.h:61
BCtype_s::rVal
float rVal
Definition bspline_base.h:61
BCtype_s::lCode
bc_code lCode
Definition bspline_base.h:60
BCtype_s::rCode
bc_code rCode
Definition bspline_base.h:60
BCtype_z
Definition bspline_base.h:77
BCtype_z::lVal_r
double lVal_r
Definition bspline_base.h:79
BCtype_z::lVal_i
double lVal_i
Definition bspline_base.h:79
BCtype_z::lCode
bc_code lCode
Definition bspline_base.h:78
BCtype_z::rVal_i
double rVal_i
Definition bspline_base.h:79
BCtype_z::rCode
bc_code rCode
Definition bspline_base.h:78
BCtype_z::rVal_r
double rVal_r
Definition bspline_base.h:79
Bspline::sp_code
spline_code sp_code
Definition bspline_base.h:94
Bspline::coefs
void *restrict coefs
Definition bspline_base.h:96
NUBspline_1d_c
Definition nubspline_structs.h:131
NUBspline_1d_c::t_code
type_code t_code
Definition nubspline_structs.h:133
NUBspline_1d_c::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:136
NUBspline_1d_c::sp_code
spline_code sp_code
Definition nubspline_structs.h:132
NUBspline_1d_c::coefs
complex_float *restrict coefs
Definition nubspline_structs.h:134
NUBspline_1d_d
Definition nubspline_structs.h:96
NUBspline_1d_d::coefs
double *restrict coefs
Definition nubspline_structs.h:99
NUBspline_1d_d::sp_code
spline_code sp_code
Definition nubspline_structs.h:97
NUBspline_1d_d::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:101
NUBspline_1d_d::t_code
type_code t_code
Definition nubspline_structs.h:98
NUBspline_1d_s
Definition nubspline_structs.h:61
NUBspline_1d_s::t_code
type_code t_code
Definition nubspline_structs.h:63
NUBspline_1d_s::coefs
float *restrict coefs
Definition nubspline_structs.h:64
NUBspline_1d_s::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:66
NUBspline_1d_s::sp_code
spline_code sp_code
Definition nubspline_structs.h:62
NUBspline_1d_z
Definition nubspline_structs.h:166
NUBspline_1d_z::coefs
complex_double *restrict coefs
Definition nubspline_structs.h:169
NUBspline_1d_z::t_code
type_code t_code
Definition nubspline_structs.h:168
NUBspline_1d_z::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:171
NUBspline_1d_z::sp_code
spline_code sp_code
Definition nubspline_structs.h:167
NUBspline_1d
Definition nubspline_structs.h:28
NUBspline_2d_c
Definition nubspline_structs.h:141
NUBspline_2d_c::t_code
type_code t_code
Definition nubspline_structs.h:143
NUBspline_2d_c::y_basis
NUBasis *restrict *restrict y_basis
Definition nubspline_structs.h:147
NUBspline_2d_c::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:147
NUBspline_2d_c::coefs
complex_float *restrict coefs
Definition nubspline_structs.h:144
NUBspline_2d_c::x_stride
int x_stride
Definition nubspline_structs.h:145
NUBspline_2d_c::sp_code
spline_code sp_code
Definition nubspline_structs.h:142
NUBspline_2d_d
Definition nubspline_structs.h:106
NUBspline_2d_d::coefs
double *restrict coefs
Definition nubspline_structs.h:109
NUBspline_2d_d::x_stride
int x_stride
Definition nubspline_structs.h:110
NUBspline_2d_d::y_basis
NUBasis *restrict *restrict y_basis
Definition nubspline_structs.h:112
NUBspline_2d_d::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:112
NUBspline_2d_d::sp_code
spline_code sp_code
Definition nubspline_structs.h:107
NUBspline_2d_d::t_code
type_code t_code
Definition nubspline_structs.h:108
NUBspline_2d_s
Definition nubspline_structs.h:71
NUBspline_2d_s::t_code
type_code t_code
Definition nubspline_structs.h:73
NUBspline_2d_s::y_basis
NUBasis *restrict *restrict y_basis
Definition nubspline_structs.h:77
NUBspline_2d_s::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:77
NUBspline_2d_s::x_stride
int x_stride
Definition nubspline_structs.h:75
NUBspline_2d_s::coefs
float *restrict coefs
Definition nubspline_structs.h:74
NUBspline_2d_s::sp_code
spline_code sp_code
Definition nubspline_structs.h:72
NUBspline_2d_z
Definition nubspline_structs.h:176
NUBspline_2d_z::x_stride
int x_stride
Definition nubspline_structs.h:180
NUBspline_2d_z::sp_code
spline_code sp_code
Definition nubspline_structs.h:177
NUBspline_2d_z::t_code
type_code t_code
Definition nubspline_structs.h:178
NUBspline_2d_z::coefs
complex_double *restrict coefs
Definition nubspline_structs.h:179
NUBspline_2d_z::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:182
NUBspline_2d_z::y_basis
NUBasis *restrict *restrict y_basis
Definition nubspline_structs.h:182
NUBspline_2d
Definition nubspline_structs.h:37
NUBspline_3d_c
Definition nubspline_structs.h:152
NUBspline_3d_c::sp_code
spline_code sp_code
Definition nubspline_structs.h:153
NUBspline_3d_c::x_stride
int x_stride
Definition nubspline_structs.h:156
NUBspline_3d_c::y_basis
NUBasis *restrict *restrict y_basis
Definition nubspline_structs.h:158
NUBspline_3d_c::y_stride
int y_stride
Definition nubspline_structs.h:156
NUBspline_3d_c::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:158
NUBspline_3d_c::t_code
type_code t_code
Definition nubspline_structs.h:154
NUBspline_3d_c::coefs
complex_float *restrict coefs
Definition nubspline_structs.h:155
NUBspline_3d_c::z_basis
NUBasis *restrict *restrict *restrict z_basis
Definition nubspline_structs.h:158
NUBspline_3d_d
Definition nubspline_structs.h:117
NUBspline_3d_d::t_code
type_code t_code
Definition nubspline_structs.h:119
NUBspline_3d_d::y_stride
int y_stride
Definition nubspline_structs.h:121
NUBspline_3d_d::z_basis
NUBasis *restrict *restrict *restrict z_basis
Definition nubspline_structs.h:123
NUBspline_3d_d::sp_code
spline_code sp_code
Definition nubspline_structs.h:118
NUBspline_3d_d::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:123
NUBspline_3d_d::y_basis
NUBasis *restrict *restrict y_basis
Definition nubspline_structs.h:123
NUBspline_3d_d::coefs
double *restrict coefs
Definition nubspline_structs.h:120
NUBspline_3d_d::x_stride
int x_stride
Definition nubspline_structs.h:121
NUBspline_3d_s
Definition nubspline_structs.h:82
NUBspline_3d_s::sp_code
spline_code sp_code
Definition nubspline_structs.h:83
NUBspline_3d_s::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:88
NUBspline_3d_s::t_code
type_code t_code
Definition nubspline_structs.h:84
NUBspline_3d_s::y_stride
int y_stride
Definition nubspline_structs.h:86
NUBspline_3d_s::x_stride
int x_stride
Definition nubspline_structs.h:86
NUBspline_3d_s::y_basis
NUBasis *restrict *restrict y_basis
Definition nubspline_structs.h:88
NUBspline_3d_s::coefs
float *restrict coefs
Definition nubspline_structs.h:85
NUBspline_3d_s::z_basis
NUBasis *restrict *restrict *restrict z_basis
Definition nubspline_structs.h:88
NUBspline_3d_z
Definition nubspline_structs.h:187
NUBspline_3d_z::x_grid
NUgrid *restrict x_grid
Definition nubspline_structs.h:192
NUBspline_3d_z::coefs
complex_double *restrict coefs
Definition nubspline_structs.h:190
NUBspline_3d_z::y_grid
NUgrid *restrict *restrict y_grid
Definition nubspline_structs.h:192
NUBspline_3d_z::x_stride
int x_stride
Definition nubspline_structs.h:191
NUBspline_3d_z::sp_code
spline_code sp_code
Definition nubspline_structs.h:188
NUBspline_3d_z::z_basis
NUBasis *restrict *restrict *restrict z_basis
Definition nubspline_structs.h:193
NUBspline_3d_z::z_grid
NUgrid *restrict *restrict *restrict z_grid
Definition nubspline_structs.h:192
NUBspline_3d_z::x_basis
NUBasis *restrict x_basis
Definition nubspline_structs.h:193
NUBspline_3d_z::y_stride
int y_stride
Definition nubspline_structs.h:191
NUBspline_3d_z::t_code
type_code t_code
Definition nubspline_structs.h:189
NUBspline_3d_z::y_basis
NUBasis *restrict *restrict y_basis
Definition nubspline_structs.h:193
NUBspline_3d
Definition nubspline_structs.h:47
NUgrid::num_points
int num_points
Definition nugrid.h:35