21 #include "base_spectral.hpp"
22 #include "base_fftw.hpp"
23 #include "headcpp.hpp"
26 #include <unordered_map>
32 std::unordered_map<int, fftw_precomp_t> fftw_precomp_map_i;
35 fftw_precomp_t& coef_i_1d_fftw(
int n) {
36 auto precomp_it = fftw_precomp_map_i.find(n);
37 if(precomp_it == fftw_precomp_map_i.end())
38 precomp_it = fftw_precomp_map_i.emplace(std::piecewise_construct,
39 std::forward_as_tuple(n),
40 std::forward_as_tuple(n,FFTW_HC2R)).first;
41 return (*precomp_it).second;
44 void coef_i_1d_pasprevu (Array<double>&) {
45 cout <<
"Coef_1d not implemented." << endl ;
49 void coef_i_1d_cheb (Array<double>& tab) {
51 assert (tab.get_ndim()==1) ;
52 int nbr = tab.get_size(0) ;
56 auto & fftw_data = coef_i_1d_fftw(nbr-1);
57 double* tab_auxi = fftw_data.buffer;
58 fftw_plan p = fftw_data.plan;
60 double* cf =
new double[nbr];
65 for (
int i=3 ; i<nbr ; i+=2) {
69 double fmoins0 = -(nbr-1)/2 *c1- somme ;
70 for (
int i=3 ; i<nbr ; i+=2)
71 tab_auxi[nbr-1-i/2] = -0.25*(cf[i]-cf[i-2]) ;
72 tab_auxi[0] = tab(0) ;
73 for (
int i=1 ; i<(nbr-1)/2 ; i++)
74 tab_auxi[i] = 0.5*tab(2*i) ;
75 tab_auxi[(nbr-1)/2] = tab(nbr-1) ;
79 for (
int i=1 ; i<(nbr-1)/2 ; i++) {
80 double fp = 0.5*(tab_auxi[i]+tab_auxi[nbr-1-i]) ;
81 double fm = 0.5*(tab_auxi[i]-tab_auxi[nbr-1-i])/sin(M_PI*i/(nbr-1)) ;
83 tab.set(nbr-i-1) = fp -fm ;
85 tab.set(0) = tab_auxi[0] + fmoins0 ;
86 tab.set(nbr-1) = tab_auxi[0] - fmoins0 ;
87 tab.set((nbr-1)/2) = tab_auxi[(nbr-1)/2] ;
93 void coef_i_1d_cheb_even (Array<double>& tab) {
95 assert (tab.get_ndim()==1) ;
96 int nbr = tab.get_size(0) ;
98 auto & fftw_data = coef_i_1d_fftw(nbr-1);
99 double* tab_auxi = fftw_data.buffer;
100 fftw_plan p = fftw_data.plan;
102 double* cf =
new double[nbr] ;
107 for (
int i=3 ; i<nbr ; i+=2) {
108 cf[i] = tab(i) - c1 ;
111 double fmoins0 = (nbr-1)/2 *c1+somme ;
112 for (
int i=3 ; i<nbr ; i+=2)
113 tab_auxi[nbr-1-i/2] = 0.25*(cf[i]-cf[i-2]) ;
114 tab_auxi[0] = tab(0) ;
115 for (
int i=1 ; i<(nbr-1)/2 ; i++)
116 tab_auxi[i] = 0.5*tab(2*i) ;
117 tab_auxi[(nbr-1)/2] = tab(nbr-1) ;
119 fftw_data.execute() ;
121 for (
int i=1 ; i<(nbr-1)/2 ; i++) {
122 double fp = 0.5*(tab_auxi[i]+tab_auxi[nbr-1-i]) ;
123 double fm = 0.5*(tab_auxi[i]-tab_auxi[nbr-1-i])/sin(M_PI*i/(nbr-1)) ;
124 tab.set(nbr-1-i) = fp+fm ;
125 tab.set(i) = fp -fm ;
127 tab.set(0) = tab_auxi[0] - fmoins0 ;
128 tab.set(nbr-1) = tab_auxi[0] + fmoins0 ;
129 tab.set((nbr-1)/2) = tab_auxi[(nbr-1)/2] ;
134 void coef_i_1d_cheb_odd (Array<double>& tab) {
136 assert (tab.get_ndim()==1) ;
137 int nbr = tab.get_size(0) ;
139 auto & fftw_data = coef_i_1d_fftw(nbr-1);
140 double* tab_auxi = fftw_data.buffer;
141 fftw_plan p = fftw_data.plan;
143 double* ti =
new double [nbr] ;
144 double* cf =
new double [nbr] ;
147 for (
int i=1 ; i<nbr-1 ; i++)
148 ti[i] = 0.5*(tab(i)+tab(i-1)) ;
149 ti[nbr-1] = 0.5*tab(nbr-2) ;
154 for (
int i=3 ; i<nbr ; i+=2) {
158 double fmoins0 = (nbr-1)/2 *c1+ somme ;
159 for (
int i=3 ; i<nbr ; i+=2)
160 tab_auxi[nbr-1-i/2] = 0.25*(cf[i]-cf[i-2]) ;
161 tab_auxi[0] = ti[0] ;
162 for (
int i=1 ; i<(nbr-1)/2 ; i++)
163 tab_auxi[i] = 0.5*ti[2*i] ;
164 tab_auxi[(nbr-1)/2]=ti[nbr-1] ;
166 fftw_data.execute() ;
168 for (
int i=1 ; i<(nbr-1)/2 ; i++) {
169 double fp = 0.5*(tab_auxi[i]+tab_auxi[nbr-1-i]) ;
170 double fm = 0.5*(tab_auxi[i]-tab_auxi[nbr-1-i])/sin(M_PI*i/(nbr-1)) ;
171 tab.set(nbr-1-i) = (fp+fm)/sin(M_PI*(nbr-1-i)/2/(nbr-1)) ;
172 tab.set(i) = (fp -fm)/sin(M_PI*i/2/(nbr-1)) ;
175 tab.set(nbr-1) = tab_auxi[0] + fmoins0 ;
176 tab.set((nbr-1)/2) = tab_auxi[(nbr-1)/2]/sin(M_PI*(nbr-1)/4/(nbr-1)) ;
182 double coloc_leg(
int,
int) ;
183 double summation_1d_leg(
double xx,
const Array<double>&) ;
184 double coloc_leg_parity(
int,
int) ;
185 double summation_1d_leg_even(
double xx,
const Array<double>&) ;
186 double summation_1d_leg_odd(
double xx,
const Array<double>&) ;
188 void coef_i_1d_leg (Array<double>& tab) {
190 assert (tab.get_ndim()==1) ;
191 int nbr = tab.get_size(0) ;
193 Array<double> res(nbr) ;
194 for (
int i=0 ; i<nbr ; i++)
195 res.set(i) = summation_1d_leg(coloc_leg(i,nbr), tab) ;
199 void coef_i_1d_leg_even (Array<double>& tab) {
201 assert (tab.get_ndim()==1) ;
202 int nbr = tab.get_size(0) ;
204 Array<double> res(nbr) ;
205 for (
int i=0 ; i<nbr ; i++)
206 res.set(i) = summation_1d_leg_even(coloc_leg_parity(i,nbr), tab) ;
210 void coef_i_1d_leg_odd (Array<double>& tab) {
212 assert (tab.get_ndim()==1) ;
213 int nbr = tab.get_size(0) ;
215 Array<double> res(nbr) ;
216 for (
int i=0 ; i<nbr ; i++)
217 res.set(i) = summation_1d_leg_odd(coloc_leg_parity(i,nbr), tab) ;
221 void coef_i_1d_cossin (Array<double>& tab) {
223 assert (tab.get_ndim()==1) ;
224 int nbr = tab.get_size(0) ;
228 auto & fftw_data = coef_i_1d_fftw(np);
229 double* cf = fftw_data.buffer;
230 fftw_plan p = fftw_data.plan;
233 for (
int i=1 ; i<np/2 ; i++) {
234 cf[i] = 0.5*tab(2*i) ;
235 cf[np-i] = -0.5*tab(2*i+1) ;
238 fftw_data.execute() ;
240 for (
int i=0 ; i<np ; i++)
245 void coef_i_1d_cos (Array<double>& tab) {
247 assert (tab.get_ndim()==1) ;
248 int nbr = tab.get_size(0) ;
250 auto & fftw_data = coef_i_1d_fftw(nbr-1);
251 double* tab_auxi = fftw_data.buffer;
252 fftw_plan p = fftw_data.plan;
254 double* cf =
new double [nbr] ;
259 for (
int i=3 ; i<nbr ; i+=2) {
260 cf[i] = tab(i) - c1 ;
263 double fmoins0 = (nbr-1)/2 *c1+ somme ;
264 for (
int i=3 ; i<nbr ; i+=2)
265 tab_auxi[nbr-1-i/2] = 0.25*(cf[i]-cf[i-2]) ;
266 tab_auxi[0] = tab(0) ;
267 for (
int i=1 ; i<(nbr-1)/2 ; i++)
268 tab_auxi[i] = 0.5*tab(2*i) ;
269 tab_auxi[(nbr-1)/2] = tab(nbr-1) ;
271 fftw_data.execute() ;
273 for (
int i=1 ; i<(nbr-1)/2 ; i++) {
274 double fp = 0.5*(tab_auxi[i]+tab_auxi[nbr-1-i]) ;
275 double fm = 0.5*(tab_auxi[i]-tab_auxi[nbr-1-i])/sin(M_PI*i/(nbr-1)) ;
277 tab.set(nbr-i-1) = fp -fm ;
279 tab.set(0) = tab_auxi[0] + fmoins0 ;
280 tab.set(nbr-1) = tab_auxi[0] - fmoins0 ;
281 tab.set((nbr-1)/2) = tab_auxi[(nbr-1)/2] ;
286 void coef_i_1d_sin (Array<double>& tab) {
288 assert (tab.get_ndim()==1) ;
289 int nbr = tab.get_size(0) ;
291 auto & fftw_data = coef_i_1d_fftw(nbr-1);
292 double* tab_auxi = fftw_data.buffer;
293 fftw_plan p = fftw_data.plan;
295 for (
int i=2 ; i<nbr-1 ; i+=2)
296 tab_auxi[nbr-1-i/2] = -0.5*tab(i) ;
297 tab_auxi[0] = 0.5*tab(1) ;
298 for (
int i=3 ; i<nbr ; i+=2)
299 tab_auxi[i/2] = 0.25*(tab(i)-tab(i-2)) ;
300 tab_auxi[(nbr-1)/2] = -0.5*tab(nbr-2) ;
302 fftw_data.execute() ;
304 for (
int i=1 ; i<(nbr-1)/2 ; i++) {
305 double fp = 0.5*(tab_auxi[i]+tab_auxi[nbr-1-i])/sin(M_PI*i/(nbr-1)) ;
306 double fm = 0.5*(tab_auxi[i]-tab_auxi[nbr-1-i]) ;
308 tab.set(nbr-i-1) = fp -fm ;
311 tab.set(nbr-1) =-2* tab_auxi[0] ;
312 tab.set((nbr-1)/2) = tab_auxi[(nbr-1)/2] ;
315 void coef_i_1d_cos_even (Array<double>& tab) {
317 assert (tab.get_ndim()==1) ;
318 int nbr = tab.get_size(0) ;
320 double* cf =
new double [nbr] ;
322 auto & fftw_data = coef_i_1d_fftw(nbr-1);
323 double* tab_auxi = fftw_data.buffer;
324 fftw_plan p = fftw_data.plan;
329 for (
int i=3 ; i<nbr ; i+=2) {
330 cf[i] = tab(i) - c1 ;
333 double fmoins0 = (nbr-1)/2 *c1+ somme ;
334 for (
int i=3 ; i<nbr ; i+=2)
335 tab_auxi[nbr-1-i/2] = 0.25*(cf[i]-cf[i-2]) ;
336 tab_auxi[0] = tab(0) ;
337 for (
int i=1 ; i<(nbr-1)/2 ; i++)
338 tab_auxi[i] = 0.5*tab(2*i) ;
339 tab_auxi[(nbr-1)/2] = tab(nbr-1) ;
341 fftw_data.execute() ;
343 for (
int i=1 ; i<(nbr-1)/2 ; i++) {
344 double fp = 0.5*(tab_auxi[i]+tab_auxi[nbr-1-i]) ;
345 double fm = 0.5*(tab_auxi[i]-tab_auxi[nbr-1-i])/sin(M_PI*i/(nbr-1)) ;
347 tab.set(nbr-i-1) = fp -fm ;
349 tab.set(0) = tab_auxi[0] + fmoins0 ;
350 tab.set(nbr-1) = tab_auxi[0] - fmoins0 ;
351 tab.set((nbr-1)/2) = tab_auxi[(nbr-1)/2] ;
357 void coef_i_1d_cos_odd (Array<double>& tab) {
359 assert (tab.get_ndim()==1) ;
360 int nbr = tab.get_size(0) ;
362 double* ti =
new double[nbr] ;
363 double* cf =
new double[nbr] ;
365 auto & fftw_data = coef_i_1d_fftw(nbr-1);
366 double* tab_auxi = fftw_data.buffer;
367 fftw_plan p = fftw_data.plan;
370 for (
int i=1 ; i<nbr-1 ; i++)
371 ti[i] = 0.5*(tab(i)+tab(i-1)) ;
372 ti[nbr-1] = 0.5*tab(nbr-2) ;
377 for (
int i=3 ; i<nbr ; i+=2) {
381 double fmoins0 = (nbr-1)/2 *c1+ somme ;
382 for (
int i=3 ; i<nbr ; i+=2)
383 tab_auxi[nbr-1-i/2] = 0.25*(cf[i]-cf[i-2]) ;
384 tab_auxi[0] = ti[0] ;
385 for (
int i=1 ; i<(nbr-1)/2 ; i++)
386 tab_auxi[i] = 0.5*ti[2*i] ;
387 tab_auxi[(nbr-1)/2] = ti[nbr-1] ;
389 fftw_data.execute() ;
391 for (
int i=1 ; i<(nbr-1)/2 ; i++) {
392 double fp = 0.5*(tab_auxi[i]+tab_auxi[nbr-1-i]) ;
393 double fm = 0.5*(tab_auxi[i]-tab_auxi[nbr-1-i])/sin(M_PI*i/(nbr-1)) ;
394 tab.set(i) = (fp+fm)/sin(M_PI*(nbr-1-i)/2/(nbr-1)) ;
395 tab.set(nbr-i-1) = (fp -fm)/sin(M_PI*i/2/(nbr-1)) ;
397 tab.set(0) = tab_auxi[0] + fmoins0 ;
399 tab.set((nbr-1)/2) = tab_auxi[(nbr-1)/2]/sin(M_PI*(nbr-1)/4/(nbr-1)) ;
406 void coef_i_1d_sin_even (Array<double>& tab) {
408 assert (tab.get_ndim()==1) ;
409 int nbr = tab.get_size(0) ;
412 auto & fftw_data = coef_i_1d_fftw(nbr-1);
413 double* tab_auxi = fftw_data.buffer;
414 fftw_plan p = fftw_data.plan;
416 for (
int i=2 ; i<nbr-1 ; i+=2)
417 tab_auxi[nbr-1-i/2] = -0.5*tab(i) ;
418 tab_auxi[0] = 0.5*tab(1) ;
419 for (
int i=1 ; i<(nbr-1)/2 ; i++)
420 tab_auxi[i] = 0.25*(tab(2*i+1)-tab(2*i-1)) ;
421 tab_auxi[(nbr-1)/2] = -0.5*tab(nbr-2) ;
423 fftw_data.execute() ;
425 for (
int i=1 ; i<(nbr-1)/2 ; i++) {
426 double fp = 0.5*(tab_auxi[i]+tab_auxi[nbr-1-i])/sin(M_PI*i/(nbr-1)) ;
427 double fm = 0.5*(tab_auxi[i]-tab_auxi[nbr-1-i]) ;
429 tab.set(nbr-i-1) = fp -fm ;
432 tab.set(nbr-1) =-2* tab_auxi[0] ;
433 tab.set((nbr-1)/2) = tab_auxi[(nbr-1)/2] ;
437 void coef_i_1d_sin_odd (Array<double>& tab) {
439 assert (tab.get_ndim()==1) ;
440 int nbr = tab.get_size(0) ;
442 double* ti =
new double[nbr] ;
443 double* cf =
new double[nbr] ;
445 auto & fftw_data = coef_i_1d_fftw(nbr-1);
446 double* tab_auxi = fftw_data.buffer;
447 fftw_plan p = fftw_data.plan;
450 for (
int i=1 ; i<nbr-1 ; i++)
451 ti[i] = 0.5*(tab(i)-tab(i-1)) ;
452 ti[nbr-1] = -0.5*tab(nbr-2) ;
457 for (
int i=3 ; i<nbr ; i+=2) {
461 double fmoins0 = (nbr-1)/2 *c1+ somme ;
462 for (
int i=3 ; i<nbr ; i+=2)
463 tab_auxi[nbr-1-i/2] = 0.25*(cf[i]-cf[i-2]) ;
464 tab_auxi[0] = ti[0] ;
465 for (
int i=1 ; i<(nbr-1)/2 ; i++)
466 tab_auxi[i] = 0.5*ti[2*i] ;
467 tab_auxi[(nbr-1)/2] = ti[nbr-1] ;
469 fftw_data.execute() ;
471 for (
int i=1 ; i<(nbr-1)/2 ; i++) {
472 double fp = 0.5*(tab_auxi[i]+tab_auxi[nbr-1-i]) ;
473 double fm = 0.5*(tab_auxi[i]-tab_auxi[nbr-1-i])/sin(M_PI*i/(nbr-1)) ;
474 tab.set(i) = (fp+fm)/sin(M_PI*i/2/(nbr-1)) ;
475 tab.set(nbr-i-1) = (fp -fm)/sin(M_PI*(nbr-1-i)/2/(nbr-1)) ;
478 tab.set(nbr-1) = tab_auxi[0] - fmoins0 ;
479 tab.set((nbr-1)/2) = tab_auxi[(nbr-1)/2]/sin(M_PI*(nbr-1)/4/(nbr-1)) ;
487 void coef_i_1d_cossin_even (Array<double>& tab) {
489 assert (tab.get_ndim()==1) ;
490 int nbr = tab.get_size(0) ;
492 Array<double> tab2 (nbr*2-2) ;
494 for (
int i=0 ; i<nbr-2 ; i+=2) {
495 tab2.set(conte) = tab(i) ;
496 tab2.set(conte+1) = tab(i+1) ;
497 tab2.set(conte+2) = 0. ;
498 tab2.set(conte+3) = 0. ;
501 tab2.set(conte) = tab(nbr-2) ;
502 tab2.set(conte+1) = tab(nbr-1) ;
504 coef_i_1d_cossin(tab2) ;
505 for (
int i=0 ; i<nbr-2 ; i++)
506 tab.set(i) = tab2(i) ;
512 void coef_i_1d_cossin_odd (Array<double>& tab) {
514 assert (tab.get_ndim()==1) ;
515 int nbr = tab.get_size(0) ;
517 Array<double> tab2 (nbr*2-2) ;
519 for (
int i=0 ; i<nbr-2 ; i+=2) {
520 tab2.set(conte) = 0. ;
521 tab2.set(conte+1) = 0. ;
522 tab2.set(conte+2) = tab(i) ;
523 tab2.set(conte+3) = tab(i+1) ;
526 tab2.set(conte) = tab(nbr-2) ;
527 tab2.set(conte+1) = tab(nbr-1) ;
529 coef_i_1d_cossin(tab2) ;
530 for (
int i=0 ; i<nbr-2 ; i++)
531 tab.set(i) = tab2(i) ;
537 void coef_i_1d (
int base, Array<double>& tab) {
538 static void (*coef_i_1d[NBR_MAX_BASE])(Array<double>&) ;
539 static bool premier_appel = true ;
543 premier_appel = false ;
545 for (
int i=0; i<NBR_MAX_BASE; i++)
546 coef_i_1d[i] = coef_i_1d_pasprevu ;
548 coef_i_1d[CHEB] = coef_i_1d_cheb ;
549 coef_i_1d[CHEB_EVEN] = coef_i_1d_cheb_even ;
550 coef_i_1d[CHEB_ODD] = coef_i_1d_cheb_odd ;
551 coef_i_1d[COSSIN] = coef_i_1d_cossin ;
552 coef_i_1d[COS_EVEN] = coef_i_1d_cos_even ;
553 coef_i_1d[COS_ODD] = coef_i_1d_cos_odd ;
554 coef_i_1d[SIN_EVEN] = coef_i_1d_sin_even ;
555 coef_i_1d[SIN_ODD] = coef_i_1d_sin_odd ;
556 coef_i_1d[COS] = coef_i_1d_cos ;
557 coef_i_1d[SIN] = coef_i_1d_sin ;
558 coef_i_1d[LEG] = coef_i_1d_leg ;
559 coef_i_1d[LEG_EVEN] = coef_i_1d_leg_even ;
560 coef_i_1d[LEG_ODD] = coef_i_1d_leg_odd;
561 coef_i_1d[COSSIN_EVEN] = coef_i_1d_cossin_even ;
562 coef_i_1d[COSSIN_ODD] = coef_i_1d_cossin_odd ;
565 coef_i_1d[base](tab) ;