OPAL (Object Oriented Parallel Accelerator Library) 2024.2
OPAL
fftpack.cpp
Go to the documentation of this file.
1//
2// IPPL FFT
3//
4// Copyright (c) 2008-2018
5// Paul Scherrer Institut, Villigen PSI, Switzerland
6// All rights reserved.
7//
8// OPAL is licensed under GNU GPL version 3.
9//
10
11/*
12 * This file is part of libfftpack.
13 *
14 * libfftpack is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License as published by
16 * the Free Software Foundation; either version 2 of the License, or
17 * (at your option) any later version.
18 *
19 * libfftpack is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with libfftpack; if not, write to the Free Software
26 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
27 */
28
29/*
30 * libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
31 * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
32 * (DLR).
33 */
34
35/*
36 fftpack.c : A set of FFT routines in C.
37 Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber
38 (Version 4, 1985).
39
40 C port by Martin Reinecke (2010)
41*/
42#include <cmath>
43#include <cstdio>
44#include <cstdlib>
45#include <cstring>
46
47#ifdef __cplusplus
48extern "C" {
49#endif
50
51/******************************************************************************
52 some helper functions and macros formerly in c_utils.c and c_utils.h
53*/
54 void *util_malloc_ (size_t sz) {
55 void *res;
56 if (sz==0) return NULL;
57 res = malloc(sz);
58 if (!(res)) {
59 fprintf (stderr, "%s, %i (%s):\n%s\n",
60 __FILE__, __LINE__, __func__,
61 "malloc() failed");
62 exit (1);
63 }
64 return res;
65 }
66
67/* #define ALLOC(ptr,type,num) \ */
68/* do { (ptr)=(type *)util_malloc_((num)*sizeof(type)); } while (0) */
69
70#define RALLOC(type,num) \
71 ((type *)util_malloc_((num)*sizeof(type)))
72
73#define DEALLOC(ptr) \
74 { if ((ptr) != 0) free (ptr); (ptr) = NULL; }
75
76 void util_free_ (void *ptr) {
77 if ((ptr) != NULL) {
78 free(ptr);
79 }
80 }
81
82#define SWAP(a,b,type) \
83 do { type tmp_=(a); (a)=(b); (b)=tmp_; } while(0)
84/******************************************************************************/
85
86 typedef struct {
87 double r,i;
88 } cmplx;
89
90#include "fftpack.h"
91
92#define WA(x,i) wa[(i)+(x)*ido]
93#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))]
94#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))]
95#define PM(a,b,c,d) { a=c+d; b=c-d; }
96#define PMC(a,b,c,d) { a.r=c.r+d.r; a.i=c.i+d.i; b.r=c.r-d.r; b.i=c.i-d.i; }
97#define ADDC(a,b,c) { a.r=b.r+c.r; a.i=b.i+c.i; }
98#define SCALEC(a,b) { a.r*=b; a.i*=b; }
99#define CONJFLIPC(a) { double tmp_=a.r; a.r=-a.i; a.i=tmp_; }
100/* (a+ib) = conj(c+id) * (e+if) */
101#define MULPM(a,b,c,d,e,f) { a=c*e+d*f; b=c*f-d*e; }
102
103#define CONCAT(a,b) a ## b
104
105#define X(arg) CONCAT(passb,arg)
106#define BACKWARD
107#include "FFT/fftpack_inc.c"
108#undef BACKWARD
109#undef X
110
111#define X(arg) CONCAT(passf,arg)
112#include "fftpack_inc.c"
113#undef X
114
115#undef CC
116#undef CH
117#define CC(a,b,c) cc[(a)+ido*((b)+l1*(c))]
118#define CH(a,b,c) ch[(a)+ido*((b)+cdim*(c))]
119
120 static void radf2 (size_t ido, size_t l1, const double *cc, double *ch,
121 const double *wa)
122 {
123 const size_t cdim=2;
124 size_t i, k, ic;
125 double ti2, tr2;
126
127 for (k=0; k<l1; k++)
128 PM (CH(0,0,k),CH(ido-1,1,k),CC(0,k,0),CC(0,k,1))
129 if ((ido&1)==0)
130 for (k=0; k<l1; k++)
131 {
132 CH( 0,1,k) = -CC(ido-1,k,1);
133 CH(ido-1,0,k) = CC(ido-1,k,0);
134 }
135 if (ido<=2) return;
136 for (k=0; k<l1; k++)
137 for (i=2; i<ido; i+=2)
138 {
139 ic=ido-i;
140 MULPM (tr2,ti2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1))
141 PM (CH(i-1,0,k),CH(ic-1,1,k),CC(i-1,k,0),tr2)
142 PM (CH(i ,0,k),CH(ic ,1,k),ti2,CC(i ,k,0))
143 }
144 }
145
146 static void radf3(size_t ido, size_t l1, const double *cc, double *ch,
147 const double *wa)
148 {
149 const size_t cdim=3;
150 static const double taur=-0.5, taui=0.86602540378443864676;
151 size_t i, k, ic;
152 double ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
153
154 for (k=0; k<l1; k++)
155 {
156 cr2=CC(0,k,1)+CC(0,k,2);
157 CH(0,0,k) = CC(0,k,0)+cr2;
158 CH(0,2,k) = taui*(CC(0,k,2)-CC(0,k,1));
159 CH(ido-1,1,k) = CC(0,k,0)+taur*cr2;
160 }
161 if (ido==1) return;
162 for (k=0; k<l1; k++)
163 for (i=2; i<ido; i+=2)
164 {
165 ic=ido-i;
166 MULPM (dr2,di2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1))
167 MULPM (dr3,di3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2))
168 cr2=dr2+dr3;
169 ci2=di2+di3;
170 CH(i-1,0,k) = CC(i-1,k,0)+cr2;
171 CH(i ,0,k) = CC(i ,k,0)+ci2;
172 tr2 = CC(i-1,k,0)+taur*cr2;
173 ti2 = CC(i ,k,0)+taur*ci2;
174 tr3 = taui*(di2-di3);
175 ti3 = taui*(dr3-dr2);
176 PM(CH(i-1,2,k),CH(ic-1,1,k),tr2,tr3)
177 PM(CH(i ,2,k),CH(ic ,1,k),ti3,ti2)
178 }
179 }
180
181 static void radf4(size_t ido, size_t l1, const double *cc, double *ch,
182 const double *wa)
183 {
184 const size_t cdim=4;
185 static const double hsqt2=0.70710678118654752440;
186 size_t i, k, ic;
187 double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
188
189 for (k=0; k<l1; k++)
190 {
191 PM (tr1,CH(0,2,k),CC(0,k,3),CC(0,k,1))
192 PM (tr2,CH(ido-1,1,k),CC(0,k,0),CC(0,k,2))
193 PM (CH(0,0,k),CH(ido-1,3,k),tr2,tr1)
194 }
195 if ((ido&1)==0)
196 for (k=0; k<l1; k++)
197 {
198 ti1=-hsqt2*(CC(ido-1,k,1)+CC(ido-1,k,3));
199 tr1= hsqt2*(CC(ido-1,k,1)-CC(ido-1,k,3));
200 PM (CH(ido-1,0,k),CH(ido-1,2,k),CC(ido-1,k,0),tr1)
201 PM (CH( 0,3,k),CH( 0,1,k),ti1,CC(ido-1,k,2))
202 }
203 if (ido<=2) return;
204 for (k=0; k<l1; k++)
205 for (i=2; i<ido; i+=2)
206 {
207 ic=ido-i;
208 MULPM(cr2,ci2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1))
209 MULPM(cr3,ci3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2))
210 MULPM(cr4,ci4,WA(2,i-2),WA(2,i-1),CC(i-1,k,3),CC(i,k,3))
211 PM(tr1,tr4,cr4,cr2)
212 PM(ti1,ti4,ci2,ci4)
213 PM(tr2,tr3,CC(i-1,k,0),cr3)
214 PM(ti2,ti3,CC(i ,k,0),ci3)
215 PM(CH(i-1,0,k),CH(ic-1,3,k),tr2,tr1)
216 PM(CH(i ,0,k),CH(ic ,3,k),ti1,ti2)
217 PM(CH(i-1,2,k),CH(ic-1,1,k),tr3,ti4)
218 PM(CH(i ,2,k),CH(ic ,1,k),tr4,ti3)
219 }
220 }
221
222 static void radf5(size_t ido, size_t l1, const double *cc, double *ch,
223 const double *wa)
224 {
225 const size_t cdim=5;
226 static const double tr11= 0.3090169943749474241, ti11=0.95105651629515357212,
227 tr12=-0.8090169943749474241, ti12=0.58778525229247312917;
228 size_t i, k, ic;
229 double ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3,
230 dr4, dr5, cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
231
232 for (k=0; k<l1; k++)
233 {
234 PM (cr2,ci5,CC(0,k,4),CC(0,k,1))
235 PM (cr3,ci4,CC(0,k,3),CC(0,k,2))
236 CH(0,0,k)=CC(0,k,0)+cr2+cr3;
237 CH(ido-1,1,k)=CC(0,k,0)+tr11*cr2+tr12*cr3;
238 CH(0,2,k)=ti11*ci5+ti12*ci4;
239 CH(ido-1,3,k)=CC(0,k,0)+tr12*cr2+tr11*cr3;
240 CH(0,4,k)=ti12*ci5-ti11*ci4;
241 }
242 if (ido==1) return;
243 for (k=0; k<l1;++k)
244 for (i=2; i<ido; i+=2)
245 {
246 ic=ido-i;
247 MULPM (dr2,di2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1))
248 MULPM (dr3,di3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2))
249 MULPM (dr4,di4,WA(2,i-2),WA(2,i-1),CC(i-1,k,3),CC(i,k,3))
250 MULPM (dr5,di5,WA(3,i-2),WA(3,i-1),CC(i-1,k,4),CC(i,k,4))
251 PM(cr2,ci5,dr5,dr2)
252 PM(ci2,cr5,di2,di5)
253 PM(cr3,ci4,dr4,dr3)
254 PM(ci3,cr4,di3,di4)
255 CH(i-1,0,k)=CC(i-1,k,0)+cr2+cr3;
256 CH(i ,0,k)=CC(i ,k,0)+ci2+ci3;
257 tr2=CC(i-1,k,0)+tr11*cr2+tr12*cr3;
258 ti2=CC(i ,k,0)+tr11*ci2+tr12*ci3;
259 tr3=CC(i-1,k,0)+tr12*cr2+tr11*cr3;
260 ti3=CC(i ,k,0)+tr12*ci2+tr11*ci3;
261 MULPM(tr5,tr4,cr5,cr4,ti11,ti12)
262 MULPM(ti5,ti4,ci5,ci4,ti11,ti12)
263 PM(CH(i-1,2,k),CH(ic-1,1,k),tr2,tr5)
264 PM(CH(i ,2,k),CH(ic ,1,k),ti5,ti2)
265 PM(CH(i-1,4,k),CH(ic-1,3,k),tr3,tr4)
266 PM(CH(i ,4,k),CH(ic ,3,k),ti4,ti3)
267 }
268 }
269
270#undef CH
271#undef CC
272#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))]
273#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))]
274#define C1(a,b,c) cc[(a)+ido*((b)+l1*(c))]
275#define C2(a,b) cc[(a)+idl1*(b)]
276#define CH2(a,b) ch[(a)+idl1*(b)]
277 static void radfg(size_t ido, size_t ip, size_t l1, size_t idl1,
278 double *cc, double *ch, const double *wa)
279 {
280 const size_t cdim=ip;
281 static const double twopi=6.28318530717958647692;
282 size_t idij, ipph, i, j, k, l, j2, ic, jc, lc, ik;
283 double ai1, ai2, ar1, ar2, arg;
284 double *csarr;
285 size_t aidx;
286
287 ipph=(ip+1)/ 2;
288 if(ido!=1)
289 {
290 std::memcpy(ch,cc,idl1*sizeof(double));
291
292 for(j=1; j<ip; j++)
293 for(k=0; k<l1; k++)
294 {
295 CH(0,k,j)=C1(0,k,j);
296 idij=(j-1)*ido+1;
297 for(i=2; i<ido; i+=2,idij+=2)
298 MULPM(CH(i-1,k,j),CH(i,k,j),wa[idij-1],wa[idij],C1(i-1,k,j),C1(i,k,j))
299 }
300
301 for(j=1,jc=ip-1; j<ipph; j++,jc--)
302 for(k=0; k<l1; k++)
303 for(i=2; i<ido; i+=2)
304 {
305 PM(C1(i-1,k,j),C1(i ,k,jc),CH(i-1,k,jc),CH(i-1,k,j ))
306 PM(C1(i ,k,j),C1(i-1,k,jc),CH(i ,k,j ),CH(i ,k,jc))
307 }
308 }
309 else
310 std::memcpy(cc,ch,idl1*sizeof(double));
311
312 for(j=1,jc=ip-1; j<ipph; j++,jc--)
313 for(k=0; k<l1; k++)
314 PM(C1(0,k,j),C1(0,k,jc),CH(0,k,jc),CH(0,k,j))
315
316 csarr=RALLOC(double,2*ip);
317 arg=twopi / ip;
318 csarr[0]=1.;
319 csarr[1]=0.;
320 csarr[2]=csarr[2*ip-2]=cos(arg);
321 csarr[3]=sin(arg); csarr[2*ip-1]=-csarr[3];
322 for (i=2; i<=ip/2; ++i)
323 {
324 csarr[2*i]=csarr[2*ip-2*i]=cos(i*arg);
325 csarr[2*i+1]=sin(i*arg);
326 csarr[2*ip-2*i+1]=-csarr[2*i+1];
327 }
328 for(l=1,lc=ip-1; l<ipph; l++,lc--)
329 {
330 ar1=csarr[2*l];
331 ai1=csarr[2*l+1];
332 for(ik=0; ik<idl1; ik++)
333 {
334 CH2(ik,l)=C2(ik,0)+ar1*C2(ik,1);
335 CH2(ik,lc)=ai1*C2(ik,ip-1);
336 }
337 aidx=2*l;
338 for(j=2,jc=ip-2; j<ipph; j++,jc--)
339 {
340 aidx+=2*l;
341 if (aidx>=2*ip) aidx-=2*ip;
342 ar2=csarr[aidx];
343 ai2=csarr[aidx+1];
344 for(ik=0; ik<idl1; ik++)
345 {
346 CH2(ik,l )+=ar2*C2(ik,j );
347 CH2(ik,lc)+=ai2*C2(ik,jc);
348 }
349 }
350 }
351 DEALLOC(csarr);
352
353 for(j=1; j<ipph; j++)
354 for(ik=0; ik<idl1; ik++)
355 CH2(ik,0)+=C2(ik,j);
356
357 for(k=0; k<l1; k++)
358 std::memcpy(&CC(0,0,k),&CH(0,k,0),ido*sizeof(double));
359 for(j=1; j<ipph; j++)
360 {
361 jc=ip-j;
362 j2=2*j;
363 for(k=0; k<l1; k++)
364 {
365 CC(ido-1,j2-1,k) = CH(0,k,j );
366 CC(0 ,j2 ,k) = CH(0,k,jc);
367 }
368 }
369 if(ido==1) return;
370
371 for(j=1; j<ipph; j++)
372 {
373 jc=ip-j;
374 j2=2*j;
375 for(k=0; k<l1; k++)
376 for(i=2; i<ido; i+=2)
377 {
378 ic=ido-i;
379 PM (CC(i-1,j2,k),CC(ic-1,j2-1,k),CH(i-1,k,j ),CH(i-1,k,jc))
380 PM (CC(i ,j2,k),CC(ic ,j2-1,k),CH(i ,k,jc),CH(i ,k,j ))
381 }
382 }
383 }
384
385#undef CC
386#undef CH
387#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))]
388#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))]
389
390 static void radb2(size_t ido, size_t l1, const double *cc, double *ch,
391 const double *wa)
392 {
393 const size_t cdim=2;
394 size_t i, k, ic;
395 double ti2, tr2;
396
397 for (k=0; k<l1; k++)
398 PM (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(ido-1,1,k))
399 if ((ido&1)==0)
400 for (k=0; k<l1; k++)
401 {
402 CH(ido-1,k,0) = 2*CC(ido-1,0,k);
403 CH(ido-1,k,1) = -2*CC(0 ,1,k);
404 }
405 if (ido<=2) return;
406 for (k=0; k<l1;++k)
407 for (i=2; i<ido; i+=2)
408 {
409 ic=ido-i;
410 PM (CH(i-1,k,0),tr2,CC(i-1,0,k),CC(ic-1,1,k))
411 PM (ti2,CH(i ,k,0),CC(i ,0,k),CC(ic ,1,k))
412 MULPM (CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),ti2,tr2)
413 }
414 }
415
416 static void radb3(size_t ido, size_t l1, const double *cc, double *ch,
417 const double *wa)
418 {
419 const size_t cdim=3;
420 static const double taur=-0.5, taui=0.86602540378443864676;
421 size_t i, k, ic;
422 double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
423
424 for (k=0; k<l1; k++)
425 {
426 tr2=2*CC(ido-1,1,k);
427 cr2=CC(0,0,k)+taur*tr2;
428 CH(0,k,0)=CC(0,0,k)+tr2;
429 ci3=2*taui*CC(0,2,k);
430 PM (CH(0,k,2),CH(0,k,1),cr2,ci3);
431 }
432 if (ido==1) return;
433 for (k=0; k<l1; k++)
434 for (i=2; i<ido; i+=2)
435 {
436 ic=ido-i;
437 tr2=CC(i-1,2,k)+CC(ic-1,1,k);
438 ti2=CC(i ,2,k)-CC(ic ,1,k);
439 cr2=CC(i-1,0,k)+taur*tr2;
440 ci2=CC(i ,0,k)+taur*ti2;
441 CH(i-1,k,0)=CC(i-1,0,k)+tr2;
442 CH(i ,k,0)=CC(i ,0,k)+ti2;
443 cr3=taui*(CC(i-1,2,k)-CC(ic-1,1,k));
444 ci3=taui*(CC(i ,2,k)+CC(ic ,1,k));
445 PM(dr3,dr2,cr2,ci3)
446 PM(di2,di3,ci2,cr3)
447 MULPM(CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),di2,dr2)
448 MULPM(CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),di3,dr3)
449 }
450 }
451
452 static void radb4(size_t ido, size_t l1, const double *cc, double *ch,
453 const double *wa)
454 {
455 const size_t cdim=4;
456 static const double sqrt2=1.41421356237309504880;
457 size_t i, k, ic;
458 double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
459
460 for (k=0; k<l1; k++)
461 {
462 PM (tr2,tr1,CC(0,0,k),CC(ido-1,3,k))
463 tr3=2*CC(ido-1,1,k);
464 tr4=2*CC(0,2,k);
465 PM (CH(0,k,0),CH(0,k,2),tr2,tr3)
466 PM (CH(0,k,3),CH(0,k,1),tr1,tr4)
467 }
468 if ((ido&1)==0)
469 for (k=0; k<l1; k++)
470 {
471 PM (ti1,ti2,CC(0 ,3,k),CC(0 ,1,k))
472 PM (tr2,tr1,CC(ido-1,0,k),CC(ido-1,2,k))
473 CH(ido-1,k,0)=tr2+tr2;
474 CH(ido-1,k,1)=sqrt2*(tr1-ti1);
475 CH(ido-1,k,2)=ti2+ti2;
476 CH(ido-1,k,3)=-sqrt2*(tr1+ti1);
477 }
478 if (ido<=2) return;
479 for (k=0; k<l1;++k)
480 for (i=2; i<ido; i+=2)
481 {
482 ic=ido-i;
483 PM (tr2,tr1,CC(i-1,0,k),CC(ic-1,3,k))
484 PM (ti1,ti2,CC(i ,0,k),CC(ic ,3,k))
485 PM (tr4,ti3,CC(i ,2,k),CC(ic ,1,k))
486 PM (tr3,ti4,CC(i-1,2,k),CC(ic-1,1,k))
487 PM (CH(i-1,k,0),cr3,tr2,tr3)
488 PM (CH(i ,k,0),ci3,ti2,ti3)
489 PM (cr4,cr2,tr1,tr4)
490 PM (ci2,ci4,ti1,ti4)
491 MULPM (CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),ci2,cr2)
492 MULPM (CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),ci3,cr3)
493 MULPM (CH(i,k,3),CH(i-1,k,3),WA(2,i-2),WA(2,i-1),ci4,cr4)
494 }
495 }
496
497 static void radb5(size_t ido, size_t l1, const double *cc, double *ch,
498 const double *wa)
499 {
500 const size_t cdim=5;
501 static const double tr11= 0.3090169943749474241, ti11=0.95105651629515357212,
502 tr12=-0.8090169943749474241, ti12=0.58778525229247312917;
503 size_t i, k, ic;
504 double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
505 ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
506
507 for (k=0; k<l1; k++)
508 {
509 ti5=2*CC(0,2,k);
510 ti4=2*CC(0,4,k);
511 tr2=2*CC(ido-1,1,k);
512 tr3=2*CC(ido-1,3,k);
513 CH(0,k,0)=CC(0,0,k)+tr2+tr3;
514 cr2=CC(0,0,k)+tr11*tr2+tr12*tr3;
515 cr3=CC(0,0,k)+tr12*tr2+tr11*tr3;
516 MULPM(ci5,ci4,ti5,ti4,ti11,ti12)
517 PM(CH(0,k,4),CH(0,k,1),cr2,ci5)
518 PM(CH(0,k,3),CH(0,k,2),cr3,ci4)
519 }
520 if (ido==1) return;
521 for (k=0; k<l1;++k)
522 for (i=2; i<ido; i+=2)
523 {
524 ic=ido-i;
525 PM(tr2,tr5,CC(i-1,2,k),CC(ic-1,1,k))
526 PM(ti5,ti2,CC(i ,2,k),CC(ic ,1,k))
527 PM(tr3,tr4,CC(i-1,4,k),CC(ic-1,3,k))
528 PM(ti4,ti3,CC(i ,4,k),CC(ic ,3,k))
529 CH(i-1,k,0)=CC(i-1,0,k)+tr2+tr3;
530 CH(i ,k,0)=CC(i ,0,k)+ti2+ti3;
531 cr2=CC(i-1,0,k)+tr11*tr2+tr12*tr3;
532 ci2=CC(i ,0,k)+tr11*ti2+tr12*ti3;
533 cr3=CC(i-1,0,k)+tr12*tr2+tr11*tr3;
534 ci3=CC(i ,0,k)+tr12*ti2+tr11*ti3;
535 MULPM(cr5,cr4,tr5,tr4,ti11,ti12)
536 MULPM(ci5,ci4,ti5,ti4,ti11,ti12)
537 PM(dr4,dr3,cr3,ci4)
538 PM(di3,di4,ci3,cr4)
539 PM(dr5,dr2,cr2,ci5)
540 PM(di2,di5,ci2,cr5)
541 MULPM(CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),di2,dr2)
542 MULPM(CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),di3,dr3)
543 MULPM(CH(i,k,3),CH(i-1,k,3),WA(2,i-2),WA(2,i-1),di4,dr4)
544 MULPM(CH(i,k,4),CH(i-1,k,4),WA(3,i-2),WA(3,i-1),di5,dr5)
545 }
546 }
547
548 static void radbg(size_t ido, size_t ip, size_t l1, size_t idl1,
549 double *cc, double *ch, const double *wa)
550 {
551 const size_t cdim=ip;
552 static const double twopi=6.28318530717958647692;
553 size_t idij, ipph, i, j, k, l, j2, ic, jc, lc, ik;
554 double ai1, ai2, ar1, ar2, arg;
555 double *csarr;
556 size_t aidx;
557
558 ipph=(ip+1)/ 2;
559 for(k=0; k<l1; k++)
560 std::memcpy(&CH(0,k,0),&CC(0,0,k),ido*sizeof(double));
561 for(j=1; j<ipph; j++)
562 {
563 jc=ip-j;
564 j2=2*j;
565 for(k=0; k<l1; k++)
566 {
567 CH(0,k,j )=2*CC(ido-1,j2-1,k);
568 CH(0,k,jc)=2*CC(0 ,j2 ,k);
569 }
570 }
571
572 if(ido!=1)
573 for(j=1,jc=ip-1; j<ipph; j++,jc--)
574 for(k=0; k<l1; k++)
575 for(i=2; i<ido; i+=2)
576 {
577 ic=ido-i;
578 PM (CH(i-1,k,j ),CH(i-1,k,jc),CC(i-1,2*j,k),CC(ic-1,2*j-1,k))
579 PM (CH(i ,k,jc),CH(i ,k,j ),CC(i ,2*j,k),CC(ic ,2*j-1,k))
580 }
581
582 csarr=RALLOC(double,2*ip);
583 arg=twopi/ip;
584 csarr[0]=1.;
585 csarr[1]=0.;
586 csarr[2]=csarr[2*ip-2]=cos(arg);
587 csarr[3]=sin(arg); csarr[2*ip-1]=-csarr[3];
588 for (i=2; i<=ip/2; ++i)
589 {
590 csarr[2*i]=csarr[2*ip-2*i]=cos(i*arg);
591 csarr[2*i+1]=sin(i*arg);
592 csarr[2*ip-2*i+1]=-csarr[2*i+1];
593 }
594 for(l=1; l<ipph; l++)
595 {
596 lc=ip-l;
597 ar1=csarr[2*l];
598 ai1=csarr[2*l+1];
599 for(ik=0; ik<idl1; ik++)
600 {
601 C2(ik,l)=CH2(ik,0)+ar1*CH2(ik,1);
602 C2(ik,lc)=ai1*CH2(ik,ip-1);
603 }
604 aidx=2*l;
605 for(j=2; j<ipph; j++)
606 {
607 jc=ip-j;
608 aidx+=2*l;
609 if (aidx>=2*ip) aidx-=2*ip;
610 ar2=csarr[aidx];
611 ai2=csarr[aidx+1];
612 for(ik=0; ik<idl1; ik++)
613 {
614 C2(ik,l )+=ar2*CH2(ik,j );
615 C2(ik,lc)+=ai2*CH2(ik,jc);
616 }
617 }
618 }
619 DEALLOC(csarr);
620
621 for(j=1; j<ipph; j++)
622 for(ik=0; ik<idl1; ik++)
623 CH2(ik,0)+=CH2(ik,j);
624
625 for(j=1,jc=ip-1; j<ipph; j++,jc--)
626 for(k=0; k<l1; k++)
627 PM (CH(0,k,jc),CH(0,k,j),C1(0,k,j),C1(0,k,jc))
628
629 if(ido==1)
630 return;
631 for(j=1,jc=ip-1; j<ipph; j++,jc--)
632 for(k=0; k<l1; k++)
633 for(i=2; i<ido; i+=2)
634 {
635 PM (CH(i-1,k,jc),CH(i-1,k,j ),C1(i-1,k,j),C1(i ,k,jc))
636 PM (CH(i ,k,j ),CH(i ,k,jc),C1(i ,k,j),C1(i-1,k,jc))
637 }
638 std::memcpy(cc,ch,idl1*sizeof(double));
639
640 for(j=1; j<ip; j++)
641 for(k=0; k<l1; k++)
642 {
643 C1(0,k,j)=CH(0,k,j);
644 idij=(j-1)*ido+1;
645 for(i=2; i<ido; i+=2,idij+=2)
646 MULPM (C1(i,k,j),C1(i-1,k,j),wa[idij-1],wa[idij],CH(i,k,j),CH(i-1,k,j))
647 }
648 }
649
650#undef CC
651#undef CH
652#undef PM
653#undef MULPM
654
655
656/*----------------------------------------------------------------------
657 cfftf1, cfftb1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
658 ----------------------------------------------------------------------*/
659
660 static void cfft1(size_t n, cmplx c[], cmplx ch[], const cmplx wa[],
661 const size_t ifac[], int isign)
662 {
663 size_t k1, l1=1, nf=ifac[1], iw=0;
664 cmplx *p1=c, *p2=ch;
665
666 for(k1=0; k1<nf; k1++)
667 {
668 size_t ip=ifac[k1+2];
669 size_t l2=ip*l1;
670 size_t ido = n/l2;
671 if(ip==4)
672 (isign>0) ? passb4(ido, l1, p1, p2, wa+iw)
673 : passf4(ido, l1, p1, p2, wa+iw);
674 else if(ip==2)
675 (isign>0) ? passb2(ido, l1, p1, p2, wa+iw)
676 : passf2(ido, l1, p1, p2, wa+iw);
677 else if(ip==3)
678 (isign>0) ? passb3(ido, l1, p1, p2, wa+iw)
679 : passf3(ido, l1, p1, p2, wa+iw);
680 else if(ip==5)
681 (isign>0) ? passb5(ido, l1, p1, p2, wa+iw)
682 : passf5(ido, l1, p1, p2, wa+iw);
683 else if(ip==6)
684 (isign>0) ? passb6(ido, l1, p1, p2, wa+iw)
685 : passf6(ido, l1, p1, p2, wa+iw);
686 else
687 (isign>0) ? passbg(ido, ip, l1, p1, p2, wa+iw)
688 : passfg(ido, ip, l1, p1, p2, wa+iw);
689 SWAP(p1,p2,cmplx *);
690 l1=l2;
691 iw+=(ip-1)*ido;
692 }
693 if (p1!=c)
694 std::memcpy (c,p1,n*sizeof(cmplx));
695 }
696
697 void cfftf(size_t n, double c[], double wsave[])
698 {
699 if (n!=1)
700 cfft1(n, (cmplx*)c, (cmplx*)wsave, (cmplx*)(wsave+2*n),
701 (size_t*)(wsave+4*n),-1);
702 }
703
704 void cfftb(size_t n, double c[], double wsave[])
705 {
706 if (n!=1)
707 cfft1(n, (cmplx*)c, (cmplx*)wsave, (cmplx*)(wsave+2*n),
708 (size_t*)(wsave+4*n),+1);
709 }
710
711 static void factorize (size_t n, const size_t *pf, size_t npf, size_t *ifac)
712 {
713 size_t nl=n, nf=0, ntry=0, j=0, i;
714
715 startloop:
716 j++;
717 ntry = (j<=npf) ? pf[j-1] : ntry+2;
718 do
719 {
720 size_t nq=nl / ntry;
721 size_t nr=nl-ntry*nq;
722 if (nr!=0)
723 goto startloop;
724 nf++;
725 ifac[nf+1]=ntry;
726 nl=nq;
727 if ((ntry==2) && (nf!=1))
728 {
729 for (i=nf+1; i>2; --i)
730 ifac[i]=ifac[i-1];
731 ifac[2]=2;
732 }
733 }
734 while(nl!=1);
735 ifac[0]=n;
736 ifac[1]=nf;
737 }
738
739 static void cffti1(size_t n, double wa[], size_t ifac[])
740 {
741 static const size_t ntryh[5]={4,6,3,2,5};
742 static const double twopi=6.28318530717958647692;
743 size_t j, k, fi;
744
745 double argh=twopi/n;
746 size_t i=0, l1=1;
747 factorize (n,ntryh,5,ifac);
748 for(k=1; k<=ifac[1]; k++)
749 {
750 size_t ip=ifac[k+1];
751 size_t ido=n/(l1*ip);
752 for(j=1; j<ip; j++)
753 {
754 size_t is = i;
755 double argld=j*l1*argh;
756 wa[i ]=1;
757 wa[i+1]=0;
758 for(fi=1; fi<=ido; fi++)
759 {
760 double arg=fi*argld;
761 i+=2;
762 wa[i ]=cos(arg);
763 wa[i+1]=sin(arg);
764 }
765 if(ip>6)
766 {
767 wa[is ]=wa[i ];
768 wa[is+1]=wa[i+1];
769 }
770 }
771 l1*=ip;
772 }
773 }
774
775 void cffti(size_t n, double wsave[])
776 { if (n!=1) cffti1(n, wsave+2*n,(size_t*)(wsave+4*n)); }
777
778
779/*----------------------------------------------------------------------
780 rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Real FFTs.
781 ----------------------------------------------------------------------*/
782
783 static void rfftf1(
784 size_t n,
785 double c[],
786 double ch[],
787 const double wa[],
788 const size_t ifac[]
789 ) {
790 size_t k1, l1=n, nf=ifac[1], iw=n-1;
791 double *p1=ch, *p2=c;
792
793 for(k1=1; k1<=nf;++k1)
794 {
795 size_t ip=ifac[nf-k1+2];
796 size_t ido=n / l1;
797 l1 /= ip;
798 iw-=(ip-1)*ido;
799 SWAP (p1,p2,double *);
800 if(ip==4)
801 radf4(ido, l1, p1, p2, wa+iw);
802 else if(ip==2)
803 radf2(ido, l1, p1, p2, wa+iw);
804 else if(ip==3)
805 radf3(ido, l1, p1, p2, wa+iw);
806 else if(ip==5)
807 radf5(ido, l1, p1, p2, wa+iw);
808 else
809 {
810 if (ido==1)
811 SWAP (p1,p2,double *);
812 radfg(ido, ip, l1, ido*l1, p1, p2, wa+iw);
813 SWAP (p1,p2,double *);
814 }
815 }
816 if (p1==c)
817 std::memcpy (c,ch,n*sizeof(double));
818 }
819
820 static void rfftb1(size_t n, double c[], double ch[], const double wa[],
821 const size_t ifac[])
822 {
823 size_t k1, l1=1, nf=ifac[1], iw=0;
824 double *p1=c, *p2=ch;
825
826 for(k1=1; k1<=nf; k1++)
827 {
828 size_t ip = ifac[k1+1],
829 ido= n/(ip*l1);
830 if(ip==4)
831 radb4(ido, l1, p1, p2, wa+iw);
832 else if(ip==2)
833 radb2(ido, l1, p1, p2, wa+iw);
834 else if(ip==3)
835 radb3(ido, l1, p1, p2, wa+iw);
836 else if(ip==5)
837 radb5(ido, l1, p1, p2, wa+iw);
838 else
839 {
840 radbg(ido, ip, l1, ido*l1, p1, p2, wa+iw);
841 if (ido!=1)
842 SWAP (p1,p2,double *);
843 }
844 SWAP (p1,p2,double *);
845 l1*=ip;
846 iw+=(ip-1)*ido;
847 }
848 if (p1!=c)
849 std::memcpy (c,ch,n*sizeof(double));
850 }
851
852 void rfftf(size_t n, double r[], double wsave[])
853 { if(n!=1) rfftf1(n, r, wsave, wsave+n,(size_t*)(wsave+2*n)); }
854
855 void rfftb(size_t n, double r[], double wsave[])
856 { if(n!=1) rfftb1(n, r, wsave, wsave+n,(size_t*)(wsave+2*n)); }
857
858 static void rffti1(size_t n, double wa[], size_t ifac[])
859 {
860 static const size_t ntryh[4]={4,2,3,5};
861 static const double twopi=6.28318530717958647692;
862 size_t i, j, k, fi;
863
864 double argh=twopi/n;
865 size_t is=0, l1=1;
866 factorize (n,ntryh,4,ifac);
867 for (k=1; k<ifac[1]; k++)
868 {
869 size_t ip=ifac[k+1],
870 ido=n/(l1*ip);
871 for (j=1; j<ip; ++j)
872 {
873 double argld=j*l1*argh;
874 for(i=is,fi=1; i<=ido+is-3; i+=2,++fi)
875 {
876 double arg=fi*argld;
877 wa[i ]=cos(arg);
878 wa[i+1]=sin(arg);
879 }
880 is+=ido;
881 }
882 l1*=ip;
883 }
884 }
885
886 void rffti(size_t n, double wsave[])
887 {
888 if (n!=1) {
889 rffti1(n, wsave+n,(size_t*)(wsave+2*n));
890 }
891 }
892
893 const double pi = 3.14159265358979;
894 void sinti (size_t n, double wsave[])
895 {
896 if (n <= 1) return;
897 size_t ns2 = n / 2;
898 double dt = pi / double (n+1);
899 for (size_t k = 0; k < ns2; k++) {
900 wsave[k] = 2.0 * sin ((k+1) * dt);
901 }
902 rffti (n+1, wsave+ns2);
903 }
904
905 /*
906 with n = 20
907 +0 -> was wsave[0..9]
908 +10 -> xh wsave[10..30]
909 +31 -> x wsave[31..51]
910 +52 -> xxifac wsave[52..]
911 */
912
913 void sint1 (size_t n, double war[], double was[], double xh[], double x[], double xxifac[])
914 {
915 double sqrt3 = 1.73205080756888;
916
917 for (size_t i = 0; i < n; i++) {
918 xh[i] = war[i];
919 war[i] = x[i];
920 }
921
922 if (n < 2) {
923 xh[0] *= 2;
924 } else if (n == 2) {
925 double xhold = sqrt3*(xh[0]+xh[1]);
926 xh[1] = sqrt3*(xh[0]-xh[1]);
927 xh[0] = xhold;
928 } else {
929 size_t np1 = n + 1;
930 size_t ns2 = n / 2;
931 x[0] = 0.0;
932 for (size_t k = 0; k < ns2; k++) {
933 size_t kc = np1 - k;
934 double t1 = xh[k] - xh[kc];
935 double t2 = was[k] * (xh[k] + xh[kc]);
936 x[k+1] = t1 + t2;
937 x[kc+1] = t2 - t1;
938 }
939 size_t modn = n % 2;
940 if (modn != 0) {
941 x[ns2+2] = 4.0 * xh[ns2+1];
942 }
943 rfftf1 (np1, x, xh, war, (size_t*)xxifac);
944 xh[0] *= 0.5;
945 for (size_t i = 2; i < n; i += 2) {
946 xh[i-1] = -x[i];
947 xh[i] = xh[i-2] + x[i-1];
948 }
949 if (modn == 0) {
950 xh[n-1] = -x[n];
951 }
952 }
953 for (size_t i = 0; i < n; i++) {
954 x[i] = war[i];
955 war[i] = xh[i];
956 }
957 }
958
959 void sint (size_t n, double x[], double wsave[])
960 {
961 sint1 (
962 n,
963 x,
964 wsave,
965 wsave + n/2,
966 wsave + n/2 + n + 1,
967 wsave + n/2 + 2*n + 2
968 );
969 }
970
971#ifdef __cplusplus
972}
973#endif
nr
const int nr
Definition ClassicRandom.h:24
cos
Tps< T > cos(const Tps< T > &x)
Cosine.
Definition TpsMath.h:129
sin
Tps< T > sin(const Tps< T > &x)
Sine.
Definition TpsMath.h:111
PM
#define PM(a, b, c, d)
Definition fftpack.cpp:95
WA
#define WA(x, i)
Definition fftpack.cpp:92
C1
#define C1(a, b, c)
Definition fftpack.cpp:274
sinti
void sinti(size_t n, double wsave[])
Definition fftpack.cpp:894
pi
const double pi
Definition fftpack.cpp:893
RALLOC
#define RALLOC(type, num)
Definition fftpack.cpp:70
SWAP
#define SWAP(a, b, type)
Definition fftpack.cpp:82
CH
#define CH(a, b, c)
Definition fftpack.cpp:93
DEALLOC
#define DEALLOC(ptr)
Definition fftpack.cpp:73
cffti
void cffti(size_t n, double wsave[])
Definition fftpack.cpp:775
cfftf
void cfftf(size_t n, double c[], double wsave[])
Definition fftpack.cpp:697
sint
void sint(size_t n, double x[], double wsave[])
Definition fftpack.cpp:959
sint1
void sint1(size_t n, double war[], double was[], double xh[], double x[], double xxifac[])
Definition fftpack.cpp:913
CH2
#define CH2(a, b)
Definition fftpack.cpp:276
util_free_
void util_free_(void *ptr)
Definition fftpack.cpp:76
rffti
void rffti(size_t n, double wsave[])
Definition fftpack.cpp:886
util_malloc_
void * util_malloc_(size_t sz)
Definition fftpack.cpp:54
C2
#define C2(a, b)
Definition fftpack.cpp:275
rfftf
void rfftf(size_t n, double r[], double wsave[])
Definition fftpack.cpp:852
cfftb
void cfftb(size_t n, double c[], double wsave[])
Definition fftpack.cpp:704
MULPM
#define MULPM(a, b, c, d, e, f)
Definition fftpack.cpp:101
rfftb
void rfftb(size_t n, double r[], double wsave[])
Definition fftpack.cpp:855
CC
#define CC(a, b, c)
Definition fftpack.cpp:94
arg
arg(a))
Physics::c
constexpr double c
The velocity of light in m/s.
Definition Physics.h:45
cmplx::i
double i
Definition fftpack.cpp:87
Generated on Wed Sep 30 2026 07:14:49 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.8