21 #include "base_spectral.hpp"
22 #include "base_fftw.hpp"
23 #include "headcpp.hpp"
27 #include <unordered_map>
33 std::unordered_map<int, fftw_precomp_t> fftw_precomp_map;
36 fftw_precomp_t& coef_1d_fftw(
int n) {
37 auto precomp_it = fftw_precomp_map.find(n);
38 if(precomp_it == fftw_precomp_map.end())
39 precomp_it = fftw_precomp_map.emplace(std::piecewise_construct,
40 std::forward_as_tuple(n),
41 std::forward_as_tuple(n,FFTW_R2HC)).first;
42 return (*precomp_it).second;
45 void coef_i_1d (
int, Array<double>&) ;
47 void coef_1d_pasprevu (Array<double>&) {
48 cout <<
"Coef_1d not implemented." << endl ;
52 void coef_1d_cheb (Array<double>& tab) {
54 assert (tab.get_ndim()==1) ;
55 int nr = tab.get_size(0) ;
57 double fmoins0 = 0.5 * ( tab(0) - tab(nr-1));
59 auto & fftw_data = coef_1d_fftw(nr-1);
60 double* tab_auxi = fftw_data.buffer;
61 fftw_plan p = fftw_data.plan;
63 for (
int i=1 ; i<(nr-1)/2 ; i++) {
64 double fp = 0.5*(tab(i)+tab(nr-1-i)) ;
65 double fms = 0.5*(tab(i)-tab(nr-1-i))*sin(M_PI*i/(nr-1)) ;
66 tab_auxi[i] = fp + fms ;
67 tab_auxi[nr-1-i] = fp-fms ;
69 tab_auxi[0] = 0.5*(tab(0) + tab(nr-1)) ;
70 tab_auxi[(nr-1)/2] = tab((nr-1)/2) ;
75 tab.set(0) = tab_auxi[0] / (nr-1) ;
76 for (
int i=2; i<nr-1; i += 2)
77 tab.set(i) = 2*tab_auxi[(i/2)]/(nr-1) ;
78 tab.set(nr-1) = tab_auxi[(nr-1)/2] / (nr-1) ;
83 for (
int i = 3; i < nr; i += 2 ) {
84 tab.set(i) = tab(i-2) -4 * tab_auxi[nr-1-(i-1)/2]/(nr-1) ;
89 double c1 = - ( fmoins0 + som ) / ((nr-1)/2) ;
93 for (
int i = 3; i < nr; i += 2 )
97 void coef_1d_cheb_even (Array<double>& tab) {
98 assert (tab.get_ndim()==1) ;
99 int nr = tab.get_size(0) ;
101 double fmoins0 = - 0.5 * ( tab(0) - tab(nr-1));
103 auto & fftw_data = coef_1d_fftw(nr-1);
104 double* tab_auxi = fftw_data.buffer;
105 fftw_plan p = fftw_data.plan;
107 for (
int i=1 ; i<(nr-1)/2 ; i++) {
108 double fp = 0.5*(tab(i)+tab(nr-1-i)) ;
109 double fms = 0.5*(-tab(i)+tab(nr-1-i))*sin(M_PI*i/(nr-1)) ;
110 tab_auxi[i] = fp + fms ;
111 tab_auxi[nr-1-i] = fp-fms ;
113 tab_auxi[0] = 0.5*(tab(0) + tab(nr-1)) ;
114 tab_auxi[(nr-1)/2] = tab((nr-1)/2) ;
116 fftw_data.execute() ;
119 tab.set(0) = tab_auxi[0] / (nr-1) ;
120 for (
int i=2; i<nr-1; i += 2)
121 tab.set(i) = 2*tab_auxi[(i/2)]/(nr-1) ;
122 tab.set(nr-1) = tab_auxi[(nr-1)/2] / (nr-1) ;
127 for (
int i = 3; i < nr; i += 2 ) {
128 tab.set(i) = tab(i-2) +4 * tab_auxi[nr-1-i/2]/(nr-1) ;
133 double c1 = ( fmoins0 - som ) / ((nr-1)/2) ;
137 for (
int i = 3; i < nr; i += 2 )
141 void coef_1d_cheb_odd (Array<double>& tab) {
142 assert (tab.get_ndim()==1) ;
143 int nr = tab.get_size(0) ;
145 double* cf =
new double[nr];
146 for (
int i=0 ; i<nr ; i++)
147 cf[i] = tab(i) * sin(M_PI/2.*i/(nr-1)) ;
149 double fmoins0 = - 0.5 * ( cf[0] - cf[nr-1] );
151 auto & fftw_data = coef_1d_fftw(nr-1);
152 double* tab_auxi = fftw_data.buffer;
153 fftw_plan p = fftw_data.plan;
155 for (
int i=1 ; i<(nr-1)/2 ; i++) {
156 double fp = 0.5*(cf[i]+cf[nr-1-i]) ;
157 double fms = 0.5*(-cf[i]+cf[nr-1-i])*sin(M_PI*i/(nr-1)) ;
158 tab_auxi[i] = fp + fms ;
159 tab_auxi[nr-1-i] = fp-fms ;
161 tab_auxi[0] = 0.5*(cf[0] + cf[nr-1]) ;
162 tab_auxi[(nr-1)/2] = cf[(nr-1)/2] ;
164 fftw_data.execute() ;
167 cf[0] = tab_auxi[0] / (nr-1) ;
168 for (
int i=2; i<nr-1; i += 2)
169 cf[i] = 2*tab_auxi[(i/2)]/(nr-1) ;
170 cf[nr-1] = tab_auxi[(nr-1)/2] / (nr-1) ;
175 for (
int i = 3; i < nr; i += 2 ) {
176 cf[i] = cf[i-2] +4 * tab_auxi[nr-1-i/2]/(nr-1) ;
181 double c1 = ( fmoins0 - som ) / ((nr-1)/2) ;
185 for (
int i = 3; i < nr; i += 2 )
189 for (
int i=1; i<nr-1; i++)
190 cf[i] = 2 * cf[i] - cf[i-1] ;
193 for (
int i=0 ; i<nr ; i++)
199 double coloc_leg(
int,
int) ;
200 double weight_leg(
int,
int) ;
201 double gamma_leg(
int,
int) ;
202 double leg(
int,
double) ;
203 double coloc_leg_parity(
int,
int) ;
204 double weight_leg_parity(
int,
int) ;
205 double gamma_leg_even(
int,
int) ;
206 double gamma_leg_odd(
int,
int) ;
208 void coef_1d_leg (Array<double>& tab) {
210 assert (tab.get_ndim()==1) ;
211 int nbr = tab.get_size(0) ;
213 Array<double> res (nbr) ;
215 for (
int i=0 ; i<nbr ; i++)
216 for (
int j=0 ; j<nbr ; j++)
217 res.set(i) += tab(j)*leg(i, coloc_leg(j,nbr))*weight_leg(j,nbr)/gamma_leg(i,nbr) ;
221 void coef_1d_leg_even (Array<double>& tab) {
223 assert (tab.get_ndim()==1) ;
224 int nbr = tab.get_size(0) ;
226 Array<double> res (nbr) ;
228 for (
int i=0 ; i<nbr ; i++)
229 for (
int j=0 ; j<nbr ; j++)
230 res.set(i) += tab(j)*leg(i*2, coloc_leg_parity(j,nbr))
231 *weight_leg_parity(j,nbr)/gamma_leg_even(i,nbr) ;
235 void coef_1d_leg_odd (Array<double>& tab) {
237 assert (tab.get_ndim()==1) ;
238 int nbr = tab.get_size(0) ;
240 Array<double> res (nbr) ;
242 for (
int i=0 ; i<nbr-1 ; i++)
243 for (
int j=0 ; j<nbr ; j++)
244 res.set(i) += tab(j)*leg(i*2+1, coloc_leg_parity(j,nbr))
245 *weight_leg_parity(j,nbr)/gamma_leg_odd(i,nbr) ;
249 void coef_1d_cossin (Array<double>& tab) {
251 assert (tab.get_ndim()==1) ;
252 int nbr = tab.get_size(0) ;
255 auto & fftw_data = coef_1d_fftw(nbr-2);
256 double* cf = fftw_data.buffer;
257 fftw_plan p = fftw_data.plan;
259 for (
int i=0 ; i<nbr-2 ; i++)
261 fftw_data.execute() ;
265 double* psin = cf+nbr-3 ;
267 tab.set(index) = *pcos/double(nbr-2) ;
272 for (
int i=1 ; i<nbr/2 - 1; i++) {
273 tab.set(index) = 2.*(*pcos)/double(nbr-2) ;
275 tab.set(index) = -2.*(*psin)/double(nbr-2) ;
279 tab.set(index) = *pcos/double(nbr-2) ;
289 void coef_1d_cos (Array<double>& tab) {
291 assert (tab.get_ndim()==1) ;
292 int nbr = tab.get_size(0) ;
294 double fmoins0 = 0.5 * ( tab(0) - tab(nbr-1) );
296 auto & fftw_data = coef_1d_fftw(nbr-1);
297 double* tab_auxi = fftw_data.buffer;
298 fftw_plan p = fftw_data.plan;
300 for (
int i = 1; i < (nbr-1)/2 ; i++ ) {
301 double fp = 0.5 * ( tab(i) + tab(nbr-1-i) ) ;
302 double fms = 0.5 * ( tab(i) - tab(nbr-1-i) ) * sin(M_PI*i/(nbr-1)) ;
303 tab_auxi[i] = fp + fms ;
304 tab_auxi[nbr-1-i] = fp - fms ;
307 tab_auxi[0] = 0.5 * ( tab(0) + tab(nbr-1) );
308 tab_auxi[(nbr-1)/2] = tab((nbr-1)/2);
310 fftw_data.execute() ;
312 tab.set(0) = tab_auxi[0] / (nbr-1) ;
313 for (
int i=2; i<nbr-1; i += 2 )
314 tab.set(i) = 2*tab_auxi[i/2]/(nbr-1) ;
315 tab.set(nbr-1) = tab_auxi[(nbr-1)/2] / (nbr-1) ;
319 for (
int i = 3; i < nbr; i += 2 ) {
320 tab.set(i) = tab(i-2) + 4 * tab_auxi[nbr-1 - i/2]/(nbr-1) ;
325 double c1 = ( fmoins0 - som ) / ((nbr-1)/2) ;
329 for (
int i = 3; i < nbr; i += 2 )
333 void coef_1d_sin (Array<double>& tab) {
335 assert (tab.get_ndim()==1) ;
336 int nbr = tab.get_size(0) ;
338 auto & fftw_data = coef_1d_fftw(nbr-1);
339 double* tab_auxi = fftw_data.buffer;
340 fftw_plan p = fftw_data.plan;
342 for (
int i = 1; i < (nbr-1)/2 ; i++ ) {
343 double fp = 0.5 * ( tab(i) + tab(nbr-1-i) ) * sin(M_PI*i/(nbr-1)) ;
344 double fms = 0.5 * ( tab(i) - tab(nbr-1-i) ) ;
345 tab_auxi[i] = fp + fms ;
346 tab_auxi[nbr-1-i] = fp - fms ;
349 tab_auxi[0] = 0.5 * ( tab(0) + tab(nbr-1) );
350 tab_auxi[(nbr-1)/2] = tab((nbr-1)/2);
352 fftw_data.execute() ;
355 for (
int i=2; i<nbr-1; i += 2 )
356 tab.set(i) = -2*tab_auxi[nbr-1-i/2]/(nbr-1) ;
359 tab.set(1) = 2*tab_auxi[0]/(nbr-1) ;
360 for (
int i = 3; i < nbr; i += 2 )
361 tab.set(i) = tab(i-2) + 4 * tab_auxi[i/2]/(nbr-1) ;
364 void coef_1d_cos_even (Array<double>& tab) {
365 assert (tab.get_ndim()==1) ;
366 int nbr = tab.get_size(0) ;
369 double fmoins0 = 0.5 * ( tab(0) - tab(nbr-1) );
371 auto & fftw_data = coef_1d_fftw(nbr-1);
372 double* tab_auxi = fftw_data.buffer;
373 fftw_plan p = fftw_data.plan;
375 for (
int i = 1; i < (nbr-1)/2 ; i++ ) {
376 double fp = 0.5 * ( tab(i) + tab(nbr-1-i) ) ;
377 double fms = 0.5 * ( tab(i) - tab(nbr-1-i) ) * sin(M_PI*i/(nbr-1)) ;
378 tab_auxi[i] = fp + fms ;
379 tab_auxi[nbr-1-i] = fp - fms ;
382 tab_auxi[0] = 0.5 * ( tab(0) + tab(nbr-1) );
383 tab_auxi[(nbr-1)/2] = tab((nbr-1)/2);
385 fftw_data.execute() ;
387 tab.set(0) = tab_auxi[0] / (nbr-1) ;
388 for (
int i=2; i<nbr-1; i += 2 )
389 tab.set(i) = 2*tab_auxi[i/2]/(nbr-1) ;
390 tab.set(nbr-1) = tab_auxi[(nbr-1)/2] / (nbr-1) ;
394 for (
int i = 3; i < nbr; i += 2 ) {
395 tab.set(i) = tab(i-2) + 4 * tab_auxi[nbr-1 - i/2]/(nbr-1) ;
400 double c1 = ( fmoins0 - som ) / ((nbr-1)/2) ;
404 for (
int i = 3; i < nbr; i += 2 )
409 void coef_1d_cos_odd (Array<double>& tab) {
410 assert (tab.get_ndim()==1) ;
411 int nbr = tab.get_size(0) ;
414 double* cf =
new double[nbr] ;
415 for (
int i=0 ; i<nbr-1 ; i++)
416 cf[i] = tab(i)*sin(M_PI*(nbr-1-i)/2./(nbr-1)) ;
418 double fmoins0 = 0.5 * (cf[0] - cf[nbr-1] );
420 auto & fftw_data = coef_1d_fftw(nbr-1);
421 double* tab_auxi = fftw_data.buffer;
422 fftw_plan p = fftw_data.plan;
424 for (
int i = 1; i < (nbr-1)/2 ; i++ ) {
425 double fp = 0.5 * ( cf[i] + cf[nbr-1-i] ) ;
426 double fms = 0.5 * ( cf[i] - cf[nbr-1-i] ) * sin(M_PI*i/(nbr-1)) ;
427 tab_auxi[i] = fp + fms ;
428 tab_auxi[nbr-1-i] = fp - fms ;
431 tab_auxi[0] = 0.5 * ( cf[0] + cf[nbr-1] );
432 tab_auxi[(nbr-1)/2] = cf[(nbr-1)/2];
434 fftw_data.execute() ;
436 cf[0] = tab_auxi[0] / (nbr-1) ;
437 for (
int i=2; i<nbr-1; i += 2 )
438 cf[i] = 2*tab_auxi[i/2]/(nbr-1) ;
439 cf[nbr-1] = tab_auxi[(nbr-1)/2] ;
443 for (
int i = 3; i < nbr; i += 2 ) {
444 cf[i] = cf[i-2] + 4 * tab_auxi[nbr-1 - i/2]/(nbr-1) ;
449 double c1 = ( fmoins0 - som ) / ((nbr-1)/2) ;
453 for (
int i = 3; i < nbr; i += 2 )
457 for (
int i=1; i<nbr-1; i++)
458 cf[i] = 2 * cf[i] - cf[i-1] ;
460 for (
int i=0 ; i<nbr ; i++)
467 void coef_1d_sin_even (Array<double>& tab) {
468 assert (tab.get_ndim()==1) ;
469 int nbr = tab.get_size(0) ;
472 auto & fftw_data = coef_1d_fftw(nbr-1);
473 double* tab_auxi = fftw_data.buffer;
474 fftw_plan p = fftw_data.plan;
476 for (
int i = 1; i < (nbr-1)/2 ; i++ ) {
477 double fp = 0.5 * ( tab(i) + tab(nbr-1-i) ) * sin(M_PI*i/(nbr-1));
478 double fms = 0.5 * ( tab(i) - tab(nbr-1-i) ) ;
479 tab_auxi[i] = fp + fms ;
480 tab_auxi[nbr-1-i] = fp - fms ;
483 tab_auxi[0] = 0.5 * ( tab(0) -tab(nbr-1) );
484 tab_auxi[(nbr-1)/2] = tab((nbr-1)/2);
486 fftw_data.execute() ;
489 for (
int i=2; i<nbr-1; i += 2 )
490 tab.set(i) = -2*tab_auxi[nbr-1-i/2]/(nbr-1) ;
493 tab.set(1) = 2*tab_auxi[0]/(nbr-1) ;
494 for (
int i = 3; i < nbr; i += 2 )
495 tab.set(i) = tab(i-2) + 4 * tab_auxi[i/2]/(nbr-1) ;
499 void coef_1d_sin_odd (Array<double>& tab) {
500 assert (tab.get_ndim()==1) ;
501 int nbr = tab.get_size(0) ;
504 double* cf =
new double[nbr] ;
506 for (
int i=1 ; i<nbr ; i++)
507 cf[i] = tab(i)*sin(M_PI*i/2./(nbr-1)) ;
508 double fmoins0 = 0.5 * (cf[0] - cf[nbr-1] );
510 auto & fftw_data = coef_1d_fftw(nbr-1);
511 double* tab_auxi = fftw_data.buffer;
512 fftw_plan p = fftw_data.plan;
514 for (
int i = 1; i < (nbr-1)/2 ; i++ ) {
515 double fp = 0.5 * ( cf[i] + cf[nbr-1-i] ) ;
516 double fms = 0.5 * ( cf[i] - cf[nbr-1-i] ) * sin(M_PI*i/(nbr-1)) ;
517 tab_auxi[i] = fp + fms ;
518 tab_auxi[nbr-1-i] = fp - fms ;
521 tab_auxi[0] = 0.5 * ( cf[0] + cf[nbr-1] );
522 tab_auxi[(nbr-1)/2] = cf[(nbr-1)/2];
524 fftw_data.execute() ;
526 cf[0] = tab_auxi[0] / (nbr-1) ;
527 for (
int i=2; i<nbr-1; i += 2 )
528 cf[i] = 2*tab_auxi[i/2]/(nbr-1) ;
529 cf[nbr-1] = tab_auxi[(nbr-1)/2] / (nbr-1) ;
533 for (
int i = 3; i < nbr; i += 2 ) {
534 cf[i] = cf[i-2] + 4 * tab_auxi[nbr-1 - i/2]/(nbr-1) ;
539 double c1 = ( fmoins0 - som ) / ((nbr-1)/2) ;
543 for (
int i = 3; i < nbr; i += 2 )
547 for (
int i=1; i<nbr-1; i++)
548 cf[i] = 2 * cf[i] + cf[i-1] ;
550 for (
int i=0 ; i<nbr ; i++)
558 void coef_1d_cossin_even (Array<double>& tab) {
561 assert (tab.get_ndim()==1) ;
562 int nbr = tab.get_size(0) ;
563 Array<double> tab2(2*nbr-2) ;
564 for (
int i=0 ; i<nbr-2 ; i++) {
565 tab2.set(i) = tab(i) ;
566 tab2.set(i+nbr-2) = tab(i) ;
569 coef_1d_cossin (tab2) ;
572 for (
int k=0 ; k<nbr-1 ; k+=2) {
573 tab.set(k) = tab2(conte) ;
574 tab.set(k+1) = tab2(conte+1) ;
580 void coef_1d_cossin_odd (Array<double>& tab) {
582 assert (tab.get_ndim()==1) ;
583 int nbr = tab.get_size(0) ;
584 Array<double> tab2(2*nbr-2) ;
586 for (
int i=0 ; i<nbr-2 ; i++) {
587 tab2.set(i) = tab(i) ;
588 tab2.set(i+nbr-2) = -tab(i) ;
590 coef_1d_cossin (tab2) ;
593 for (
int k=0 ; k<nbr-3 ; k+=2) {
594 tab.set(k) = tab2(conte) ;
595 tab.set(k+1) = tab2(conte+1) ;
598 tab.set(nbr-2) = 0. ;
599 tab.set(nbr-1) = 0. ;
602 void coef_1d (
int base, Array<double>& tab) {
603 static void (*coef_1d[NBR_MAX_BASE])(Array<double>&) ;
604 static bool premier_appel = true ;
608 premier_appel = false ;
610 for (
int i=0; i<NBR_MAX_BASE; i++)
611 coef_1d[i] = coef_1d_pasprevu ;
613 coef_1d[CHEB] = coef_1d_cheb ;
614 coef_1d[CHEB_EVEN] = coef_1d_cheb_even ;
615 coef_1d[CHEB_ODD] = coef_1d_cheb_odd ;
616 coef_1d[COSSIN] = coef_1d_cossin ;
617 coef_1d[COS_EVEN] = coef_1d_cos_even ;
618 coef_1d[COS_ODD] = coef_1d_cos_odd ;
619 coef_1d[SIN_EVEN] = coef_1d_sin_even ;
620 coef_1d[SIN_ODD] = coef_1d_sin_odd ;
621 coef_1d[COS] = coef_1d_cos ;
622 coef_1d[SIN] = coef_1d_sin ;
623 coef_1d[LEG] = coef_1d_leg ;
624 coef_1d[LEG_EVEN] = coef_1d_leg_even ;
625 coef_1d[LEG_ODD] = coef_1d_leg_odd;
626 coef_1d[COSSIN_EVEN] = coef_1d_cossin_even ;
627 coef_1d[COSSIN_ODD] = coef_1d_cossin_odd ;