Krita Source Code Documentation
Loading...
Searching...
No Matches
multi_bspline_eval_std_z.h
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#ifndef MULTI_BSPLINE_EVAL_STD_Z_IMPL_H
22#define MULTI_BSPLINE_EVAL_STD_Z_IMPL_H
23
24#include <math.h>
25#include <stdio.h>
26#include "bspline_base.h"
27#include "multi_bspline_structs.h"
28
29extern const double* restrict Ad;
30extern const double* restrict dAd;
31extern const double* restrict d2Ad;
32extern const double* restrict d3Ad;
33
34/************************************************************/
35/* 1D double-precision, complex evaulation functions */
36/************************************************************/
37void
38eval_multi_UBspline_1d_z (multi_UBspline_1d_z *spline,
39 double x,
40 complex_double* restrict vals)
41{
42 x -= spline->x_grid.start;
43 double ux = x*spline->x_grid.delta_inv;
44 double ipartx, tx;
45 tx = modf (ux, &ipartx); int ix = (int) ipartx;
46
47 double tpx[4], a[4];
48 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
49 complex_double* restrict coefs = spline->coefs;
50
51 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
52 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
53 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
54 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
55
56 intptr_t xs = spline->x_stride;
57
58 for (int n=0; n<spline->num_splines; n++)
59 vals[n] = 0.0;
60
61 for (int i=0; i<4; i++) {
62 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs);
63 for (int n=0; n<spline->num_splines; n++)
64 vals[n] += a[i] * coefs[n];
65 }
66}
67
68
69
70void
71eval_multi_UBspline_1d_z_vg (multi_UBspline_1d_z *spline,
72 double x,
73 complex_double* restrict vals,
74 complex_double* restrict grads)
75{
76 x -= spline->x_grid.start;
77 double ux = x*spline->x_grid.delta_inv;
78 double ipartx, tx;
79 tx = modf (ux, &ipartx); int ix = (int) ipartx;
80
81 double tpx[4], a[4], da[4];
82 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
83 complex_double* restrict coefs = spline->coefs;
84
85 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
86 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
87 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
88 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
89 da[0] = (dAd[ 0]*tpx[0] + dAd[ 1]*tpx[1] + dAd[ 2]*tpx[2] + dAd[ 3]*tpx[3]);
90 da[1] = (dAd[ 4]*tpx[0] + dAd[ 5]*tpx[1] + dAd[ 6]*tpx[2] + dAd[ 7]*tpx[3]);
91 da[2] = (dAd[ 8]*tpx[0] + dAd[ 9]*tpx[1] + dAd[10]*tpx[2] + dAd[11]*tpx[3]);
92 da[3] = (dAd[12]*tpx[0] + dAd[13]*tpx[1] + dAd[14]*tpx[2] + dAd[15]*tpx[3]);
93
94 intptr_t xs = spline->x_stride;
95
96 for (int n=0; n<spline->num_splines; n++) {
97 vals[n] = 0.0;
98 grads[n] = 0.0;
99 }
100
101 for (int i=0; i<4; i++) {
102 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs);
103 for (int n=0; n<spline->num_splines; n++) {
104 vals[n] += a[i] * coefs[n];
105 grads[n] += da[i] * coefs[n];
106 }
107 }
108
109 double dxInv = spline->x_grid.delta_inv;
110 for (int n=0; n<spline->num_splines; n++)
111 grads[n] *= dxInv;
112}
113
114
115void
116eval_multi_UBspline_1d_z_vgl (multi_UBspline_1d_z *spline,
117 double x,
118 complex_double* restrict vals,
119 complex_double* restrict grads,
120 complex_double* restrict lapl)
121{
122 x -= spline->x_grid.start;
123 double ux = x*spline->x_grid.delta_inv;
124 double ipartx, tx;
125 tx = modf (ux, &ipartx); int ix = (int) ipartx;
126
127 double tpx[4], a[4], da[4], d2a[4];
128 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
129 complex_double* restrict coefs = spline->coefs;
130
131 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
132 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
133 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
134 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
135 da[0] = (dAd[ 0]*tpx[0] + dAd[ 1]*tpx[1] + dAd[ 2]*tpx[2] + dAd[ 3]*tpx[3]);
136 da[1] = (dAd[ 4]*tpx[0] + dAd[ 5]*tpx[1] + dAd[ 6]*tpx[2] + dAd[ 7]*tpx[3]);
137 da[2] = (dAd[ 8]*tpx[0] + dAd[ 9]*tpx[1] + dAd[10]*tpx[2] + dAd[11]*tpx[3]);
138 da[3] = (dAd[12]*tpx[0] + dAd[13]*tpx[1] + dAd[14]*tpx[2] + dAd[15]*tpx[3]);
139 d2a[0] = (d2Ad[ 0]*tpx[0] + d2Ad[ 1]*tpx[1] + d2Ad[ 2]*tpx[2] + d2Ad[ 3]*tpx[3]);
140 d2a[1] = (d2Ad[ 4]*tpx[0] + d2Ad[ 5]*tpx[1] + d2Ad[ 6]*tpx[2] + d2Ad[ 7]*tpx[3]);
141 d2a[2] = (d2Ad[ 8]*tpx[0] + d2Ad[ 9]*tpx[1] + d2Ad[10]*tpx[2] + d2Ad[11]*tpx[3]);
142 d2a[3] = (d2Ad[12]*tpx[0] + d2Ad[13]*tpx[1] + d2Ad[14]*tpx[2] + d2Ad[15]*tpx[3]);
143
144 intptr_t xs = spline->x_stride;
145
146 for (int n=0; n<spline->num_splines; n++) {
147 vals[n] = 0.0;
148 grads[n] = 0.0;
149 lapl[n] = 0.0;
150 }
151
152 for (int i=0; i<4; i++) {
153 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs);
154 for (int n=0; n<spline->num_splines; n++) {
155 vals[n] += a[i] * coefs[n];
156 grads[n] += da[i] * coefs[n];
157 lapl[n] += d2a[i] * coefs[n];
158 }
159 }
160
161 double dxInv = spline->x_grid.delta_inv;
162 for (int n=0; n<spline->num_splines; n++) {
163 grads[n] *= dxInv;
164 lapl [n] *= dxInv*dxInv;
165 }
166}
167
168
169void
170eval_multi_UBspline_1d_z_vgh (multi_UBspline_1d_z *spline,
171 double x,
172 complex_double* restrict vals,
173 complex_double* restrict grads,
174 complex_double* restrict hess)
175{
176 eval_multi_UBspline_1d_z_vgl (spline, x, vals, grads, hess);
177}
178
179
180/************************************************************/
181/* 2D double-precision, complex evaulation functions */
182/************************************************************/
183void
184eval_multi_UBspline_2d_z (multi_UBspline_2d_z *spline,
185 double x, double y,
186 complex_double* restrict vals)
187{
188 x -= spline->x_grid.start;
189 y -= spline->y_grid.start;
190 double ux = x*spline->x_grid.delta_inv;
191 double uy = y*spline->y_grid.delta_inv;
192 double ipartx, iparty, tx, ty;
193 tx = modf (ux, &ipartx); int ix = (int) ipartx;
194 ty = modf (uy, &iparty); int iy = (int) iparty;
195
196 double tpx[4], tpy[4], a[4], b[4];
197 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
198 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
199 complex_double* restrict coefs = spline->coefs;
200
201 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
202 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
203 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
204 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
205
206 b[0] = (Ad[ 0]*tpy[0] + Ad[ 1]*tpy[1] + Ad[ 2]*tpy[2] + Ad[ 3]*tpy[3]);
207 b[1] = (Ad[ 4]*tpy[0] + Ad[ 5]*tpy[1] + Ad[ 6]*tpy[2] + Ad[ 7]*tpy[3]);
208 b[2] = (Ad[ 8]*tpy[0] + Ad[ 9]*tpy[1] + Ad[10]*tpy[2] + Ad[11]*tpy[3]);
209 b[3] = (Ad[12]*tpy[0] + Ad[13]*tpy[1] + Ad[14]*tpy[2] + Ad[15]*tpy[3]);
210
211 intptr_t xs = spline->x_stride;
212 intptr_t ys = spline->y_stride;
213
214 for (int n=0; n<spline->num_splines; n++)
215 vals[n] = 0.0;
216
217 for (int i=0; i<4; i++)
218 for (int j=0; j<4; j++) {
219 double prefactor = a[i]*b[j];
220 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys);
221 for (int n=0; n<spline->num_splines; n++)
222 vals[n] += prefactor*coefs[n];
223 }
224}
225
226
227void
228eval_multi_UBspline_2d_z_vg (multi_UBspline_2d_z *spline,
229 double x, double y,
230 complex_double* restrict vals,
231 complex_double* restrict grads)
232{
233 x -= spline->x_grid.start;
234 y -= spline->y_grid.start;
235 double ux = x*spline->x_grid.delta_inv;
236 double uy = y*spline->y_grid.delta_inv;
237 double ipartx, iparty, tx, ty;
238 tx = modf (ux, &ipartx); int ix = (int) ipartx;
239 ty = modf (uy, &iparty); int iy = (int) iparty;
240
241 double tpx[4], tpy[4], a[4], b[4], da[4], db[4];
242 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
243 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
244 complex_double* restrict coefs = spline->coefs;
245
246 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
247 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
248 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
249 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
250 da[0] = (dAd[ 0]*tpx[0] + dAd[ 1]*tpx[1] + dAd[ 2]*tpx[2] + dAd[ 3]*tpx[3]);
251 da[1] = (dAd[ 4]*tpx[0] + dAd[ 5]*tpx[1] + dAd[ 6]*tpx[2] + dAd[ 7]*tpx[3]);
252 da[2] = (dAd[ 8]*tpx[0] + dAd[ 9]*tpx[1] + dAd[10]*tpx[2] + dAd[11]*tpx[3]);
253 da[3] = (dAd[12]*tpx[0] + dAd[13]*tpx[1] + dAd[14]*tpx[2] + dAd[15]*tpx[3]);
254
255 b[0] = (Ad[ 0]*tpy[0] + Ad[ 1]*tpy[1] + Ad[ 2]*tpy[2] + Ad[ 3]*tpy[3]);
256 b[1] = (Ad[ 4]*tpy[0] + Ad[ 5]*tpy[1] + Ad[ 6]*tpy[2] + Ad[ 7]*tpy[3]);
257 b[2] = (Ad[ 8]*tpy[0] + Ad[ 9]*tpy[1] + Ad[10]*tpy[2] + Ad[11]*tpy[3]);
258 b[3] = (Ad[12]*tpy[0] + Ad[13]*tpy[1] + Ad[14]*tpy[2] + Ad[15]*tpy[3]);
259 db[0] = (dAd[ 0]*tpy[0] + dAd[ 1]*tpy[1] + dAd[ 2]*tpy[2] + dAd[ 3]*tpy[3]);
260 db[1] = (dAd[ 4]*tpy[0] + dAd[ 5]*tpy[1] + dAd[ 6]*tpy[2] + dAd[ 7]*tpy[3]);
261 db[2] = (dAd[ 8]*tpy[0] + dAd[ 9]*tpy[1] + dAd[10]*tpy[2] + dAd[11]*tpy[3]);
262 db[3] = (dAd[12]*tpy[0] + dAd[13]*tpy[1] + dAd[14]*tpy[2] + dAd[15]*tpy[3]);
263
264 intptr_t xs = spline->x_stride;
265 intptr_t ys = spline->y_stride;
266
267 for (int n=0; n<spline->num_splines; n++) {
268 vals[n] = 0.0;
269 grads[2*n+0] = grads[2*n+1] = grads[2*n+2] = 0.0;
270 }
271
272 for (int i=0; i<4; i++)
273 for (int j=0; j<4; j++) {
274 double ab = a[i]*b[j];
275 double dab[2];
276 dab[0] = da[i]* b[j];
277 dab[1] = a[i]*db[j];
278
279 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys);
280 for (int n=0; n<spline->num_splines; n++) {
281 vals [n] += ab *coefs[n];
282 grads[2*n+0] += dab[0]*coefs[n];
283 grads[2*n+1] += dab[1]*coefs[n];
284 }
285 }
286
287 double dxInv = spline->x_grid.delta_inv;
288 double dyInv = spline->y_grid.delta_inv;
289 for (int n=0; n<spline->num_splines; n++) {
290 grads[2*n+0] *= dxInv;
291 grads[2*n+1] *= dyInv;
292 }
293}
294
295void
296eval_multi_UBspline_2d_z_vgl (multi_UBspline_2d_z *spline,
297 double x, double y,
298 complex_double* restrict vals,
299 complex_double* restrict grads,
300 complex_double* restrict lapl)
301{
302 x -= spline->x_grid.start;
303 y -= spline->y_grid.start;
304 double ux = x*spline->x_grid.delta_inv;
305 double uy = y*spline->y_grid.delta_inv;
306 double ipartx, iparty, tx, ty;
307 tx = modf (ux, &ipartx); int ix = (int) ipartx;
308 ty = modf (uy, &iparty); int iy = (int) iparty;
309
310 double tpx[4], tpy[4], a[4], b[4], da[4], db[4], d2a[4], d2b[4];
311 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
312 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
313 complex_double* restrict coefs = spline->coefs;
314
315 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
316 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
317 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
318 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
319 da[0] = (dAd[ 0]*tpx[0] + dAd[ 1]*tpx[1] + dAd[ 2]*tpx[2] + dAd[ 3]*tpx[3]);
320 da[1] = (dAd[ 4]*tpx[0] + dAd[ 5]*tpx[1] + dAd[ 6]*tpx[2] + dAd[ 7]*tpx[3]);
321 da[2] = (dAd[ 8]*tpx[0] + dAd[ 9]*tpx[1] + dAd[10]*tpx[2] + dAd[11]*tpx[3]);
322 da[3] = (dAd[12]*tpx[0] + dAd[13]*tpx[1] + dAd[14]*tpx[2] + dAd[15]*tpx[3]);
323 d2a[0] = (d2Ad[ 0]*tpx[0] + d2Ad[ 1]*tpx[1] + d2Ad[ 2]*tpx[2] + d2Ad[ 3]*tpx[3]);
324 d2a[1] = (d2Ad[ 4]*tpx[0] + d2Ad[ 5]*tpx[1] + d2Ad[ 6]*tpx[2] + d2Ad[ 7]*tpx[3]);
325 d2a[2] = (d2Ad[ 8]*tpx[0] + d2Ad[ 9]*tpx[1] + d2Ad[10]*tpx[2] + d2Ad[11]*tpx[3]);
326 d2a[3] = (d2Ad[12]*tpx[0] + d2Ad[13]*tpx[1] + d2Ad[14]*tpx[2] + d2Ad[15]*tpx[3]);
327
328 b[0] = (Ad[ 0]*tpy[0] + Ad[ 1]*tpy[1] + Ad[ 2]*tpy[2] + Ad[ 3]*tpy[3]);
329 b[1] = (Ad[ 4]*tpy[0] + Ad[ 5]*tpy[1] + Ad[ 6]*tpy[2] + Ad[ 7]*tpy[3]);
330 b[2] = (Ad[ 8]*tpy[0] + Ad[ 9]*tpy[1] + Ad[10]*tpy[2] + Ad[11]*tpy[3]);
331 b[3] = (Ad[12]*tpy[0] + Ad[13]*tpy[1] + Ad[14]*tpy[2] + Ad[15]*tpy[3]);
332 db[0] = (dAd[ 0]*tpy[0] + dAd[ 1]*tpy[1] + dAd[ 2]*tpy[2] + dAd[ 3]*tpy[3]);
333 db[1] = (dAd[ 4]*tpy[0] + dAd[ 5]*tpy[1] + dAd[ 6]*tpy[2] + dAd[ 7]*tpy[3]);
334 db[2] = (dAd[ 8]*tpy[0] + dAd[ 9]*tpy[1] + dAd[10]*tpy[2] + dAd[11]*tpy[3]);
335 db[3] = (dAd[12]*tpy[0] + dAd[13]*tpy[1] + dAd[14]*tpy[2] + dAd[15]*tpy[3]);
336 d2b[0] = (d2Ad[ 0]*tpy[0] + d2Ad[ 1]*tpy[1] + d2Ad[ 2]*tpy[2] + d2Ad[ 3]*tpy[3]);
337 d2b[1] = (d2Ad[ 4]*tpy[0] + d2Ad[ 5]*tpy[1] + d2Ad[ 6]*tpy[2] + d2Ad[ 7]*tpy[3]);
338 d2b[2] = (d2Ad[ 8]*tpy[0] + d2Ad[ 9]*tpy[1] + d2Ad[10]*tpy[2] + d2Ad[11]*tpy[3]);
339 d2b[3] = (d2Ad[12]*tpy[0] + d2Ad[13]*tpy[1] + d2Ad[14]*tpy[2] + d2Ad[15]*tpy[3]);
340
341 intptr_t xs = spline->x_stride;
342 intptr_t ys = spline->y_stride;
343
344 //complex_double lapl2[2*spline->num_splines];
345 complex_double* restrict lapl2 = spline->lapl2;
346 for (int n=0; n<spline->num_splines; n++) {
347 vals[n] = 0.0;
348 grads[2*n+0] = grads[2*n+1] = 0.0;
349 lapl2[2*n+0] = lapl2[2*n+1] = 0.0;
350 }
351
352 for (int i=0; i<4; i++)
353 for (int j=0; j<4; j++) {
354 double ab = a[i]*b[j];
355 double dab[2], d2ab[2];
356 dab[0] = da[i]* b[j];
357 dab[1] = a[i]*db[j];
358 d2ab[0] = d2a[i]* b[j];
359 d2ab[1] = a[i]*d2b[j];
360
361 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys);
362 for (int n=0; n<spline->num_splines; n++) {
363 vals[n] += ab *coefs[n];
364 grads[2*n+0] += dab[0]*coefs[n];
365 grads[2*n+1] += dab[1]*coefs[n];
366 lapl2[2*n+0] += d2ab[0]*coefs[n];
367 lapl2[2*n+1] += d2ab[1]*coefs[n];
368 }
369 }
370
371 double dxInv = spline->x_grid.delta_inv;
372 double dyInv = spline->y_grid.delta_inv;
373 for (int n=0; n<spline->num_splines; n++) {
374 grads[2*n+0] *= dxInv;
375 grads[2*n+1] *= dyInv;
376 lapl2[2*n+0] *= dxInv*dxInv;
377 lapl2[2*n+1] *= dyInv*dyInv;
378 lapl[n] = lapl2[2*n+0] + lapl2[2*n+1];
379 }
380}
381
382
383
384
385void
386eval_multi_UBspline_2d_z_vgh (multi_UBspline_2d_z *spline,
387 double x, double y,
388 complex_double* restrict vals,
389 complex_double* restrict grads,
390 complex_double* restrict hess)
391{
392 x -= spline->x_grid.start;
393 y -= spline->y_grid.start;
394 double ux = x*spline->x_grid.delta_inv;
395 double uy = y*spline->y_grid.delta_inv;
396 double ipartx, iparty, tx, ty;
397 tx = modf (ux, &ipartx); int ix = (int) ipartx;
398 ty = modf (uy, &iparty); int iy = (int) iparty;
399
400 double tpx[4], tpy[4], a[4], b[4], da[4], db[4], d2a[4], d2b[4];
401 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
402 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
403 complex_double* restrict coefs = spline->coefs;
404
405 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
406 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
407 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
408 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
409 da[0] = (dAd[ 0]*tpx[0] + dAd[ 1]*tpx[1] + dAd[ 2]*tpx[2] + dAd[ 3]*tpx[3]);
410 da[1] = (dAd[ 4]*tpx[0] + dAd[ 5]*tpx[1] + dAd[ 6]*tpx[2] + dAd[ 7]*tpx[3]);
411 da[2] = (dAd[ 8]*tpx[0] + dAd[ 9]*tpx[1] + dAd[10]*tpx[2] + dAd[11]*tpx[3]);
412 da[3] = (dAd[12]*tpx[0] + dAd[13]*tpx[1] + dAd[14]*tpx[2] + dAd[15]*tpx[3]);
413 d2a[0] = (d2Ad[ 0]*tpx[0] + d2Ad[ 1]*tpx[1] + d2Ad[ 2]*tpx[2] + d2Ad[ 3]*tpx[3]);
414 d2a[1] = (d2Ad[ 4]*tpx[0] + d2Ad[ 5]*tpx[1] + d2Ad[ 6]*tpx[2] + d2Ad[ 7]*tpx[3]);
415 d2a[2] = (d2Ad[ 8]*tpx[0] + d2Ad[ 9]*tpx[1] + d2Ad[10]*tpx[2] + d2Ad[11]*tpx[3]);
416 d2a[3] = (d2Ad[12]*tpx[0] + d2Ad[13]*tpx[1] + d2Ad[14]*tpx[2] + d2Ad[15]*tpx[3]);
417
418 b[0] = (Ad[ 0]*tpy[0] + Ad[ 1]*tpy[1] + Ad[ 2]*tpy[2] + Ad[ 3]*tpy[3]);
419 b[1] = (Ad[ 4]*tpy[0] + Ad[ 5]*tpy[1] + Ad[ 6]*tpy[2] + Ad[ 7]*tpy[3]);
420 b[2] = (Ad[ 8]*tpy[0] + Ad[ 9]*tpy[1] + Ad[10]*tpy[2] + Ad[11]*tpy[3]);
421 b[3] = (Ad[12]*tpy[0] + Ad[13]*tpy[1] + Ad[14]*tpy[2] + Ad[15]*tpy[3]);
422 db[0] = (dAd[ 0]*tpy[0] + dAd[ 1]*tpy[1] + dAd[ 2]*tpy[2] + dAd[ 3]*tpy[3]);
423 db[1] = (dAd[ 4]*tpy[0] + dAd[ 5]*tpy[1] + dAd[ 6]*tpy[2] + dAd[ 7]*tpy[3]);
424 db[2] = (dAd[ 8]*tpy[0] + dAd[ 9]*tpy[1] + dAd[10]*tpy[2] + dAd[11]*tpy[3]);
425 db[3] = (dAd[12]*tpy[0] + dAd[13]*tpy[1] + dAd[14]*tpy[2] + dAd[15]*tpy[3]);
426 d2b[0] = (d2Ad[ 0]*tpy[0] + d2Ad[ 1]*tpy[1] + d2Ad[ 2]*tpy[2] + d2Ad[ 3]*tpy[3]);
427 d2b[1] = (d2Ad[ 4]*tpy[0] + d2Ad[ 5]*tpy[1] + d2Ad[ 6]*tpy[2] + d2Ad[ 7]*tpy[3]);
428 d2b[2] = (d2Ad[ 8]*tpy[0] + d2Ad[ 9]*tpy[1] + d2Ad[10]*tpy[2] + d2Ad[11]*tpy[3]);
429 d2b[3] = (d2Ad[12]*tpy[0] + d2Ad[13]*tpy[1] + d2Ad[14]*tpy[2] + d2Ad[15]*tpy[3]);
430
431 intptr_t xs = spline->x_stride;
432 intptr_t ys = spline->y_stride;
433
434 for (int n=0; n<spline->num_splines; n++) {
435 vals[n] = 0.0;
436 grads[2*n+0] = grads[2*n+1] = 0.0;
437 for (int i=0; i<4; i++)
438 hess[4*n+i] = 0.0;
439 }
440
441 for (int i=0; i<4; i++)
442 for (int j=0; j<4; j++){
443 double ab = a[i]*b[j];
444 double dab[2], d2ab[3];
445 dab[0] = da[i]* b[j];
446 dab[1] = a[i]*db[j];
447 d2ab[0] = d2a[i] * b[j];
448 d2ab[1] = da[i] * db[j];
449 d2ab[2] = a[i] * d2b[j];
450
451 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys);
452 for (int n=0; n<spline->num_splines; n++) {
453 vals[n] += ab *coefs[n];
454 grads[2*n+0] += dab[0]*coefs[n];
455 grads[2*n+1] += dab[1]*coefs[n];
456 hess [4*n+0] += d2ab[0]*coefs[n];
457 hess [4*n+1] += d2ab[1]*coefs[n];
458 hess [4*n+3] += d2ab[2]*coefs[n];
459 }
460 }
461
462 double dxInv = spline->x_grid.delta_inv;
463 double dyInv = spline->y_grid.delta_inv;
464 for (int n=0; n<spline->num_splines; n++) {
465 grads[2*n+0] *= dxInv;
466 grads[2*n+1] *= dyInv;
467 hess[4*n+0] *= dxInv*dxInv;
468 hess[4*n+1] *= dxInv*dyInv;
469 hess[4*n+3] *= dyInv*dyInv;
470 // Copy hessian elements into lower half of 3x3 matrix
471 hess[4*n+2] = hess[4*n+1];
472 }
473}
474
475
476
477/************************************************************/
478/* 3D double-precision, complex evaulation functions */
479/************************************************************/
480void
481eval_multi_UBspline_3d_z (multi_UBspline_3d_z *spline,
482 double x, double y, double z,
483 complex_double* restrict vals)
484{
485 x -= spline->x_grid.start;
486 y -= spline->y_grid.start;
487 z -= spline->z_grid.start;
488 double ux = x*spline->x_grid.delta_inv;
489 double uy = y*spline->y_grid.delta_inv;
490 double uz = z*spline->z_grid.delta_inv;
491 double ipartx, iparty, ipartz, tx, ty, tz;
492 tx = modf (ux, &ipartx); int ix = (int) ipartx;
493 ty = modf (uy, &iparty); int iy = (int) iparty;
494 tz = modf (uz, &ipartz); int iz = (int) ipartz;
495
496 double tpx[4], tpy[4], tpz[4], a[4], b[4], c[4];
497 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
498 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
499 tpz[0] = tz*tz*tz; tpz[1] = tz*tz; tpz[2] = tz; tpz[3] = 1.0;
500 complex_double* restrict coefs = spline->coefs;
501
502 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
503 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
504 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
505 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
506
507 b[0] = (Ad[ 0]*tpy[0] + Ad[ 1]*tpy[1] + Ad[ 2]*tpy[2] + Ad[ 3]*tpy[3]);
508 b[1] = (Ad[ 4]*tpy[0] + Ad[ 5]*tpy[1] + Ad[ 6]*tpy[2] + Ad[ 7]*tpy[3]);
509 b[2] = (Ad[ 8]*tpy[0] + Ad[ 9]*tpy[1] + Ad[10]*tpy[2] + Ad[11]*tpy[3]);
510 b[3] = (Ad[12]*tpy[0] + Ad[13]*tpy[1] + Ad[14]*tpy[2] + Ad[15]*tpy[3]);
511
512 c[0] = (Ad[ 0]*tpz[0] + Ad[ 1]*tpz[1] + Ad[ 2]*tpz[2] + Ad[ 3]*tpz[3]);
513 c[1] = (Ad[ 4]*tpz[0] + Ad[ 5]*tpz[1] + Ad[ 6]*tpz[2] + Ad[ 7]*tpz[3]);
514 c[2] = (Ad[ 8]*tpz[0] + Ad[ 9]*tpz[1] + Ad[10]*tpz[2] + Ad[11]*tpz[3]);
515 c[3] = (Ad[12]*tpz[0] + Ad[13]*tpz[1] + Ad[14]*tpz[2] + Ad[15]*tpz[3]);
516
517 intptr_t xs = spline->x_stride;
518 intptr_t ys = spline->y_stride;
519 intptr_t zs = spline->z_stride;
520
521 for (int n=0; n<spline->num_splines; n++)
522 vals[n] = 0.0;
523
524 for (int i=0; i<4; i++)
525 for (int j=0; j<4; j++)
526 for (int k=0; k<4; k++) {
527 double prefactor = a[i]*b[j]*c[k];
528 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys + (iz+k)*zs);
529 for (int n=0; n<spline->num_splines; n++)
530 vals[n] += prefactor*coefs[n];
531 }
532}
533
534
535void
536eval_multi_UBspline_3d_z_vg (multi_UBspline_3d_z *spline,
537 double x, double y, double z,
538 complex_double* restrict vals,
539 complex_double* restrict grads)
540{
541 x -= spline->x_grid.start;
542 y -= spline->y_grid.start;
543 z -= spline->z_grid.start;
544 double ux = x*spline->x_grid.delta_inv;
545 double uy = y*spline->y_grid.delta_inv;
546 double uz = z*spline->z_grid.delta_inv;
547 double ipartx, iparty, ipartz, tx, ty, tz;
548 tx = modf (ux, &ipartx); int ix = (int) ipartx;
549 ty = modf (uy, &iparty); int iy = (int) iparty;
550 tz = modf (uz, &ipartz); int iz = (int) ipartz;
551
552 double tpx[4], tpy[4], tpz[4], a[4], b[4], c[4],
553 da[4], db[4], dc[4];
554 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
555 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
556 tpz[0] = tz*tz*tz; tpz[1] = tz*tz; tpz[2] = tz; tpz[3] = 1.0;
557 complex_double* restrict coefs = spline->coefs;
558
559 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
560 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
561 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
562 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
563 da[0] = (dAd[ 0]*tpx[0] + dAd[ 1]*tpx[1] + dAd[ 2]*tpx[2] + dAd[ 3]*tpx[3]);
564 da[1] = (dAd[ 4]*tpx[0] + dAd[ 5]*tpx[1] + dAd[ 6]*tpx[2] + dAd[ 7]*tpx[3]);
565 da[2] = (dAd[ 8]*tpx[0] + dAd[ 9]*tpx[1] + dAd[10]*tpx[2] + dAd[11]*tpx[3]);
566 da[3] = (dAd[12]*tpx[0] + dAd[13]*tpx[1] + dAd[14]*tpx[2] + dAd[15]*tpx[3]);
567
568 b[0] = (Ad[ 0]*tpy[0] + Ad[ 1]*tpy[1] + Ad[ 2]*tpy[2] + Ad[ 3]*tpy[3]);
569 b[1] = (Ad[ 4]*tpy[0] + Ad[ 5]*tpy[1] + Ad[ 6]*tpy[2] + Ad[ 7]*tpy[3]);
570 b[2] = (Ad[ 8]*tpy[0] + Ad[ 9]*tpy[1] + Ad[10]*tpy[2] + Ad[11]*tpy[3]);
571 b[3] = (Ad[12]*tpy[0] + Ad[13]*tpy[1] + Ad[14]*tpy[2] + Ad[15]*tpy[3]);
572 db[0] = (dAd[ 0]*tpy[0] + dAd[ 1]*tpy[1] + dAd[ 2]*tpy[2] + dAd[ 3]*tpy[3]);
573 db[1] = (dAd[ 4]*tpy[0] + dAd[ 5]*tpy[1] + dAd[ 6]*tpy[2] + dAd[ 7]*tpy[3]);
574 db[2] = (dAd[ 8]*tpy[0] + dAd[ 9]*tpy[1] + dAd[10]*tpy[2] + dAd[11]*tpy[3]);
575 db[3] = (dAd[12]*tpy[0] + dAd[13]*tpy[1] + dAd[14]*tpy[2] + dAd[15]*tpy[3]);
576
577 c[0] = (Ad[ 0]*tpz[0] + Ad[ 1]*tpz[1] + Ad[ 2]*tpz[2] + Ad[ 3]*tpz[3]);
578 c[1] = (Ad[ 4]*tpz[0] + Ad[ 5]*tpz[1] + Ad[ 6]*tpz[2] + Ad[ 7]*tpz[3]);
579 c[2] = (Ad[ 8]*tpz[0] + Ad[ 9]*tpz[1] + Ad[10]*tpz[2] + Ad[11]*tpz[3]);
580 c[3] = (Ad[12]*tpz[0] + Ad[13]*tpz[1] + Ad[14]*tpz[2] + Ad[15]*tpz[3]);
581 dc[0] = (dAd[ 0]*tpz[0] + dAd[ 1]*tpz[1] + dAd[ 2]*tpz[2] + dAd[ 3]*tpz[3]);
582 dc[1] = (dAd[ 4]*tpz[0] + dAd[ 5]*tpz[1] + dAd[ 6]*tpz[2] + dAd[ 7]*tpz[3]);
583 dc[2] = (dAd[ 8]*tpz[0] + dAd[ 9]*tpz[1] + dAd[10]*tpz[2] + dAd[11]*tpz[3]);
584 dc[3] = (dAd[12]*tpz[0] + dAd[13]*tpz[1] + dAd[14]*tpz[2] + dAd[15]*tpz[3]);
585
586 intptr_t xs = spline->x_stride;
587 intptr_t ys = spline->y_stride;
588 intptr_t zs = spline->z_stride;
589
590 for (int n=0; n<spline->num_splines; n++) {
591 vals[n] = 0.0;
592 grads[3*n+0] = grads[3*n+1] = grads[3*n+2] = 0.0;
593 }
594
595 for (int i=0; i<4; i++)
596 for (int j=0; j<4; j++)
597 for (int k=0; k<4; k++) {
598 double abc = a[i]*b[j]*c[k];
599 double dabc[3];
600 dabc[0] = da[i]* b[j]* c[k];
601 dabc[1] = a[i]*db[j]* c[k];
602 dabc[2] = a[i]* b[j]*dc[k];
603
604 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys + (iz+k)*zs);
605 for (int n=0; n<spline->num_splines; n++) {
606 vals[n] += abc *coefs[n];
607 grads[3*n+0] += dabc[0]*coefs[n];
608 grads[3*n+1] += dabc[1]*coefs[n];
609 grads[3*n+2] += dabc[2]*coefs[n];
610 }
611 }
612
613 double dxInv = spline->x_grid.delta_inv;
614 double dyInv = spline->y_grid.delta_inv;
615 double dzInv = spline->z_grid.delta_inv;
616 for (int n=0; n<spline->num_splines; n++) {
617 grads[3*n+0] *= dxInv;
618 grads[3*n+1] *= dyInv;
619 grads[3*n+2] *= dzInv;
620 }
621}
622
623
624
625void
626eval_multi_UBspline_3d_z_vgl (multi_UBspline_3d_z *spline,
627 double x, double y, double z,
628 complex_double* restrict vals,
629 complex_double* restrict grads,
630 complex_double* restrict lapl)
631{
632 x -= spline->x_grid.start;
633 y -= spline->y_grid.start;
634 z -= spline->z_grid.start;
635 double ux = x*spline->x_grid.delta_inv;
636 double uy = y*spline->y_grid.delta_inv;
637 double uz = z*spline->z_grid.delta_inv;
638 double ipartx, iparty, ipartz, tx, ty, tz;
639 tx = modf (ux, &ipartx); int ix = (int) ipartx;
640 ty = modf (uy, &iparty); int iy = (int) iparty;
641 tz = modf (uz, &ipartz); int iz = (int) ipartz;
642
643 double tpx[4], tpy[4], tpz[4], a[4], b[4], c[4],
644 da[4], db[4], dc[4], d2a[4], d2b[4], d2c[4];
645 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
646 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
647 tpz[0] = tz*tz*tz; tpz[1] = tz*tz; tpz[2] = tz; tpz[3] = 1.0;
648 complex_double* restrict coefs = spline->coefs;
649
650 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
651 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
652 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
653 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
654 da[0] = (dAd[ 0]*tpx[0] + dAd[ 1]*tpx[1] + dAd[ 2]*tpx[2] + dAd[ 3]*tpx[3]);
655 da[1] = (dAd[ 4]*tpx[0] + dAd[ 5]*tpx[1] + dAd[ 6]*tpx[2] + dAd[ 7]*tpx[3]);
656 da[2] = (dAd[ 8]*tpx[0] + dAd[ 9]*tpx[1] + dAd[10]*tpx[2] + dAd[11]*tpx[3]);
657 da[3] = (dAd[12]*tpx[0] + dAd[13]*tpx[1] + dAd[14]*tpx[2] + dAd[15]*tpx[3]);
658 d2a[0] = (d2Ad[ 0]*tpx[0] + d2Ad[ 1]*tpx[1] + d2Ad[ 2]*tpx[2] + d2Ad[ 3]*tpx[3]);
659 d2a[1] = (d2Ad[ 4]*tpx[0] + d2Ad[ 5]*tpx[1] + d2Ad[ 6]*tpx[2] + d2Ad[ 7]*tpx[3]);
660 d2a[2] = (d2Ad[ 8]*tpx[0] + d2Ad[ 9]*tpx[1] + d2Ad[10]*tpx[2] + d2Ad[11]*tpx[3]);
661 d2a[3] = (d2Ad[12]*tpx[0] + d2Ad[13]*tpx[1] + d2Ad[14]*tpx[2] + d2Ad[15]*tpx[3]);
662
663 b[0] = (Ad[ 0]*tpy[0] + Ad[ 1]*tpy[1] + Ad[ 2]*tpy[2] + Ad[ 3]*tpy[3]);
664 b[1] = (Ad[ 4]*tpy[0] + Ad[ 5]*tpy[1] + Ad[ 6]*tpy[2] + Ad[ 7]*tpy[3]);
665 b[2] = (Ad[ 8]*tpy[0] + Ad[ 9]*tpy[1] + Ad[10]*tpy[2] + Ad[11]*tpy[3]);
666 b[3] = (Ad[12]*tpy[0] + Ad[13]*tpy[1] + Ad[14]*tpy[2] + Ad[15]*tpy[3]);
667 db[0] = (dAd[ 0]*tpy[0] + dAd[ 1]*tpy[1] + dAd[ 2]*tpy[2] + dAd[ 3]*tpy[3]);
668 db[1] = (dAd[ 4]*tpy[0] + dAd[ 5]*tpy[1] + dAd[ 6]*tpy[2] + dAd[ 7]*tpy[3]);
669 db[2] = (dAd[ 8]*tpy[0] + dAd[ 9]*tpy[1] + dAd[10]*tpy[2] + dAd[11]*tpy[3]);
670 db[3] = (dAd[12]*tpy[0] + dAd[13]*tpy[1] + dAd[14]*tpy[2] + dAd[15]*tpy[3]);
671 d2b[0] = (d2Ad[ 0]*tpy[0] + d2Ad[ 1]*tpy[1] + d2Ad[ 2]*tpy[2] + d2Ad[ 3]*tpy[3]);
672 d2b[1] = (d2Ad[ 4]*tpy[0] + d2Ad[ 5]*tpy[1] + d2Ad[ 6]*tpy[2] + d2Ad[ 7]*tpy[3]);
673 d2b[2] = (d2Ad[ 8]*tpy[0] + d2Ad[ 9]*tpy[1] + d2Ad[10]*tpy[2] + d2Ad[11]*tpy[3]);
674 d2b[3] = (d2Ad[12]*tpy[0] + d2Ad[13]*tpy[1] + d2Ad[14]*tpy[2] + d2Ad[15]*tpy[3]);
675
676 c[0] = (Ad[ 0]*tpz[0] + Ad[ 1]*tpz[1] + Ad[ 2]*tpz[2] + Ad[ 3]*tpz[3]);
677 c[1] = (Ad[ 4]*tpz[0] + Ad[ 5]*tpz[1] + Ad[ 6]*tpz[2] + Ad[ 7]*tpz[3]);
678 c[2] = (Ad[ 8]*tpz[0] + Ad[ 9]*tpz[1] + Ad[10]*tpz[2] + Ad[11]*tpz[3]);
679 c[3] = (Ad[12]*tpz[0] + Ad[13]*tpz[1] + Ad[14]*tpz[2] + Ad[15]*tpz[3]);
680 dc[0] = (dAd[ 0]*tpz[0] + dAd[ 1]*tpz[1] + dAd[ 2]*tpz[2] + dAd[ 3]*tpz[3]);
681 dc[1] = (dAd[ 4]*tpz[0] + dAd[ 5]*tpz[1] + dAd[ 6]*tpz[2] + dAd[ 7]*tpz[3]);
682 dc[2] = (dAd[ 8]*tpz[0] + dAd[ 9]*tpz[1] + dAd[10]*tpz[2] + dAd[11]*tpz[3]);
683 dc[3] = (dAd[12]*tpz[0] + dAd[13]*tpz[1] + dAd[14]*tpz[2] + dAd[15]*tpz[3]);
684 d2c[0] = (d2Ad[ 0]*tpz[0] + d2Ad[ 1]*tpz[1] + d2Ad[ 2]*tpz[2] + d2Ad[ 3]*tpz[3]);
685 d2c[1] = (d2Ad[ 4]*tpz[0] + d2Ad[ 5]*tpz[1] + d2Ad[ 6]*tpz[2] + d2Ad[ 7]*tpz[3]);
686 d2c[2] = (d2Ad[ 8]*tpz[0] + d2Ad[ 9]*tpz[1] + d2Ad[10]*tpz[2] + d2Ad[11]*tpz[3]);
687 d2c[3] = (d2Ad[12]*tpz[0] + d2Ad[13]*tpz[1] + d2Ad[14]*tpz[2] + d2Ad[15]*tpz[3]);
688
689 intptr_t xs = spline->x_stride;
690 intptr_t ys = spline->y_stride;
691 intptr_t zs = spline->z_stride;
692
693 //complex_double lapl3[3*spline->num_splines];
694 complex_double* restrict lapl3 = spline->lapl3;
695 for (int n=0; n<spline->num_splines; n++) {
696 vals[n] = 0.0;
697 grads[3*n+0] = grads[3*n+1] = grads[3*n+2] = 0.0;
698 lapl3[3*n+0] = lapl3[3*n+1] = lapl3[3*n+2] = 0.0;
699 }
700
701 for (int i=0; i<4; i++)
702 for (int j=0; j<4; j++)
703 for (int k=0; k<4; k++) {
704 double abc = a[i]*b[j]*c[k];
705 double dabc[3], d2abc[3];
706 dabc[0] = da[i]* b[j]* c[k];
707 dabc[1] = a[i]*db[j]* c[k];
708 dabc[2] = a[i]* b[j]*dc[k];
709 d2abc[0] = d2a[i]* b[j]* c[k];
710 d2abc[1] = a[i]*d2b[j]* c[k];
711 d2abc[2] = a[i]* b[j]*d2c[k];
712
713 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys + (iz+k)*zs);
714 for (int n=0; n<spline->num_splines; n++) {
715 vals[n] += abc *coefs[n];
716 grads[3*n+0] += dabc[0]*coefs[n];
717 grads[3*n+1] += dabc[1]*coefs[n];
718 grads[3*n+2] += dabc[2]*coefs[n];
719 lapl3[3*n+0] += d2abc[0]*coefs[n];
720 lapl3[3*n+1] += d2abc[1]*coefs[n];
721 lapl3[3*n+2] += d2abc[2]*coefs[n];
722 }
723 }
724
725 double dxInv = spline->x_grid.delta_inv;
726 double dyInv = spline->y_grid.delta_inv;
727 double dzInv = spline->z_grid.delta_inv;
728 for (int n=0; n<spline->num_splines; n++) {
729 grads[3*n+0] *= dxInv;
730 grads[3*n+1] *= dyInv;
731 grads[3*n+2] *= dzInv;
732 lapl3[3*n+0] *= dxInv*dxInv;
733 lapl3[3*n+1] *= dyInv*dyInv;
734 lapl3[3*n+2] *= dzInv*dzInv;
735 lapl[n] = lapl3[3*n+0] + lapl3[3*n+1] + lapl3[3*n+2];
736 }
737}
738
739
740
741void
742eval_multi_UBspline_3d_z_vgh (multi_UBspline_3d_z *spline,
743 double x, double y, double z,
744 complex_double* restrict vals,
745 complex_double* restrict grads,
746 complex_double* restrict hess)
747{
748 x -= spline->x_grid.start;
749 y -= spline->y_grid.start;
750 z -= spline->z_grid.start;
751 double ux = x*spline->x_grid.delta_inv;
752 double uy = y*spline->y_grid.delta_inv;
753 double uz = z*spline->z_grid.delta_inv;
754 double ipartx, iparty, ipartz, tx, ty, tz;
755 tx = modf (ux, &ipartx); int ix = (int) ipartx;
756 ty = modf (uy, &iparty); int iy = (int) iparty;
757 tz = modf (uz, &ipartz); int iz = (int) ipartz;
758
759 double tpx[4], tpy[4], tpz[4], a[4], b[4], c[4],
760 da[4], db[4], dc[4], d2a[4], d2b[4], d2c[4];
761 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
762 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
763 tpz[0] = tz*tz*tz; tpz[1] = tz*tz; tpz[2] = tz; tpz[3] = 1.0;
764 complex_double* restrict coefs = spline->coefs;
765
766 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
767 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
768 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
769 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
770 da[0] = (dAd[ 0]*tpx[0] + dAd[ 1]*tpx[1] + dAd[ 2]*tpx[2] + dAd[ 3]*tpx[3]);
771 da[1] = (dAd[ 4]*tpx[0] + dAd[ 5]*tpx[1] + dAd[ 6]*tpx[2] + dAd[ 7]*tpx[3]);
772 da[2] = (dAd[ 8]*tpx[0] + dAd[ 9]*tpx[1] + dAd[10]*tpx[2] + dAd[11]*tpx[3]);
773 da[3] = (dAd[12]*tpx[0] + dAd[13]*tpx[1] + dAd[14]*tpx[2] + dAd[15]*tpx[3]);
774 d2a[0] = (d2Ad[ 0]*tpx[0] + d2Ad[ 1]*tpx[1] + d2Ad[ 2]*tpx[2] + d2Ad[ 3]*tpx[3]);
775 d2a[1] = (d2Ad[ 4]*tpx[0] + d2Ad[ 5]*tpx[1] + d2Ad[ 6]*tpx[2] + d2Ad[ 7]*tpx[3]);
776 d2a[2] = (d2Ad[ 8]*tpx[0] + d2Ad[ 9]*tpx[1] + d2Ad[10]*tpx[2] + d2Ad[11]*tpx[3]);
777 d2a[3] = (d2Ad[12]*tpx[0] + d2Ad[13]*tpx[1] + d2Ad[14]*tpx[2] + d2Ad[15]*tpx[3]);
778
779 b[0] = (Ad[ 0]*tpy[0] + Ad[ 1]*tpy[1] + Ad[ 2]*tpy[2] + Ad[ 3]*tpy[3]);
780 b[1] = (Ad[ 4]*tpy[0] + Ad[ 5]*tpy[1] + Ad[ 6]*tpy[2] + Ad[ 7]*tpy[3]);
781 b[2] = (Ad[ 8]*tpy[0] + Ad[ 9]*tpy[1] + Ad[10]*tpy[2] + Ad[11]*tpy[3]);
782 b[3] = (Ad[12]*tpy[0] + Ad[13]*tpy[1] + Ad[14]*tpy[2] + Ad[15]*tpy[3]);
783 db[0] = (dAd[ 0]*tpy[0] + dAd[ 1]*tpy[1] + dAd[ 2]*tpy[2] + dAd[ 3]*tpy[3]);
784 db[1] = (dAd[ 4]*tpy[0] + dAd[ 5]*tpy[1] + dAd[ 6]*tpy[2] + dAd[ 7]*tpy[3]);
785 db[2] = (dAd[ 8]*tpy[0] + dAd[ 9]*tpy[1] + dAd[10]*tpy[2] + dAd[11]*tpy[3]);
786 db[3] = (dAd[12]*tpy[0] + dAd[13]*tpy[1] + dAd[14]*tpy[2] + dAd[15]*tpy[3]);
787 d2b[0] = (d2Ad[ 0]*tpy[0] + d2Ad[ 1]*tpy[1] + d2Ad[ 2]*tpy[2] + d2Ad[ 3]*tpy[3]);
788 d2b[1] = (d2Ad[ 4]*tpy[0] + d2Ad[ 5]*tpy[1] + d2Ad[ 6]*tpy[2] + d2Ad[ 7]*tpy[3]);
789 d2b[2] = (d2Ad[ 8]*tpy[0] + d2Ad[ 9]*tpy[1] + d2Ad[10]*tpy[2] + d2Ad[11]*tpy[3]);
790 d2b[3] = (d2Ad[12]*tpy[0] + d2Ad[13]*tpy[1] + d2Ad[14]*tpy[2] + d2Ad[15]*tpy[3]);
791
792 c[0] = (Ad[ 0]*tpz[0] + Ad[ 1]*tpz[1] + Ad[ 2]*tpz[2] + Ad[ 3]*tpz[3]);
793 c[1] = (Ad[ 4]*tpz[0] + Ad[ 5]*tpz[1] + Ad[ 6]*tpz[2] + Ad[ 7]*tpz[3]);
794 c[2] = (Ad[ 8]*tpz[0] + Ad[ 9]*tpz[1] + Ad[10]*tpz[2] + Ad[11]*tpz[3]);
795 c[3] = (Ad[12]*tpz[0] + Ad[13]*tpz[1] + Ad[14]*tpz[2] + Ad[15]*tpz[3]);
796 dc[0] = (dAd[ 0]*tpz[0] + dAd[ 1]*tpz[1] + dAd[ 2]*tpz[2] + dAd[ 3]*tpz[3]);
797 dc[1] = (dAd[ 4]*tpz[0] + dAd[ 5]*tpz[1] + dAd[ 6]*tpz[2] + dAd[ 7]*tpz[3]);
798 dc[2] = (dAd[ 8]*tpz[0] + dAd[ 9]*tpz[1] + dAd[10]*tpz[2] + dAd[11]*tpz[3]);
799 dc[3] = (dAd[12]*tpz[0] + dAd[13]*tpz[1] + dAd[14]*tpz[2] + dAd[15]*tpz[3]);
800 d2c[0] = (d2Ad[ 0]*tpz[0] + d2Ad[ 1]*tpz[1] + d2Ad[ 2]*tpz[2] + d2Ad[ 3]*tpz[3]);
801 d2c[1] = (d2Ad[ 4]*tpz[0] + d2Ad[ 5]*tpz[1] + d2Ad[ 6]*tpz[2] + d2Ad[ 7]*tpz[3]);
802 d2c[2] = (d2Ad[ 8]*tpz[0] + d2Ad[ 9]*tpz[1] + d2Ad[10]*tpz[2] + d2Ad[11]*tpz[3]);
803 d2c[3] = (d2Ad[12]*tpz[0] + d2Ad[13]*tpz[1] + d2Ad[14]*tpz[2] + d2Ad[15]*tpz[3]);
804
805 intptr_t xs = spline->x_stride;
806 intptr_t ys = spline->y_stride;
807 intptr_t zs = spline->z_stride;
808
809 for (int n=0; n<spline->num_splines; n++) {
810 vals[n] = 0.0;
811 grads[3*n+0] = grads[3*n+1] = grads[3*n+2] = 0.0;
812 for (int i=0; i<9; i++)
813 hess[9*n+i] = 0.0;
814 }
815
816 for (int i=0; i<4; i++)
817 for (int j=0; j<4; j++)
818 for (int k=0; k<4; k++) {
819 double abc = a[i]*b[j]*c[k];
820 double dabc[3], d2abc[6];
821 dabc[0] = da[i]* b[j]* c[k];
822 dabc[1] = a[i]*db[j]* c[k];
823 dabc[2] = a[i]* b[j]*dc[k];
824 d2abc[0] = d2a[i]* b[j]* c[k];
825 d2abc[1] = da[i]* db[j]* c[k];
826 d2abc[2] = da[i]* b[j]* dc[k];
827 d2abc[3] = a[i]*d2b[j]* c[k];
828 d2abc[4] = a[i]* db[j]* dc[k];
829 d2abc[5] = a[i]* b[j]*d2c[k];
830
831 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys + (iz+k)*zs);
832 for (int n=0; n<spline->num_splines; n++) {
833 vals[n] += abc *coefs[n];
834 grads[3*n+0] += dabc[0]*coefs[n];
835 grads[3*n+1] += dabc[1]*coefs[n];
836 grads[3*n+2] += dabc[2]*coefs[n];
837 hess [9*n+0] += d2abc[0]*coefs[n];
838 hess [9*n+1] += d2abc[1]*coefs[n];
839 hess [9*n+2] += d2abc[2]*coefs[n];
840 hess [9*n+4] += d2abc[3]*coefs[n];
841 hess [9*n+5] += d2abc[4]*coefs[n];
842 hess [9*n+8] += d2abc[5]*coefs[n];
843 }
844 }
845
846 double dxInv = spline->x_grid.delta_inv;
847 double dyInv = spline->y_grid.delta_inv;
848 double dzInv = spline->z_grid.delta_inv;
849 for (int n=0; n<spline->num_splines; n++) {
850 grads[3*n+0] *= dxInv;
851 grads[3*n+1] *= dyInv;
852 grads[3*n+2] *= dzInv;
853 hess[9*n+0] *= dxInv*dxInv;
854 hess[9*n+4] *= dyInv*dyInv;
855 hess[9*n+8] *= dzInv*dzInv;
856 hess[9*n+1] *= dxInv*dyInv;
857 hess[9*n+2] *= dxInv*dzInv;
858 hess[9*n+5] *= dyInv*dzInv;
859 // Copy hessian elements into lower half of 3x3 matrix
860 hess[9*n+3] = hess[9*n+1];
861 hess[9*n+6] = hess[9*n+2];
862 hess[9*n+7] = hess[9*n+5];
863 }
864}
865
866
867void
868eval_multi_UBspline_3d_z_vghgh (multi_UBspline_3d_z *spline,
869 double x, double y, double z,
870 complex_double* restrict vals,
871 complex_double* restrict grads,
872 complex_double* restrict hess,
873 complex_double* restrict gradhess)
874{
875 x -= spline->x_grid.start;
876 y -= spline->y_grid.start;
877 z -= spline->z_grid.start;
878 double ux = x*spline->x_grid.delta_inv;
879 double uy = y*spline->y_grid.delta_inv;
880 double uz = z*spline->z_grid.delta_inv;
881 double ipartx, iparty, ipartz, tx, ty, tz;
882 tx = modf (ux, &ipartx); int ix = (int) ipartx;
883 ty = modf (uy, &iparty); int iy = (int) iparty;
884 tz = modf (uz, &ipartz); int iz = (int) ipartz;
885
886 double tpx[4], tpy[4], tpz[4], a[4], b[4], c[4],
887 da[4], db[4], dc[4], d2a[4], d2b[4], d2c[4],
888 d3a[4], d3b[4], d3c[4];
889 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
890 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
891 tpz[0] = tz*tz*tz; tpz[1] = tz*tz; tpz[2] = tz; tpz[3] = 1.0;
892 complex_double* restrict coefs = spline->coefs;
893
894 a[0] = (Ad[ 0]*tpx[0] + Ad[ 1]*tpx[1] + Ad[ 2]*tpx[2] + Ad[ 3]*tpx[3]);
895 a[1] = (Ad[ 4]*tpx[0] + Ad[ 5]*tpx[1] + Ad[ 6]*tpx[2] + Ad[ 7]*tpx[3]);
896 a[2] = (Ad[ 8]*tpx[0] + Ad[ 9]*tpx[1] + Ad[10]*tpx[2] + Ad[11]*tpx[3]);
897 a[3] = (Ad[12]*tpx[0] + Ad[13]*tpx[1] + Ad[14]*tpx[2] + Ad[15]*tpx[3]);
898 da[0] = (dAd[ 0]*tpx[0] + dAd[ 1]*tpx[1] + dAd[ 2]*tpx[2] + dAd[ 3]*tpx[3]);
899 da[1] = (dAd[ 4]*tpx[0] + dAd[ 5]*tpx[1] + dAd[ 6]*tpx[2] + dAd[ 7]*tpx[3]);
900 da[2] = (dAd[ 8]*tpx[0] + dAd[ 9]*tpx[1] + dAd[10]*tpx[2] + dAd[11]*tpx[3]);
901 da[3] = (dAd[12]*tpx[0] + dAd[13]*tpx[1] + dAd[14]*tpx[2] + dAd[15]*tpx[3]);
902 d2a[0] = (d2Ad[ 0]*tpx[0] + d2Ad[ 1]*tpx[1] + d2Ad[ 2]*tpx[2] + d2Ad[ 3]*tpx[3]);
903 d2a[1] = (d2Ad[ 4]*tpx[0] + d2Ad[ 5]*tpx[1] + d2Ad[ 6]*tpx[2] + d2Ad[ 7]*tpx[3]);
904 d2a[2] = (d2Ad[ 8]*tpx[0] + d2Ad[ 9]*tpx[1] + d2Ad[10]*tpx[2] + d2Ad[11]*tpx[3]);
905 d2a[3] = (d2Ad[12]*tpx[0] + d2Ad[13]*tpx[1] + d2Ad[14]*tpx[2] + d2Ad[15]*tpx[3]);
906 d3a[0] = (/*d2Ad[ 0]*tpx[0] + d2Ad[ 1]*tpx[1] + d2Ad[ 2]*tpx[2] +*/ d3Ad[ 3]*tpx[3]);
907 d3a[1] = (/*d2Ad[ 4]*tpx[0] + d2Ad[ 5]*tpx[1] + d2Ad[ 6]*tpx[2] +*/ d3Ad[ 7]*tpx[3]);
908 d3a[2] = (/*d2Ad[ 8]*tpx[0] + d2Ad[ 9]*tpx[1] + d2Ad[10]*tpx[2] +*/ d3Ad[11]*tpx[3]);
909 d3a[3] = (/*d2Ad[12]*tpx[0] + d2Ad[13]*tpx[1] + d2Ad[14]*tpx[2] +*/ d3Ad[15]*tpx[3]);
910
911 b[0] = (Ad[ 0]*tpy[0] + Ad[ 1]*tpy[1] + Ad[ 2]*tpy[2] + Ad[ 3]*tpy[3]);
912 b[1] = (Ad[ 4]*tpy[0] + Ad[ 5]*tpy[1] + Ad[ 6]*tpy[2] + Ad[ 7]*tpy[3]);
913 b[2] = (Ad[ 8]*tpy[0] + Ad[ 9]*tpy[1] + Ad[10]*tpy[2] + Ad[11]*tpy[3]);
914 b[3] = (Ad[12]*tpy[0] + Ad[13]*tpy[1] + Ad[14]*tpy[2] + Ad[15]*tpy[3]);
915 db[0] = (dAd[ 0]*tpy[0] + dAd[ 1]*tpy[1] + dAd[ 2]*tpy[2] + dAd[ 3]*tpy[3]);
916 db[1] = (dAd[ 4]*tpy[0] + dAd[ 5]*tpy[1] + dAd[ 6]*tpy[2] + dAd[ 7]*tpy[3]);
917 db[2] = (dAd[ 8]*tpy[0] + dAd[ 9]*tpy[1] + dAd[10]*tpy[2] + dAd[11]*tpy[3]);
918 db[3] = (dAd[12]*tpy[0] + dAd[13]*tpy[1] + dAd[14]*tpy[2] + dAd[15]*tpy[3]);
919 d2b[0] = (d2Ad[ 0]*tpy[0] + d2Ad[ 1]*tpy[1] + d2Ad[ 2]*tpy[2] + d2Ad[ 3]*tpy[3]);
920 d2b[1] = (d2Ad[ 4]*tpy[0] + d2Ad[ 5]*tpy[1] + d2Ad[ 6]*tpy[2] + d2Ad[ 7]*tpy[3]);
921 d2b[2] = (d2Ad[ 8]*tpy[0] + d2Ad[ 9]*tpy[1] + d2Ad[10]*tpy[2] + d2Ad[11]*tpy[3]);
922 d2b[3] = (d2Ad[12]*tpy[0] + d2Ad[13]*tpy[1] + d2Ad[14]*tpy[2] + d2Ad[15]*tpy[3]);
923 d3b[0] = (/*d2Ad[ 0]*tpx[0] + d2Ad[ 1]*tpx[1] + d2Ad[ 2]*tpx[2] +*/ d3Ad[ 3]*tpy[3]);
924 d3b[1] = (/*d2Ad[ 4]*tpx[0] + d2Ad[ 5]*tpx[1] + d2Ad[ 6]*tpx[2] +*/ d3Ad[ 7]*tpy[3]);
925 d3b[2] = (/*d2Ad[ 8]*tpx[0] + d2Ad[ 9]*tpx[1] + d2Ad[10]*tpx[2] +*/ d3Ad[11]*tpy[3]);
926 d3b[3] = (/*d2Ad[12]*tpx[0] + d2Ad[13]*tpx[1] + d2Ad[14]*tpx[2] +*/ d3Ad[15]*tpy[3]);
927
928 c[0] = (Ad[ 0]*tpz[0] + Ad[ 1]*tpz[1] + Ad[ 2]*tpz[2] + Ad[ 3]*tpz[3]);
929 c[1] = (Ad[ 4]*tpz[0] + Ad[ 5]*tpz[1] + Ad[ 6]*tpz[2] + Ad[ 7]*tpz[3]);
930 c[2] = (Ad[ 8]*tpz[0] + Ad[ 9]*tpz[1] + Ad[10]*tpz[2] + Ad[11]*tpz[3]);
931 c[3] = (Ad[12]*tpz[0] + Ad[13]*tpz[1] + Ad[14]*tpz[2] + Ad[15]*tpz[3]);
932 dc[0] = (dAd[ 0]*tpz[0] + dAd[ 1]*tpz[1] + dAd[ 2]*tpz[2] + dAd[ 3]*tpz[3]);
933 dc[1] = (dAd[ 4]*tpz[0] + dAd[ 5]*tpz[1] + dAd[ 6]*tpz[2] + dAd[ 7]*tpz[3]);
934 dc[2] = (dAd[ 8]*tpz[0] + dAd[ 9]*tpz[1] + dAd[10]*tpz[2] + dAd[11]*tpz[3]);
935 dc[3] = (dAd[12]*tpz[0] + dAd[13]*tpz[1] + dAd[14]*tpz[2] + dAd[15]*tpz[3]);
936 d2c[0] = (d2Ad[ 0]*tpz[0] + d2Ad[ 1]*tpz[1] + d2Ad[ 2]*tpz[2] + d2Ad[ 3]*tpz[3]);
937 d2c[1] = (d2Ad[ 4]*tpz[0] + d2Ad[ 5]*tpz[1] + d2Ad[ 6]*tpz[2] + d2Ad[ 7]*tpz[3]);
938 d2c[2] = (d2Ad[ 8]*tpz[0] + d2Ad[ 9]*tpz[1] + d2Ad[10]*tpz[2] + d2Ad[11]*tpz[3]);
939 d2c[3] = (d2Ad[12]*tpz[0] + d2Ad[13]*tpz[1] + d2Ad[14]*tpz[2] + d2Ad[15]*tpz[3]);
940 d3c[0] = (/*d2Ad[ 0]*tpx[0] + d2Ad[ 1]*tpx[1] + d2Ad[ 2]*tpx[2] +*/ d3Ad[ 3]*tpz[3]);
941 d3c[1] = (/*d2Ad[ 4]*tpx[0] + d2Ad[ 5]*tpx[1] + d2Ad[ 6]*tpx[2] +*/ d3Ad[ 7]*tpz[3]);
942 d3c[2] = (/*d2Ad[ 8]*tpx[0] + d2Ad[ 9]*tpx[1] + d2Ad[10]*tpx[2] +*/ d3Ad[11]*tpz[3]);
943 d3c[3] = (/*d2Ad[12]*tpx[0] + d2Ad[13]*tpx[1] + d2Ad[14]*tpx[2] +*/ d3Ad[15]*tpz[3]);
944
945 intptr_t xs = spline->x_stride;
946 intptr_t ys = spline->y_stride;
947 intptr_t zs = spline->z_stride;
948
949 for (int n=0; n<spline->num_splines; n++) {
950 vals[n] = 0.0;
951 grads[3*n+0] = grads[3*n+1] = grads[3*n+2] = 0.0;
952 for (int i=0; i<9; i++)
953 hess[9*n+i] = 0.0;
954 for (int i=0; i<27; i++)
955 gradhess[27*n+i] = 0.0;
956 }
957
958 for (int i=0; i<4; i++)
959 for (int j=0; j<4; j++)
960 for (int k=0; k<4; k++) {
961 double abc = a[i]*b[j]*c[k];
962 double dabc[3], d2abc[6], d3abc[10];
963 dabc[0] = da[i]* b[j]* c[k];
964 dabc[1] = a[i]*db[j]* c[k];
965 dabc[2] = a[i]* b[j]*dc[k];
966 d2abc[0] = d2a[i]* b[j]* c[k];
967 d2abc[1] = da[i]* db[j]* c[k];
968 d2abc[2] = da[i]* b[j]* dc[k];
969 d2abc[3] = a[i]*d2b[j]* c[k];
970 d2abc[4] = a[i]* db[j]* dc[k];
971 d2abc[5] = a[i]* b[j]*d2c[k];
972
973 d3abc[0] = d3a[i]* b[j]* c[k];
974 d3abc[1] = d2a[i]* db[j]* c[k];
975 d3abc[2] = d2a[i]* b[j]* dc[k];
976 d3abc[3] = da[i]*d2b[j]* c[k];
977 d3abc[4] = da[i]* db[j]* dc[k];
978 d3abc[5] = da[i]* b[j]*d2c[k];
979 d3abc[6] = a[i]*d3b[j]* c[k];
980 d3abc[7] = a[i]*d2b[j]* dc[k];
981 d3abc[8] = a[i]* db[j]*d2c[k];
982 d3abc[9] = a[i]* b[j]*d3c[k];
983
984 complex_double* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys + (iz+k)*zs);
985 for (int n=0; n<spline->num_splines; n++) {
986 vals[n] += abc *coefs[n];
987 grads[3*n+0] += dabc[0]*coefs[n];
988 grads[3*n+1] += dabc[1]*coefs[n];
989 grads[3*n+2] += dabc[2]*coefs[n];
990 hess [9*n+0] += d2abc[0]*coefs[n];
991 hess [9*n+1] += d2abc[1]*coefs[n];
992 hess [9*n+2] += d2abc[2]*coefs[n];
993 hess [9*n+4] += d2abc[3]*coefs[n];
994 hess [9*n+5] += d2abc[4]*coefs[n];
995 hess [9*n+8] += d2abc[5]*coefs[n];
996
997 gradhess [27*n+0 ] += d3abc[0]*coefs[n];
998 gradhess [27*n+1 ] += d3abc[1]*coefs[n];
999 gradhess [27*n+2 ] += d3abc[2]*coefs[n];
1000 gradhess [27*n+4 ] += d3abc[3]*coefs[n];
1001 gradhess [27*n+5 ] += d3abc[4]*coefs[n];
1002 gradhess [27*n+8 ] += d3abc[5]*coefs[n];
1003 gradhess [27*n+13] += d3abc[6]*coefs[n];
1004 gradhess [27*n+14] += d3abc[7]*coefs[n];
1005 gradhess [27*n+17] += d3abc[8]*coefs[n];
1006 gradhess [27*n+26] += d3abc[9]*coefs[n];
1007 }
1008 }
1009
1010 double dxInv = spline->x_grid.delta_inv;
1011 double dyInv = spline->y_grid.delta_inv;
1012 double dzInv = spline->z_grid.delta_inv;
1013 for (int n=0; n<spline->num_splines; n++) {
1014 grads[3*n+0] *= dxInv;
1015 grads[3*n+1] *= dyInv;
1016 grads[3*n+2] *= dzInv;
1017 hess [9*n+0] *= dxInv*dxInv;
1018 hess [9*n+4] *= dyInv*dyInv;
1019 hess [9*n+8] *= dzInv*dzInv;
1020 hess [9*n+1] *= dxInv*dyInv;
1021 hess [9*n+2] *= dxInv*dzInv;
1022 hess [9*n+5] *= dyInv*dzInv;
1023 // Copy hessian elements into lower half of 3x3 matrix
1024 hess [9*n+3] = hess[9*n+1];
1025 hess [9*n+6] = hess[9*n+2];
1026 hess [9*n+7] = hess[9*n+5];
1027
1028 gradhess [27*n+0 ] *= dxInv*dxInv*dxInv;
1029 gradhess [27*n+1 ] *= dxInv*dxInv*dyInv;
1030 gradhess [27*n+2 ] *= dxInv*dxInv*dzInv;
1031 gradhess [27*n+4 ] *= dxInv*dyInv*dyInv;
1032 gradhess [27*n+5 ] *= dxInv*dyInv*dzInv;
1033 gradhess [27*n+8 ] *= dxInv*dzInv*dzInv;
1034 gradhess [27*n+13] *= dyInv*dyInv*dyInv;
1035 gradhess [27*n+14] *= dyInv*dyInv*dzInv;
1036 gradhess [27*n+17] *= dyInv*dzInv*dzInv;
1037 gradhess [27*n+26] *= dzInv*dzInv*dzInv;
1038
1039 // Copy gradhess elements into rest of tensor
1040 gradhess [27*n+9 ] = gradhess [27*n+3 ] = gradhess [27*n+1 ];
1041 gradhess [27*n+18 ] = gradhess [27*n+6 ] = gradhess [27*n+2 ];
1042 gradhess [27*n+22 ] = gradhess [27*n+16 ] = gradhess [27*n+14];
1043 gradhess [27*n+12 ] = gradhess [27*n+10 ] = gradhess [27*n+4 ];
1044 gradhess [27*n+24 ] = gradhess [27*n+20 ] = gradhess [27*n+8 ];
1045 gradhess [27*n+25 ] = gradhess [27*n+23 ] = gradhess [27*n+17];
1046 gradhess [27*n+21 ] = gradhess [27*n+19 ] = gradhess [27*n+15] = gradhess [27*n+11 ] = gradhess [27*n+7 ] = gradhess [27*n+5];
1047
1048 }
1049}
1050
1051#endif
bspline_base.h
complex_double
complex double complex_double
Definition bspline_base.h:34
restrict
#define restrict
Definition local_definitions.h:4
eval_multi_UBspline_2d_z
void eval_multi_UBspline_2d_z(multi_UBspline_2d_z *spline, double x, double y, complex_double *restrict vals)
Definition multi_bspline_eval_std_z.h:184
eval_multi_UBspline_2d_z_vgh
void eval_multi_UBspline_2d_z_vgh(multi_UBspline_2d_z *spline, double x, double y, complex_double *restrict vals, complex_double *restrict grads, complex_double *restrict hess)
Definition multi_bspline_eval_std_z.h:386
eval_multi_UBspline_2d_z_vg
void eval_multi_UBspline_2d_z_vg(multi_UBspline_2d_z *spline, double x, double y, complex_double *restrict vals, complex_double *restrict grads)
Definition multi_bspline_eval_std_z.h:228
eval_multi_UBspline_1d_z_vgl
void eval_multi_UBspline_1d_z_vgl(multi_UBspline_1d_z *spline, double x, complex_double *restrict vals, complex_double *restrict grads, complex_double *restrict lapl)
Definition multi_bspline_eval_std_z.h:116
eval_multi_UBspline_3d_z_vgh
void eval_multi_UBspline_3d_z_vgh(multi_UBspline_3d_z *spline, double x, double y, double z, complex_double *restrict vals, complex_double *restrict grads, complex_double *restrict hess)
Definition multi_bspline_eval_std_z.h:742
Ad
const double *restrict Ad
Definition bspline_data.cpp:181
d2Ad
const double *restrict d2Ad
Definition bspline_data.cpp:195
eval_multi_UBspline_2d_z_vgl
void eval_multi_UBspline_2d_z_vgl(multi_UBspline_2d_z *spline, double x, double y, complex_double *restrict vals, complex_double *restrict grads, complex_double *restrict lapl)
Definition multi_bspline_eval_std_z.h:296
eval_multi_UBspline_1d_z
void eval_multi_UBspline_1d_z(multi_UBspline_1d_z *spline, double x, complex_double *restrict vals)
Definition multi_bspline_eval_std_z.h:38
eval_multi_UBspline_3d_z_vgl
void eval_multi_UBspline_3d_z_vgl(multi_UBspline_3d_z *spline, double x, double y, double z, complex_double *restrict vals, complex_double *restrict grads, complex_double *restrict lapl)
Definition multi_bspline_eval_std_z.h:626
dAd
const double *restrict dAd
Definition bspline_data.cpp:188
d3Ad
const double *restrict d3Ad
Definition bspline_data.cpp:202
eval_multi_UBspline_3d_z_vg
void eval_multi_UBspline_3d_z_vg(multi_UBspline_3d_z *spline, double x, double y, double z, complex_double *restrict vals, complex_double *restrict grads)
Definition multi_bspline_eval_std_z.h:536
eval_multi_UBspline_3d_z_vghgh
void eval_multi_UBspline_3d_z_vghgh(multi_UBspline_3d_z *spline, double x, double y, double z, complex_double *restrict vals, complex_double *restrict grads, complex_double *restrict hess, complex_double *restrict gradhess)
Definition multi_bspline_eval_std_z.h:868
eval_multi_UBspline_1d_z_vg
void eval_multi_UBspline_1d_z_vg(multi_UBspline_1d_z *spline, double x, complex_double *restrict vals, complex_double *restrict grads)
Definition multi_bspline_eval_std_z.h:71
eval_multi_UBspline_1d_z_vgh
void eval_multi_UBspline_1d_z_vgh(multi_UBspline_1d_z *spline, double x, complex_double *restrict vals, complex_double *restrict grads, complex_double *restrict hess)
Definition multi_bspline_eval_std_z.h:170
eval_multi_UBspline_3d_z
void eval_multi_UBspline_3d_z(multi_UBspline_3d_z *spline, double x, double y, double z, complex_double *restrict vals)
Definition multi_bspline_eval_std_z.h:481
multi_bspline_structs.h
Ugrid::start
double start
Definition bspline_base.h:85
Ugrid::delta_inv
double delta_inv
Definition bspline_base.h:89
multi_UBspline_1d_z
Definition multi_bspline_structs.h:145
multi_UBspline_1d_z::coefs
complex_double *restrict coefs
Definition multi_bspline_structs.h:148
multi_UBspline_1d_z::x_stride
intptr_t x_stride
Definition multi_bspline_structs.h:149
multi_UBspline_1d_z::x_grid
Ugrid x_grid
Definition multi_bspline_structs.h:150
multi_UBspline_1d_z::num_splines
int num_splines
Definition multi_bspline_structs.h:152
multi_UBspline_2d_z
Definition multi_bspline_structs.h:156
multi_UBspline_2d_z::x_stride
intptr_t x_stride
Definition multi_bspline_structs.h:160
multi_UBspline_2d_z::x_grid
Ugrid x_grid
Definition multi_bspline_structs.h:161
multi_UBspline_2d_z::y_stride
intptr_t y_stride
Definition multi_bspline_structs.h:160
multi_UBspline_2d_z::coefs
complex_double *restrict coefs
Definition multi_bspline_structs.h:159
multi_UBspline_2d_z::num_splines
int num_splines
Definition multi_bspline_structs.h:163
multi_UBspline_2d_z::lapl2
complex_double *restrict lapl2
Definition multi_bspline_structs.h:165
multi_UBspline_2d_z::y_grid
Ugrid y_grid
Definition multi_bspline_structs.h:161
multi_UBspline_3d_z
Definition multi_bspline_structs.h:169
multi_UBspline_3d_z::coefs
complex_double *restrict coefs
Definition multi_bspline_structs.h:172
multi_UBspline_3d_z::num_splines
int num_splines
Definition multi_bspline_structs.h:176
multi_UBspline_3d_z::x_grid
Ugrid x_grid
Definition multi_bspline_structs.h:174
multi_UBspline_3d_z::x_stride
intptr_t x_stride
Definition multi_bspline_structs.h:173
multi_UBspline_3d_z::z_stride
intptr_t z_stride
Definition multi_bspline_structs.h:173
multi_UBspline_3d_z::lapl3
complex_double *restrict lapl3
Definition multi_bspline_structs.h:178
multi_UBspline_3d_z::y_grid
Ugrid y_grid
Definition multi_bspline_structs.h:174
multi_UBspline_3d_z::z_grid
Ugrid z_grid
Definition multi_bspline_structs.h:174
multi_UBspline_3d_z::y_stride
intptr_t y_stride
Definition multi_bspline_structs.h:173