?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed include statements such that double precision version of genborn.c
[gromacs/rigid-bodies.git] / src / mdlib / gmx_fft_fftpack.c
bloba834e0c7d24ed2210cfd34dbd4922c984c90b4cf
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * Gromacs 4.0 Copyright (c) 1991-2003
5 * David van der Spoel, Erik Lindahl, University of Groningen.
6 *
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
11 *
12 * To help us fund GROMACS development, we humbly ask that you cite
13 * the research papers on the package. Check out http://www.gromacs.org
14 *
15 * And Hey:
16 * Gnomes, ROck Monsters And Chili Sauce
17 */
18 #ifdef HAVE_CONFIG_H
19 #include <config.h>
20 #endif
21
22 #ifdef GMX_FFT_FFTPACK
23
24 #include <math.h>
25 #include <stdlib.h>
26 #include <string.h>
27 #include <errno.h>
28
29
30 #include "gmx_fft.h"
31 #include "gmx_fatal.h"
32
33
34 /** Contents of the FFTPACK fft datatype.
35 *
36 * FFTPACK only does 1d transforms, so we use a pointers to another fft for
37 * the transform in the next dimension.
38 * Thus, a 3d-structure contains a pointer to a 2d one, which in turns contains
39 * a pointer to a 1d. The 1d structure has next==NULL.
40 */
41 struct gmx_fft
42 {
43 int ndim; /**< Dimensions, including our subdimensions. */
44 int n; /**< Number of points in this dimension. */
45 int ifac[15]; /**< 15 bytes needed for cfft and rfft */
46 struct gmx_fft *next; /**< Pointer to next dimension, or NULL. */
47 real * work; /**< 1st 4n reserved for cfft, 1st 2n for rfft */
48 };
49
50 #include <math.h>
51 #include <stdio.h>
52
53
54
55 static void
56 fftpack_passf2(int ido,
57 int l1,
58 real cc[],
59 real ch[],
60 real wa1[],
61 int isign)
62 {
63 int i, k, ah, ac;
64 real ti2, tr2;
65
66 if (ido <= 2)
67 {
68 for (k=0; k<l1; k++)
69 {
70 ah = k*ido;
71 ac = 2*k*ido;
72 ch[ah] = cc[ac] + cc[ac + ido];
73 ch[ah + ido*l1] = cc[ac] - cc[ac + ido];
74 ch[ah+1] = cc[ac+1] + cc[ac + ido + 1];
75 ch[ah + ido*l1 + 1] = cc[ac+1] - cc[ac + ido + 1];
76 }
77 }
78 else
79 {
80 for (k=0; k<l1; k++)
81 {
82 for (i=0; i<ido-1; i+=2)
83 {
84 ah = i + k*ido;
85 ac = i + 2*k*ido;
86 ch[ah] = cc[ac] + cc[ac + ido];
87 tr2 = cc[ac] - cc[ac + ido];
88 ch[ah+1] = cc[ac+1] + cc[ac + 1 + ido];
89 ti2 = cc[ac+1] - cc[ac + 1 + ido];
90 ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2;
91 ch[ah+l1*ido] = wa1[i]*tr2 - isign*wa1[i+1]*ti2;
92 }
93 }
94 }
95 }
96
97
98
99 static void
100 fftpack_passf3(int ido,
101 int l1,
102 real cc[],
103 real ch[],
104 real wa1[],
105 real wa2[],
106 int isign)
107 {
108 const real taur = -0.5;
109 const real taui = 0.866025403784439;
110
111 int i, k, ac, ah;
112 real ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
113
114 if (ido == 2)
115 {
116 for (k=1; k<=l1; k++)
117 {
118 ac = (3*k - 2)*ido;
119 tr2 = cc[ac] + cc[ac + ido];
120 cr2 = cc[ac - ido] + taur*tr2;
121 ah = (k - 1)*ido;
122 ch[ah] = cc[ac - ido] + tr2;
123
124 ti2 = cc[ac + 1] + cc[ac + ido + 1];
125 ci2 = cc[ac - ido + 1] + taur*ti2;
126 ch[ah + 1] = cc[ac - ido + 1] + ti2;
127
128 cr3 = isign*taui*(cc[ac] - cc[ac + ido]);
129 ci3 = isign*taui*(cc[ac + 1] - cc[ac + ido + 1]);
130 ch[ah + l1*ido] = cr2 - ci3;
131 ch[ah + 2*l1*ido] = cr2 + ci3;
132 ch[ah + l1*ido + 1] = ci2 + cr3;
133 ch[ah + 2*l1*ido + 1] = ci2 - cr3;
134 }
135 }
136 else
137 {
138 for (k=1; k<=l1; k++)
139 {
140 for (i=0; i<ido-1; i+=2)
141 {
142 ac = i + (3*k - 2)*ido;
143 tr2 = cc[ac] + cc[ac + ido];
144 cr2 = cc[ac - ido] + taur*tr2;
145 ah = i + (k-1)*ido;
146 ch[ah] = cc[ac - ido] + tr2;
147 ti2 = cc[ac + 1] + cc[ac + ido + 1];
148 ci2 = cc[ac - ido + 1] + taur*ti2;
149 ch[ah + 1] = cc[ac - ido + 1] + ti2;
150 cr3 = isign*taui*(cc[ac] - cc[ac + ido]);
151 ci3 = isign*taui*(cc[ac + 1] - cc[ac + ido + 1]);
152 dr2 = cr2 - ci3;
153 dr3 = cr2 + ci3;
154 di2 = ci2 + cr3;
155 di3 = ci2 - cr3;
156 ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
157 ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
158 ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
159 ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
160 }
161 }
162 }
163 }
164
165
166 static void
167 fftpack_passf4(int ido,
168 int l1,
169 real cc[],
170 real ch[],
171 real wa1[],
172 real wa2[],
173 real wa3[],
174 int isign)
175 {
176 int i, k, ac, ah;
177 real ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
178
179 if (ido == 2)
180 {
181 for (k=0; k<l1; k++)
182 {
183 ac = 4*k*ido + 1;
184 ti1 = cc[ac] - cc[ac + 2*ido];
185 ti2 = cc[ac] + cc[ac + 2*ido];
186 tr4 = cc[ac + 3*ido] - cc[ac + ido];
187 ti3 = cc[ac + ido] + cc[ac + 3*ido];
188 tr1 = cc[ac - 1] - cc[ac + 2*ido - 1];
189 tr2 = cc[ac - 1] + cc[ac + 2*ido - 1];
190 ti4 = cc[ac + ido - 1] - cc[ac + 3*ido - 1];
191 tr3 = cc[ac + ido - 1] + cc[ac + 3*ido - 1];
192 ah = k*ido;
193 ch[ah] = tr2 + tr3;
194 ch[ah + 2*l1*ido] = tr2 - tr3;
195 ch[ah + 1] = ti2 + ti3;
196 ch[ah + 2*l1*ido + 1] = ti2 - ti3;
197 ch[ah + l1*ido] = tr1 + isign*tr4;
198 ch[ah + 3*l1*ido] = tr1 - isign*tr4;
199 ch[ah + l1*ido + 1] = ti1 + isign*ti4;
200 ch[ah + 3*l1*ido + 1] = ti1 - isign*ti4;
201 }
202 }
203 else
204 {
205 for (k=0; k<l1; k++)
206 {
207 for (i=0; i<ido-1; i+=2)
208 {
209 ac = i + 1 + 4*k*ido;
210 ti1 = cc[ac] - cc[ac + 2*ido];
211 ti2 = cc[ac] + cc[ac + 2*ido];
212 ti3 = cc[ac + ido] + cc[ac + 3*ido];
213 tr4 = cc[ac + 3*ido] - cc[ac + ido];
214 tr1 = cc[ac - 1] - cc[ac + 2*ido - 1];
215 tr2 = cc[ac - 1] + cc[ac + 2*ido - 1];
216 ti4 = cc[ac + ido - 1] - cc[ac + 3*ido - 1];
217 tr3 = cc[ac + ido - 1] + cc[ac + 3*ido - 1];
218 ah = i + k*ido;
219 ch[ah] = tr2 + tr3;
220 cr3 = tr2 - tr3;
221 ch[ah + 1] = ti2 + ti3;
222 ci3 = ti2 - ti3;
223 cr2 = tr1 + isign*tr4;
224 cr4 = tr1 - isign*tr4;
225 ci2 = ti1 + isign*ti4;
226 ci4 = ti1 - isign*ti4;
227 ch[ah + l1*ido] = wa1[i]*cr2 - isign*wa1[i + 1]*ci2;
228 ch[ah + l1*ido + 1] = wa1[i]*ci2 + isign*wa1[i + 1]*cr2;
229 ch[ah + 2*l1*ido] = wa2[i]*cr3 - isign*wa2[i + 1]*ci3;
230 ch[ah + 2*l1*ido + 1] = wa2[i]*ci3 + isign*wa2[i + 1]*cr3;
231 ch[ah + 3*l1*ido] = wa3[i]*cr4 -isign*wa3[i + 1]*ci4;
232 ch[ah + 3*l1*ido + 1] = wa3[i]*ci4 + isign*wa3[i + 1]*cr4;
233 }
234 }
235 }
236 }
237
238
239 static void
240 fftpack_passf5(int ido,
241 int l1,
242 real cc[],
243 real ch[],
244 real wa1[],
245 real wa2[],
246 real wa3[],
247 real wa4[],
248 int isign)
249 {
250 const real tr11 = 0.309016994374947;
251 const real ti11 = 0.951056516295154;
252 const real tr12 = -0.809016994374947;
253 const real ti12 = 0.587785252292473;
254
255 int i, k, ac, ah;
256 real ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
257 ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
258
259 if (ido == 2)
260 {
261 for (k = 1; k <= l1; ++k)
262 {
263 ac = (5*k - 4)*ido + 1;
264 ti5 = cc[ac] - cc[ac + 3*ido];
265 ti2 = cc[ac] + cc[ac + 3*ido];
266 ti4 = cc[ac + ido] - cc[ac + 2*ido];
267 ti3 = cc[ac + ido] + cc[ac + 2*ido];
268 tr5 = cc[ac - 1] - cc[ac + 3*ido - 1];
269 tr2 = cc[ac - 1] + cc[ac + 3*ido - 1];
270 tr4 = cc[ac + ido - 1] - cc[ac + 2*ido - 1];
271 tr3 = cc[ac + ido - 1] + cc[ac + 2*ido - 1];
272 ah = (k - 1)*ido;
273 ch[ah] = cc[ac - ido - 1] + tr2 + tr3;
274 ch[ah + 1] = cc[ac - ido] + ti2 + ti3;
275 cr2 = cc[ac - ido - 1] + tr11*tr2 + tr12*tr3;
276 ci2 = cc[ac - ido] + tr11*ti2 + tr12*ti3;
277 cr3 = cc[ac - ido - 1] + tr12*tr2 + tr11*tr3;
278 ci3 = cc[ac - ido] + tr12*ti2 + tr11*ti3;
279 cr5 = isign*(ti11*tr5 + ti12*tr4);
280 ci5 = isign*(ti11*ti5 + ti12*ti4);
281 cr4 = isign*(ti12*tr5 - ti11*tr4);
282 ci4 = isign*(ti12*ti5 - ti11*ti4);
283 ch[ah + l1*ido] = cr2 - ci5;
284 ch[ah + 4*l1*ido] = cr2 + ci5;
285 ch[ah + l1*ido + 1] = ci2 + cr5;
286 ch[ah + 2*l1*ido + 1] = ci3 + cr4;
287 ch[ah + 2*l1*ido] = cr3 - ci4;
288 ch[ah + 3*l1*ido] = cr3 + ci4;
289 ch[ah + 3*l1*ido + 1] = ci3 - cr4;
290 ch[ah + 4*l1*ido + 1] = ci2 - cr5;
291 }
292 }
293 else
294 {
295 for (k=1; k<=l1; k++)
296 {
297 for (i=0; i<ido-1; i+=2)
298 {
299 ac = i + 1 + (k*5 - 4)*ido;
300 ti5 = cc[ac] - cc[ac + 3*ido];
301 ti2 = cc[ac] + cc[ac + 3*ido];
302 ti4 = cc[ac + ido] - cc[ac + 2*ido];
303 ti3 = cc[ac + ido] + cc[ac + 2*ido];
304 tr5 = cc[ac - 1] - cc[ac + 3*ido - 1];
305 tr2 = cc[ac - 1] + cc[ac + 3*ido - 1];
306 tr4 = cc[ac + ido - 1] - cc[ac + 2*ido - 1];
307 tr3 = cc[ac + ido - 1] + cc[ac + 2*ido - 1];
308 ah = i + (k - 1)*ido;
309 ch[ah] = cc[ac - ido - 1] + tr2 + tr3;
310 ch[ah + 1] = cc[ac - ido] + ti2 + ti3;
311 cr2 = cc[ac - ido - 1] + tr11*tr2 + tr12*tr3;
312 ci2 = cc[ac - ido] + tr11*ti2 + tr12*ti3;
313 cr3 = cc[ac - ido - 1] + tr12*tr2 + tr11*tr3;
314 ci3 = cc[ac - ido] + tr12*ti2 + tr11*ti3;
315 cr5 = isign*(ti11*tr5 + ti12*tr4);
316 ci5 = isign*(ti11*ti5 + ti12*ti4);
317 cr4 = isign*(ti12*tr5 - ti11*tr4);
318 ci4 = isign*(ti12*ti5 - ti11*ti4);
319 dr3 = cr3 - ci4;
320 dr4 = cr3 + ci4;
321 di3 = ci3 + cr4;
322 di4 = ci3 - cr4;
323 dr5 = cr2 + ci5;
324 dr2 = cr2 - ci5;
325 di5 = ci2 - cr5;
326 di2 = ci2 + cr5;
327 ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
328 ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
329 ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
330 ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
331 ch[ah + 3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4;
332 ch[ah + 3*l1*ido + 1] = wa3[i]*di4 + isign*wa3[i+1]*dr4;
333 ch[ah + 4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5;
334 ch[ah + 4*l1*ido + 1] = wa4[i]*di5 + isign*wa4[i+1]*dr5;
335 }
336 }
337 }
338 }
339
340
341 static void
342 fftpack_passf(int * nac,
343 int ido,
344 int ip,
345 int l1,
346 int idl1,
347 real cc[],
348 real ch[],
349 real wa[],
350 int isign)
351 {
352 int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj, idl, inc,idp;
353 real wai, war;
354
355 idot = ido / 2;
356 nt = ip*idl1;
357 ipph = (ip + 1) / 2;
358 idp = ip*ido;
359 if (ido >= l1)
360 {
361 for (j=1; j<ipph; j++)
362 {
363 jc = ip - j;
364 for (k=0; k<l1; k++)
365 {
366 for (i=0; i<ido; i++)
367 {
368 ch[i + (k + j*l1)*ido] = cc[i + (j + k*ip)*ido] + cc[i + (jc + k*ip)*ido];
369 ch[i + (k + jc*l1)*ido] = cc[i + (j + k*ip)*ido] - cc[i + (jc + k*ip)*ido];
370 }
371 }
372 }
373 for (k=0; k<l1; k++)
374 for (i=0; i<ido; i++)
375 ch[i + k*ido] = cc[i + k*ip*ido];
376 }
377 else
378 {
379 for (j=1; j<ipph; j++)
380 {
381 jc = ip - j;
382 for (i=0; i<ido; i++)
383 {
384 for (k=0; k<l1; k++)
385 {
386 ch[i + (k + j*l1)*ido] = cc[i + (j + k*ip)*ido] + cc[i + (jc + k*ip)*ido];
387 ch[i + (k + jc*l1)*ido] = cc[i + (j + k*ip)*ido] - cc[i + (jc + k*ip)*ido];
388 }
389 }
390 }
391 for (i=0; i<ido; i++)
392 for (k=0; k<l1; k++)
393 ch[i + k*ido] = cc[i + k*ip*ido];
394 }
395
396 idl = 2 - ido;
397 inc = 0;
398 for (l=1; l<ipph; l++)
399 {
400 lc = ip - l;
401 idl += ido;
402 for (ik=0; ik<idl1; ik++)
403 {
404 cc[ik + l*idl1] = ch[ik] + wa[idl - 2]*ch[ik + idl1];
405 cc[ik + lc*idl1] = isign*wa[idl-1]*ch[ik + (ip-1)*idl1];
406 }
407 idlj = idl;
408 inc += ido;
409 for (j=2; j<ipph; j++)
410 {
411 jc = ip - j;
412 idlj += inc;
413 if (idlj > idp) idlj -= idp;
414 war = wa[idlj - 2];
415 wai = wa[idlj-1];
416 for (ik=0; ik<idl1; ik++)
417 {
418 cc[ik + l*idl1] += war*ch[ik + j*idl1];
419 cc[ik + lc*idl1] += isign*wai*ch[ik + jc*idl1];
420 }
421 }
422 }
423 for (j=1; j<ipph; j++)
424 for (ik=0; ik<idl1; ik++)
425 ch[ik] += ch[ik + j*idl1];
426 for (j=1; j<ipph; j++)
427 {
428 jc = ip - j;
429 for (ik=1; ik<idl1; ik+=2)
430 {
431 ch[ik - 1 + j*idl1] = cc[ik - 1 + j*idl1] - cc[ik + jc*idl1];
432 ch[ik - 1 + jc*idl1] = cc[ik - 1 + j*idl1] + cc[ik + jc*idl1];
433 ch[ik + j*idl1] = cc[ik + j*idl1] + cc[ik - 1 + jc*idl1];
434 ch[ik + jc*idl1] = cc[ik + j*idl1] - cc[ik - 1 + jc*idl1];
435 }
436 }
437 *nac = 1;
438 if (ido == 2)
439 return;
440 *nac = 0;
441 for (ik=0; ik<idl1; ik++)
442 {
443 cc[ik] = ch[ik];
444 }
445 for (j=1; j<ip; j++)
446 {
447 for (k=0; k<l1; k++)
448 {
449 cc[(k + j*l1)*ido + 0] = ch[(k + j*l1)*ido + 0];
450 cc[(k + j*l1)*ido + 1] = ch[(k + j*l1)*ido + 1];
451 }
452 }
453 if (idot <= l1)
454 {
455 idij = 0;
456 for (j=1; j<ip; j++)
457 {
458 idij += 2;
459 for (i=3; i<ido; i+=2)
460 {
461 idij += 2;
462 for (k=0; k<l1; k++)
463 {
464 cc[i - 1 + (k + j*l1)*ido] =
465 wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
466 isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
467 cc[i + (k + j*l1)*ido] =
468 wa[idij - 2]*ch[i + (k + j*l1)*ido] +
469 isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
470 }
471 }
472 }
473 }
474 else
475 {
476 idj = 2 - ido;
477 for (j=1; j<ip; j++)
478 {
479 idj += ido;
480 for (k = 0; k < l1; k++)
481 {
482 idij = idj;
483 for (i=3; i<ido; i+=2)
484 {
485 idij += 2;
486 cc[i - 1 + (k + j*l1)*ido] =
487 wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
488 isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
489 cc[i + (k + j*l1)*ido] =
490 wa[idij - 2]*ch[i + (k + j*l1)*ido] +
491 isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
492 }
493 }
494 }
495 }
496 }
497
498
499
500 static void
501 fftpack_radf2(int ido,
502 int l1,
503 real cc[],
504 real ch[],
505 real wa1[])
506 {
507 int i, k, ic;
508 real ti2, tr2;
509 for (k=0; k<l1; k++)
510 {
511 ch[2*k*ido] = cc[k*ido] + cc[(k + l1)*ido];
512 ch[(2*k+1)*ido + ido-1] = cc[k*ido] - cc[(k + l1)*ido];
513 }
514 if (ido < 2)
515 return;
516 if (ido != 2)
517 {
518 for (k=0; k<l1; k++)
519 {
520 for (i=2; i<ido; i+=2)
521 {
522 ic = ido - i;
523 tr2 = wa1[i - 2]*cc[i-1 + (k + l1)*ido] + wa1[i - 1]*cc[i + (k + l1)*ido];
524 ti2 = wa1[i - 2]*cc[i + (k + l1)*ido] - wa1[i - 1]*cc[i-1 + (k + l1)*ido];
525 ch[i + 2*k*ido] = cc[i + k*ido] + ti2;
526 ch[ic + (2*k+1)*ido] = ti2 - cc[i + k*ido];
527 ch[i - 1 + 2*k*ido] = cc[i - 1 + k*ido] + tr2;
528 ch[ic - 1 + (2*k+1)*ido] = cc[i - 1 + k*ido] - tr2;
529 }
530 }
531 if (ido % 2 == 1)
532 return;
533 }
534 for (k=0; k<l1; k++)
535 {
536 ch[(2*k+1)*ido] = -cc[ido-1 + (k + l1)*ido];
537 ch[ido-1 + 2*k*ido] = cc[ido-1 + k*ido];
538 }
539 }
540
541
542 static void
543 fftpack_radb2(int ido,
544 int l1,
545 real cc[],
546 real ch[],
547 real wa1[])
548 {
549 int i, k, ic;
550 real ti2, tr2;
551 for (k=0; k<l1; k++)
552 {
553 ch[k*ido] = cc[2*k*ido] + cc[ido-1 + (2*k+1)*ido];
554 ch[(k + l1)*ido] = cc[2*k*ido] - cc[ido-1 + (2*k+1)*ido];
555 }
556 if (ido < 2)
557 return;
558 if (ido != 2)
559 {
560 for (k = 0; k < l1; ++k)
561 {
562 for (i = 2; i < ido; i += 2)
563 {
564 ic = ido - i;
565 ch[i-1 + k*ido] = cc[i-1 + 2*k*ido] + cc[ic-1 + (2*k+1)*ido];
566 tr2 = cc[i-1 + 2*k*ido] - cc[ic-1 + (2*k+1)*ido];
567 ch[i + k*ido] = cc[i + 2*k*ido] - cc[ic + (2*k+1)*ido];
568 ti2 = cc[i + (2*k)*ido] + cc[ic + (2*k+1)*ido];
569 ch[i-1 + (k + l1)*ido] = wa1[i - 2]*tr2 - wa1[i - 1]*ti2;
570 ch[i + (k + l1)*ido] = wa1[i - 2]*ti2 + wa1[i - 1]*tr2;
571 }
572 }
573 if (ido % 2 == 1)
574 return;
575 }
576 for (k = 0; k < l1; k++)
577 {
578 ch[ido-1 + k*ido] = 2*cc[ido-1 + 2*k*ido];
579 ch[ido-1 + (k + l1)*ido] = -2*cc[(2*k+1)*ido];
580 }
581 }
582
583
584 static void
585 fftpack_radf3(int ido,
586 int l1,
587 real cc[],
588 real ch[],
589 real wa1[],
590 real wa2[])
591 {
592 const real taur = -0.5;
593 const real taui = 0.866025403784439;
594 int i, k, ic;
595 real ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
596
597 for (k=0; k<l1; k++)
598 {
599 cr2 = cc[(k + l1)*ido] + cc[(k + 2*l1)*ido];
600 ch[3*k*ido] = cc[k*ido] + cr2;
601 ch[(3*k+2)*ido] = taui*(cc[(k + l1*2)*ido] - cc[(k + l1)*ido]);
602 ch[ido-1 + (3*k + 1)*ido] = cc[k*ido] + taur*cr2;
603 }
604 if (ido == 1)
605 return;
606 for (k=0; k<l1; k++)
607 {
608 for (i=2; i<ido; i+=2)
609 {
610 ic = ido - i;
611 dr2 = wa1[i - 2]*cc[i - 1 + (k + l1)*ido] +wa1[i - 1]*cc[i + (k + l1)*ido];
612 di2 = wa1[i - 2]*cc[i + (k + l1)*ido] - wa1[i - 1]*cc[i - 1 + (k + l1)*ido];
613 dr3 = wa2[i - 2]*cc[i - 1 + (k + l1*2)*ido] + wa2[i - 1]*cc[i + (k + l1*2)*ido];
614 di3 = wa2[i - 2]*cc[i + (k + l1*2)*ido] - wa2[i - 1]*cc[i - 1 + (k + l1*2)*ido];
615 cr2 = dr2 + dr3;
616 ci2 = di2 + di3;
617 ch[i - 1 + 3*k*ido] = cc[i - 1 + k*ido] + cr2;
618 ch[i + 3*k*ido] = cc[i + k*ido] + ci2;
619 tr2 = cc[i - 1 + k*ido] + taur*cr2;
620 ti2 = cc[i + k*ido] + taur*ci2;
621 tr3 = taui*(di2 - di3);
622 ti3 = taui*(dr3 - dr2);
623 ch[i - 1 + (3*k + 2)*ido] = tr2 + tr3;
624 ch[ic - 1 + (3*k + 1)*ido] = tr2 - tr3;
625 ch[i + (3*k + 2)*ido] = ti2 + ti3;
626 ch[ic + (3*k + 1)*ido] = ti3 - ti2;
627 }
628 }
629 }
630
631
632 static void
633 fftpack_radb3(int ido,
634 int l1,
635 real cc[],
636 real ch[],
637 real wa1[],
638 real wa2[])
639 {
640 const real taur = -0.5;
641 const real taui = 0.866025403784439;
642 int i, k, ic;
643 real ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
644
645 for (k=0; k<l1; k++)
646 {
647 tr2 = 2*cc[ido-1 + (3*k + 1)*ido];
648 cr2 = cc[3*k*ido] + taur*tr2;
649 ch[k*ido] = cc[3*k*ido] + tr2;
650 ci3 = 2*taui*cc[(3*k + 2)*ido];
651 ch[(k + l1)*ido] = cr2 - ci3;
652 ch[(k + 2*l1)*ido] = cr2 + ci3;
653 }
654 if (ido == 1)
655 return;
656
657 for (k=0; k<l1; k++)
658 {
659 for (i=2; i<ido; i+=2)
660 {
661 ic = ido - i;
662 tr2 = cc[i - 1 + (3*k + 2)*ido] + cc[ic - 1 + (3*k + 1)*ido];
663 cr2 = cc[i - 1 + 3*k*ido] + taur*tr2;
664 ch[i - 1 + k*ido] = cc[i - 1 + 3*k*ido] + tr2;
665 ti2 = cc[i + (3*k + 2)*ido]- cc[ic + (3*k + 1)*ido];
666 ci2 = cc[i + 3*k*ido] + taur*ti2;
667 ch[i + k*ido] = cc[i + 3*k*ido] + ti2;
668 cr3 = taui*(cc[i - 1 + (3*k + 2)*ido] - cc[ic - 1 + (3*k + 1)*ido]);
669 ci3 = taui*(cc[i + (3*k + 2)*ido] + cc[ic + (3*k + 1)*ido]);
670 dr2 = cr2 - ci3;
671 dr3 = cr2 + ci3;
672 di2 = ci2 + cr3;
673 di3 = ci2 - cr3;
674 ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
675 ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
676 ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
677 ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
678 }
679 }
680 }
681
682
683 static void
684 fftpack_radf4(int ido,
685 int l1,
686 real cc[],
687 real ch[],
688 real wa1[],
689 real wa2[],
690 real wa3[])
691 {
692 const real hsqt2 = 0.7071067811865475;
693 int i, k, ic;
694 real ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
695
696 for (k=0; k<l1; k++)
697 {
698 tr1 = cc[(k + l1)*ido] + cc[(k + 3*l1)*ido];
699 tr2 = cc[k*ido] + cc[(k + 2*l1)*ido];
700 ch[4*k*ido] = tr1 + tr2;
701 ch[ido-1 + (4*k + 3)*ido] = tr2 - tr1;
702 ch[ido-1 + (4*k + 1)*ido] = cc[k*ido] - cc[(k + 2*l1)*ido];
703 ch[(4*k + 2)*ido] = cc[(k + 3*l1)*ido] - cc[(k + l1)*ido];
704 }
705 if (ido < 2)
706 return;
707 if (ido != 2)
708 {
709 for (k=0; k<l1; k++)
710 {
711 for (i=2; i<ido; i += 2)
712 {
713 ic = ido - i;
714 cr2 = wa1[i - 2]*cc[i - 1 + (k + l1)*ido] + wa1[i - 1]*cc[i + (k + l1)*ido];
715 ci2 = wa1[i - 2]*cc[i + (k + l1)*ido] - wa1[i - 1]*cc[i - 1 + (k + l1)*ido];
716 cr3 = wa2[i - 2]*cc[i - 1 + (k + 2*l1)*ido] + wa2[i - 1]*cc[i + (k + 2*l1)*ido];
717 ci3 = wa2[i - 2]*cc[i + (k + 2*l1)*ido] - wa2[i - 1]*cc[i - 1 + (k + 2*l1)*ido];
718 cr4 = wa3[i - 2]*cc[i - 1 + (k + 3*l1)*ido] + wa3[i - 1]*cc[i + (k + 3*l1)*ido];
719 ci4 = wa3[i - 2]*cc[i + (k + 3*l1)*ido] - wa3[i - 1]*cc[i - 1 + (k + 3*l1)*ido];
720 tr1 = cr2 + cr4;
721 tr4 = cr4 - cr2;
722 ti1 = ci2 + ci4;
723 ti4 = ci2 - ci4;
724 ti2 = cc[i + k*ido] + ci3;
725 ti3 = cc[i + k*ido] - ci3;
726 tr2 = cc[i - 1 + k*ido] + cr3;
727 tr3 = cc[i - 1 + k*ido] - cr3;
728 ch[i - 1 + 4*k*ido] = tr1 + tr2;
729 ch[ic - 1 + (4*k + 3)*ido] = tr2 - tr1;
730 ch[i + 4*k*ido] = ti1 + ti2;
731 ch[ic + (4*k + 3)*ido] = ti1 - ti2;
732 ch[i - 1 + (4*k + 2)*ido] = ti4 + tr3;
733 ch[ic - 1 + (4*k + 1)*ido] = tr3 - ti4;
734 ch[i + (4*k + 2)*ido] = tr4 + ti3;
735 ch[ic + (4*k + 1)*ido] = tr4 - ti3;
736 }
737 }
738 if (ido % 2 == 1)
739 return;
740 }
741 for (k=0; k<l1; k++)
742 {
743 ti1 = -hsqt2*(cc[ido-1 + (k + l1)*ido] + cc[ido-1 + (k + 3*l1)*ido]);
744 tr1 = hsqt2*(cc[ido-1 + (k + l1)*ido] - cc[ido-1 + (k + 3*l1)*ido]);
745 ch[ido-1 + 4*k*ido] = tr1 + cc[ido-1 + k*ido];
746 ch[ido-1 + (4*k + 2)*ido] = cc[ido-1 + k*ido] - tr1;
747 ch[(4*k + 1)*ido] = ti1 - cc[ido-1 + (k + 2*l1)*ido];
748 ch[(4*k + 3)*ido] = ti1 + cc[ido-1 + (k + 2*l1)*ido];
749 }
750 }
751
752
753 static void
754 fftpack_radb4(int ido,
755 int l1,
756 real cc[],
757 real ch[],
758 real wa1[],
759 real wa2[],
760 real wa3[])
761 {
762 const real sqrt2 = 1.414213562373095;
763 int i, k, ic;
764 real ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
765 for (k = 0; k < l1; k++)
766 {
767 tr1 = cc[4*k*ido] - cc[ido-1 + (4*k + 3)*ido];
768 tr2 = cc[4*k*ido] + cc[ido-1 + (4*k + 3)*ido];
769 tr3 = cc[ido-1 + (4*k + 1)*ido] + cc[ido-1 + (4*k + 1)*ido];
770 tr4 = cc[(4*k + 2)*ido] + cc[(4*k + 2)*ido];
771 ch[k*ido] = tr2 + tr3;
772 ch[(k + l1)*ido] = tr1 - tr4;
773 ch[(k + 2*l1)*ido] = tr2 - tr3;
774 ch[(k + 3*l1)*ido] = tr1 + tr4;
775 }
776 if (ido < 2)
777 return;
778 if (ido != 2)
779 {
780 for (k = 0; k < l1; ++k)
781 {
782 for (i = 2; i < ido; i += 2)
783 {
784 ic = ido - i;
785 ti1 = cc[i + 4*k*ido] + cc[ic + (4*k + 3)*ido];
786 ti2 = cc[i + 4*k*ido] - cc[ic + (4*k + 3)*ido];
787 ti3 = cc[i + (4*k + 2)*ido] - cc[ic + (4*k + 1)*ido];
788 tr4 = cc[i + (4*k + 2)*ido] + cc[ic + (4*k + 1)*ido];
789 tr1 = cc[i - 1 + 4*k*ido] - cc[ic - 1 + (4*k + 3)*ido];
790 tr2 = cc[i - 1 + 4*k*ido] + cc[ic - 1 + (4*k + 3)*ido];
791 ti4 = cc[i - 1 + (4*k + 2)*ido] - cc[ic - 1 + (4*k + 1)*ido];
792 tr3 = cc[i - 1 + (4*k + 2)*ido] + cc[ic - 1 + (4*k + 1)*ido];
793 ch[i - 1 + k*ido] = tr2 + tr3;
794 cr3 = tr2 - tr3;
795 ch[i + k*ido] = ti2 + ti3;
796 ci3 = ti2 - ti3;
797 cr2 = tr1 - tr4;
798 cr4 = tr1 + tr4;
799 ci2 = ti1 + ti4;
800 ci4 = ti1 - ti4;
801 ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*cr2 - wa1[i - 1]*ci2;
802 ch[i + (k + l1)*ido] = wa1[i - 2]*ci2 + wa1[i - 1]*cr2;
803 ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*cr3 - wa2[i - 1]*ci3;
804 ch[i + (k + 2*l1)*ido] = wa2[i - 2]*ci3 + wa2[i - 1]*cr3;
805 ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*cr4 - wa3[i - 1]*ci4;
806 ch[i + (k + 3*l1)*ido] = wa3[i - 2]*ci4 + wa3[i - 1]*cr4;
807 }
808 }
809 if (ido % 2 == 1)
810 return;
811 }
812 for (k = 0; k < l1; k++)
813 {
814 ti1 = cc[(4*k + 1)*ido] + cc[(4*k + 3)*ido];
815 ti2 = cc[(4*k + 3)*ido] - cc[(4*k + 1)*ido];
816 tr1 = cc[ido-1 + 4*k*ido] - cc[ido-1 + (4*k + 2)*ido];
817 tr2 = cc[ido-1 + 4*k*ido] + cc[ido-1 + (4*k + 2)*ido];
818 ch[ido-1 + k*ido] = tr2 + tr2;
819 ch[ido-1 + (k + l1)*ido] = sqrt2*(tr1 - ti1);
820 ch[ido-1 + (k + 2*l1)*ido] = ti2 + ti2;
821 ch[ido-1 + (k + 3*l1)*ido] = -sqrt2*(tr1 + ti1);
822 }
823 }
824
825
826 static void
827 fftpack_radf5(int ido,
828 int l1,
829 real cc[],
830 real ch[],
831 real wa1[],
832 real wa2[],
833 real wa3[],
834 real wa4[])
835 {
836 const real tr11 = 0.309016994374947;
837 const real ti11 = 0.951056516295154;
838 const real tr12 = -0.809016994374947;
839 const real ti12 = 0.587785252292473;
840 int i, k, ic;
841 real ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,
842 cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
843
844 for (k = 0; k < l1; k++)
845 {
846 cr2 = cc[(k + 4*l1)*ido] + cc[(k + l1)*ido];
847 ci5 = cc[(k + 4*l1)*ido] - cc[(k + l1)*ido];
848 cr3 = cc[(k + 3*l1)*ido] + cc[(k + 2*l1)*ido];
849 ci4 = cc[(k + 3*l1)*ido] - cc[(k + 2*l1)*ido];
850 ch[5*k*ido] = cc[k*ido] + cr2 + cr3;
851 ch[ido-1 + (5*k + 1)*ido] = cc[k*ido] + tr11*cr2 + tr12*cr3;
852 ch[(5*k + 2)*ido] = ti11*ci5 + ti12*ci4;
853 ch[ido-1 + (5*k + 3)*ido] = cc[k*ido] + tr12*cr2 + tr11*cr3;
854 ch[(5*k + 4)*ido] = ti12*ci5 - ti11*ci4;
855 }
856 if (ido == 1)
857 return;
858 for (k = 0; k < l1; ++k)
859 {
860 for (i = 2; i < ido; i += 2)
861 {
862 ic = ido - i;
863 dr2 = wa1[i - 2]*cc[i - 1 + (k + l1)*ido] + wa1[i - 1]*cc[i + (k + l1)*ido];
864 di2 = wa1[i - 2]*cc[i + (k + l1)*ido] - wa1[i - 1]*cc[i - 1 + (k + l1)*ido];
865 dr3 = wa2[i - 2]*cc[i - 1 + (k + 2*l1)*ido] + wa2[i - 1]*cc[i + (k + 2*l1)*ido];
866 di3 = wa2[i - 2]*cc[i + (k + 2*l1)*ido] - wa2[i - 1]*cc[i - 1 + (k + 2*l1)*ido];
867 dr4 = wa3[i - 2]*cc[i - 1 + (k + 3*l1)*ido] + wa3[i - 1]*cc[i + (k + 3*l1)*ido];
868 di4 = wa3[i - 2]*cc[i + (k + 3*l1)*ido] - wa3[i - 1]*cc[i - 1 + (k + 3*l1)*ido];
869 dr5 = wa4[i - 2]*cc[i - 1 + (k + 4*l1)*ido] + wa4[i - 1]*cc[i + (k + 4*l1)*ido];
870 di5 = wa4[i - 2]*cc[i + (k + 4*l1)*ido] - wa4[i - 1]*cc[i - 1 + (k + 4*l1)*ido];
871 cr2 = dr2 + dr5;
872 ci5 = dr5 - dr2;
873 cr5 = di2 - di5;
874 ci2 = di2 + di5;
875 cr3 = dr3 + dr4;
876 ci4 = dr4 - dr3;
877 cr4 = di3 - di4;
878 ci3 = di3 + di4;
879 ch[i - 1 + 5*k*ido] = cc[i - 1 + k*ido] + cr2 + cr3;
880 ch[i + 5*k*ido] = cc[i + k*ido] + ci2 + ci3;
881 tr2 = cc[i - 1 + k*ido] + tr11*cr2 + tr12*cr3;
882 ti2 = cc[i + k*ido] + tr11*ci2 + tr12*ci3;
883 tr3 = cc[i - 1 + k*ido] + tr12*cr2 + tr11*cr3;
884 ti3 = cc[i + k*ido] + tr12*ci2 + tr11*ci3;
885 tr5 = ti11*cr5 + ti12*cr4;
886 ti5 = ti11*ci5 + ti12*ci4;
887 tr4 = ti12*cr5 - ti11*cr4;
888 ti4 = ti12*ci5 - ti11*ci4;
889 ch[i - 1 + (5*k + 2)*ido] = tr2 + tr5;
890 ch[ic - 1 + (5*k + 1)*ido] = tr2 - tr5;
891 ch[i + (5*k + 2)*ido] = ti2 + ti5;
892 ch[ic + (5*k + 1)*ido] = ti5 - ti2;
893 ch[i - 1 + (5*k + 4)*ido] = tr3 + tr4;
894 ch[ic - 1 + (5*k + 3)*ido] = tr3 - tr4;
895 ch[i + (5*k + 4)*ido] = ti3 + ti4;
896 ch[ic + (5*k + 3)*ido] = ti4 - ti3;
897 }
898 }
899 }
900
901
902 static void
903 fftpack_radb5(int ido,
904 int l1,
905 real cc[],
906 real ch[],
907 real wa1[],
908 real wa2[],
909 real wa3[],
910 real wa4[])
911 {
912 const real tr11 = 0.309016994374947;
913 const real ti11 = 0.951056516295154;
914 const real tr12 = -0.809016994374947;
915 const real ti12 = 0.587785252292473;
916
917 int i, k, ic;
918 real ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
919 ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
920
921 for (k = 0; k < l1; k++)
922 {
923 ti5 = 2*cc[(5*k + 2)*ido];
924 ti4 = 2*cc[(5*k + 4)*ido];
925 tr2 = 2*cc[ido-1 + (5*k + 1)*ido];
926 tr3 = 2*cc[ido-1 + (5*k + 3)*ido];
927 ch[k*ido] = cc[5*k*ido] + tr2 + tr3;
928 cr2 = cc[5*k*ido] + tr11*tr2 + tr12*tr3;
929 cr3 = cc[5*k*ido] + tr12*tr2 + tr11*tr3;
930 ci5 = ti11*ti5 + ti12*ti4;
931 ci4 = ti12*ti5 - ti11*ti4;
932 ch[(k + l1)*ido] = cr2 - ci5;
933 ch[(k + 2*l1)*ido] = cr3 - ci4;
934 ch[(k + 3*l1)*ido] = cr3 + ci4;
935 ch[(k + 4*l1)*ido] = cr2 + ci5;
936 }
937 if (ido == 1) return;
938 for (k = 0; k < l1; ++k)
939 {
940 for (i = 2; i < ido; i += 2)
941 {
942 ic = ido - i;
943 ti5 = cc[i + (5*k + 2)*ido] + cc[ic + (5*k + 1)*ido];
944 ti2 = cc[i + (5*k + 2)*ido] - cc[ic + (5*k + 1)*ido];
945 ti4 = cc[i + (5*k + 4)*ido] + cc[ic + (5*k + 3)*ido];
946 ti3 = cc[i + (5*k + 4)*ido] - cc[ic + (5*k + 3)*ido];
947 tr5 = cc[i - 1 + (5*k + 2)*ido] - cc[ic - 1 + (5*k + 1)*ido];
948 tr2 = cc[i - 1 + (5*k + 2)*ido] + cc[ic - 1 + (5*k + 1)*ido];
949 tr4 = cc[i - 1 + (5*k + 4)*ido] - cc[ic - 1 + (5*k + 3)*ido];
950 tr3 = cc[i - 1 + (5*k + 4)*ido] + cc[ic - 1 + (5*k + 3)*ido];
951 ch[i - 1 + k*ido] = cc[i - 1 + 5*k*ido] + tr2 + tr3;
952 ch[i + k*ido] = cc[i + 5*k*ido] + ti2 + ti3;
953 cr2 = cc[i - 1 + 5*k*ido] + tr11*tr2 + tr12*tr3;
954 ci2 = cc[i + 5*k*ido] + tr11*ti2 + tr12*ti3;
955 cr3 = cc[i - 1 + 5*k*ido] + tr12*tr2 + tr11*tr3;
956 ci3 = cc[i + 5*k*ido] + tr12*ti2 + tr11*ti3;
957 cr5 = ti11*tr5 + ti12*tr4;
958 ci5 = ti11*ti5 + ti12*ti4;
959 cr4 = ti12*tr5 - ti11*tr4;
960 ci4 = ti12*ti5 - ti11*ti4;
961 dr3 = cr3 - ci4;
962 dr4 = cr3 + ci4;
963 di3 = ci3 + cr4;
964 di4 = ci3 - cr4;
965 dr5 = cr2 + ci5;
966 dr2 = cr2 - ci5;
967 di5 = ci2 - cr5;
968 di2 = ci2 + cr5;
969 ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
970 ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
971 ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
972 ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
973 ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*dr4 - wa3[i - 1]*di4;
974 ch[i + (k + 3*l1)*ido] = wa3[i - 2]*di4 + wa3[i - 1]*dr4;
975 ch[i - 1 + (k + 4*l1)*ido] = wa4[i - 2]*dr5 - wa4[i - 1]*di5;
976 ch[i + (k + 4*l1)*ido] = wa4[i - 2]*di5 + wa4[i - 1]*dr5;
977 }
978 }
979 }
980
981
982 static void
983 fftpack_radfg(int ido,
984 int ip,
985 int l1,
986 int idl1,
987 real cc[],
988 real ch[],
989 real wa[])
990 {
991 const real twopi = 6.28318530717959;
992 int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
993 real dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
994 arg = twopi / ip;
995 dcp = cos(arg);
996 dsp = sin(arg);
997 ipph = (ip + 1) / 2;
998 nbd = (ido - 1) / 2;
999 if (ido != 1)
1000 {
1001 for (ik=0; ik<idl1; ik++) ch[ik] = cc[ik];
1002 for (j=1; j<ip; j++)
1003 for (k=0; k<l1; k++)
1004 ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido];
1005 if (nbd <= l1)
1006 {
1007 is = -ido;
1008 for (j=1; j<ip; j++)
1009 {
1010 is += ido;
1011 idij = is-1;
1012 for (i=2; i<ido; i+=2)
1013 {
1014 idij += 2;
1015 for (k=0; k<l1; k++)
1016 {
1017 ch[i - 1 + (k + j*l1)*ido] =
1018 wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
1019 ch[i + (k + j*l1)*ido] =
1020 wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
1021 }
1022 }
1023 }
1024 }
1025 else
1026 {
1027 is = -ido;
1028 for (j=1; j<ip; j++)
1029 {
1030 is += ido;
1031 for (k=0; k<l1; k++)
1032 {
1033 idij = is-1;
1034 for (i=2; i<ido; i+=2)
1035 {
1036 idij += 2;
1037 ch[i - 1 + (k + j*l1)*ido] =
1038 wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
1039 ch[i + (k + j*l1)*ido] =
1040 wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
1041 }
1042 }
1043 }
1044 }
1045 if (nbd >= l1)
1046 {
1047 for (j=1; j<ipph; j++)
1048 {
1049 jc = ip - j;
1050 for (k=0; k<l1; k++)
1051 {
1052 for (i=2; i<ido; i+=2)
1053 {
1054 cc[i - 1 + (k + j*l1)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
1055 cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
1056 cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
1057 cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
1058 }
1059 }
1060 }
1061 }
1062 else
1063 {
1064 for (j=1; j<ipph; j++)
1065 {
1066 jc = ip - j;
1067 for (i=2; i<ido; i+=2)
1068 {
1069 for (k=0; k<l1; k++)
1070 {
1071 cc[i - 1 + (k + j*l1)*ido] =
1072 ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
1073 cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
1074 cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
1075 cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
1076 }
1077 }
1078 }
1079 }
1080 }
1081 else
1082 {
1083 for (ik=0; ik<idl1; ik++)
1084 cc[ik] = ch[ik];
1085 }
1086 for (j=1; j<ipph; j++)
1087 {
1088 jc = ip - j;
1089 for (k=0; k<l1; k++)
1090 {
1091 cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido] + ch[(k + jc*l1)*ido];
1092 cc[(k + jc*l1)*ido] = ch[(k + jc*l1)*ido] - ch[(k + j*l1)*ido];
1093 }
1094 }
1095
1096 ar1 = 1;
1097 ai1 = 0;
1098 for (l=1; l<ipph; l++)
1099 {
1100 lc = ip - l;
1101 ar1h = dcp*ar1 - dsp*ai1;
1102 ai1 = dcp*ai1 + dsp*ar1;
1103 ar1 = ar1h;
1104 for (ik=0; ik<idl1; ik++)
1105 {
1106 ch[ik + l*idl1] = cc[ik] + ar1*cc[ik + idl1];
1107 ch[ik + lc*idl1] = ai1*cc[ik + (ip-1)*idl1];
1108 }
1109 dc2 = ar1;
1110 ds2 = ai1;
1111 ar2 = ar1;
1112 ai2 = ai1;
1113 for (j=2; j<ipph; j++)
1114 {
1115 jc = ip - j;
1116 ar2h = dc2*ar2 - ds2*ai2;
1117 ai2 = dc2*ai2 + ds2*ar2;
1118 ar2 = ar2h;
1119 for (ik=0; ik<idl1; ik++)
1120 {
1121 ch[ik + l*idl1] += ar2*cc[ik + j*idl1];
1122 ch[ik + lc*idl1] += ai2*cc[ik + jc*idl1];
1123 }
1124 }
1125 }
1126 for (j=1; j<ipph; j++)
1127 for (ik=0; ik<idl1; ik++)
1128 ch[ik] += cc[ik + j*idl1];
1129
1130 if (ido >= l1)
1131 {
1132 for (k=0; k<l1; k++)
1133 {
1134 for (i=0; i<ido; i++)
1135 {
1136 cc[i + k*ip*ido] = ch[i + k*ido];
1137 }
1138 }
1139 }
1140 else
1141 {
1142 for (i=0; i<ido; i++)
1143 {
1144 for (k=0; k<l1; k++)
1145 {
1146 cc[i + k*ip*ido] = ch[i + k*ido];
1147 }
1148 }
1149 }
1150 for (j=1; j<ipph; j++)
1151 {
1152 jc = ip - j;
1153 j2 = 2*j;
1154 for (k=0; k<l1; k++)
1155 {
1156 cc[ido-1 + (j2 - 1 + k*ip)*ido] = ch[(k + j*l1)*ido];
1157 cc[(j2 + k*ip)*ido] = ch[(k + jc*l1)*ido];
1158 }
1159 }
1160 if (ido == 1) return;
1161 if (nbd >= l1)
1162 {
1163 for (j=1; j<ipph; j++)
1164 {
1165 jc = ip - j;
1166 j2 = 2*j;
1167 for (k=0; k<l1; k++)
1168 {
1169 for (i=2; i<ido; i+=2)
1170 {
1171 ic = ido - i;
1172 cc[i - 1 + (j2 + k*ip)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
1173 cc[ic - 1 + (j2 - 1 + k*ip)*ido] = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
1174 cc[i + (j2 + k*ip)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
1175 cc[ic + (j2 - 1 + k*ip)*ido] = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
1176 }
1177 }
1178 }
1179 }
1180 else
1181 {
1182 for (j=1; j<ipph; j++)
1183 {
1184 jc = ip - j;
1185 j2 = 2*j;
1186 for (i=2; i<ido; i+=2)
1187 {
1188 ic = ido - i;
1189 for (k=0; k<l1; k++)
1190 {
1191 cc[i - 1 + (j2 + k*ip)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
1192 cc[ic - 1 + (j2 - 1 + k*ip)*ido] = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
1193 cc[i + (j2 + k*ip)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
1194 cc[ic + (j2 - 1 + k*ip)*ido] = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
1195 }
1196 }
1197 }
1198 }
1199 }
1200
1201
1202 static void
1203 fftpack_radbg(int ido,
1204 int ip,
1205 int l1,
1206 int idl1,
1207 real cc[],
1208 real ch[],
1209 real wa[])
1210 {
1211 const real twopi = 6.28318530717959;
1212 int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
1213 real dc2, ai1, ai2, ar1, ar2, ds2;
1214 int nbd;
1215 real dcp, arg, dsp, ar1h, ar2h;
1216 arg = twopi / ip;
1217 dcp = cos(arg);
1218 dsp = sin(arg);
1219 nbd = (ido - 1) / 2;
1220 ipph = (ip + 1) / 2;
1221
1222 if (ido >= l1)
1223 {
1224 for (k=0; k<l1; k++)
1225 {
1226 for (i=0; i<ido; i++)
1227 {
1228 ch[i + k*ido] = cc[i + k*ip*ido];
1229 }
1230 }
1231 }
1232 else
1233 {
1234 for (i=0; i<ido; i++)
1235 {
1236 for (k=0; k<l1; k++)
1237 {
1238 ch[i + k*ido] = cc[i + k*ip*ido];
1239 }
1240 }
1241 }
1242 for (j=1; j<ipph; j++)
1243 {
1244 jc = ip - j;
1245 j2 = 2*j;
1246 for (k=0; k<l1; k++)
1247 {
1248 ch[(k + j*l1)*ido] = cc[ido-1 + (j2 - 1 + k*ip)*ido] + cc[ido-1 + (j2 - 1 + k*ip)*ido];
1249 ch[(k + jc*l1)*ido] = cc[(j2 + k*ip)*ido] + cc[(j2 + k*ip)*ido];
1250 }
1251 }
1252
1253 if (ido != 1)
1254 {
1255 if (nbd >= l1)
1256 {
1257 for (j=1; j<ipph; j++)
1258 {
1259 jc = ip - j;
1260 for (k=0; k<l1; k++)
1261 {
1262 for (i=2; i<ido; i+=2)
1263 {
1264 ic = ido - i;
1265 ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (2*j + k*ip)*ido] + cc[ic - 1 + (2*j - 1 + k*ip)*ido];
1266 ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (2*j + k*ip)*ido] - cc[ic - 1 + (2*j - 1 + k*ip)*ido];
1267 ch[i + (k + j*l1)*ido] = cc[i + (2*j + k*ip)*ido] - cc[ic + (2*j - 1 + k*ip)*ido];
1268 ch[i + (k + jc*l1)*ido] = cc[i + (2*j + k*ip)*ido] + cc[ic + (2*j - 1 + k*ip)*ido];
1269 }
1270 }
1271 }
1272 }
1273 else
1274 {
1275 for (j=1; j<ipph; j++)
1276 {
1277 jc = ip - j;
1278 for (i=2; i<ido; i+=2)
1279 {
1280 ic = ido - i;
1281 for (k=0; k<l1; k++)
1282 {
1283 ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (2*j + k*ip)*ido] + cc[ic - 1 + (2*j - 1 + k*ip)*ido];
1284 ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (2*j + k*ip)*ido] - cc[ic - 1 + (2*j - 1 + k*ip)*ido];
1285 ch[i + (k + j*l1)*ido] = cc[i + (2*j + k*ip)*ido] - cc[ic + (2*j - 1 + k*ip)*ido];
1286 ch[i + (k + jc*l1)*ido] = cc[i + (2*j + k*ip)*ido] + cc[ic + (2*j - 1 + k*ip)*ido];
1287 }
1288 }
1289 }
1290 }
1291 }
1292
1293 ar1 = 1;
1294 ai1 = 0;
1295 for (l=1; l<ipph; l++)
1296 {
1297 lc = ip - l;
1298 ar1h = dcp*ar1 - dsp*ai1;
1299 ai1 = dcp*ai1 + dsp*ar1;
1300 ar1 = ar1h;
1301 for (ik=0; ik<idl1; ik++)
1302 {
1303 cc[ik + l*idl1] = ch[ik] + ar1*ch[ik + idl1];
1304 cc[ik + lc*idl1] = ai1*ch[ik + (ip-1)*idl1];
1305 }
1306 dc2 = ar1;
1307 ds2 = ai1;
1308 ar2 = ar1;
1309 ai2 = ai1;
1310 for (j=2; j<ipph; j++)
1311 {
1312 jc = ip - j;
1313 ar2h = dc2*ar2 - ds2*ai2;
1314 ai2 = dc2*ai2 + ds2*ar2;
1315 ar2 = ar2h;
1316 for (ik=0; ik<idl1; ik++)
1317 {
1318 cc[ik + l*idl1] += ar2*ch[ik + j*idl1];
1319 cc[ik + lc*idl1] += ai2*ch[ik + jc*idl1];
1320 }
1321 }
1322 }
1323 for (j=1; j<ipph; j++)
1324 {
1325 for (ik=0; ik<idl1; ik++)
1326 {
1327 ch[ik] += ch[ik + j*idl1];
1328 }
1329 }
1330 for (j=1; j<ipph; j++)
1331 {
1332 jc = ip - j;
1333 for (k=0; k<l1; k++)
1334 {
1335 ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido] - cc[(k + jc*l1)*ido];
1336 ch[(k + jc*l1)*ido] = cc[(k + j*l1)*ido] + cc[(k + jc*l1)*ido];
1337 }
1338 }
1339
1340 if (ido == 1) return;
1341 if (nbd >= l1)
1342 {
1343 for (j=1; j<ipph; j++)
1344 {
1345 jc = ip - j;
1346 for (k=0; k<l1; k++)
1347 {
1348 for (i=2; i<ido; i+=2)
1349 {
1350 ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
1351 ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] + cc[i + (k + jc*l1)*ido];
1352 ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
1353 ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
1354 }
1355 }
1356 }
1357 }
1358 else
1359 {
1360 for (j=1; j<ipph; j++)
1361 {
1362 jc = ip - j;
1363 for (i=2; i<ido; i+=2)
1364 {
1365 for (k=0; k<l1; k++)
1366 {
1367 ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
1368 ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j *l1)*ido] + cc[i + (k + jc*l1)*ido];
1369 ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
1370 ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
1371 }
1372 }
1373 }
1374 }
1375 for (ik=0; ik<idl1; ik++)
1376 {
1377 cc[ik] = ch[ik];
1378 }
1379 for (j=1; j<ip; j++)
1380 for (k=0; k<l1; k++)
1381 cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido];
1382
1383 if (nbd <= l1)
1384 {
1385 is = -ido;
1386 for (j=1; j<ip; j++)
1387 {
1388 is += ido;
1389 idij = is-1;
1390 for (i=2; i<ido; i+=2)
1391 {
1392 idij += 2;
1393 for (k=0; k<l1; k++)
1394 {
1395 cc[i - 1 + (k + j*l1)*ido] = wa[idij - 1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*ch[i + (k + j*l1)*ido];
1396 cc[i + (k + j*l1)*ido] = wa[idij - 1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
1397 }
1398 }
1399 }
1400 }
1401 else
1402 {
1403 is = -ido;
1404 for (j=1; j<ip; j++)
1405 {
1406 is += ido;
1407 for (k=0; k<l1; k++)
1408 {
1409 idij = is;
1410 for (i=2; i<ido; i+=2)
1411 {
1412 idij += 2;
1413 cc[i - 1 + (k + j*l1)*ido] = wa[idij-1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*ch[i + (k + j*l1)*ido];
1414 cc[i + (k + j*l1)*ido] = wa[idij-1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
1415 }
1416 }
1417 }
1418 }
1419 }
1420
1421
1422
1423 static void
1424 fftpack_cfftf1(int n,
1425 real c[],
1426 real ch[],
1427 real wa[],
1428 int ifac[15],
1429 int isign)
1430 {
1431 int idot, i;
1432 int k1, l1, l2;
1433 int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
1434 real *cinput, *coutput;
1435 nf = ifac[1];
1436 na = 0;
1437 l1 = 1;
1438 iw = 0;
1439
1440 for (k1=2; k1<=nf+1; k1++)
1441 {
1442 ip = ifac[k1];
1443 l2 = ip*l1;
1444 ido = n / l2;
1445 idot = ido + ido;
1446 idl1 = idot*l1;
1447 if (na)
1448 {
1449 cinput = ch;
1450 coutput = c;
1451 }
1452 else
1453 {
1454 cinput = c;
1455 coutput = ch;
1456 }
1457 switch (ip)
1458 {
1459 case 4:
1460 ix2 = iw + idot;
1461 ix3 = ix2 + idot;
1462 fftpack_passf4(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], isign);
1463 na = !na;
1464 break;
1465 case 2:
1466 fftpack_passf2(idot, l1, cinput, coutput, &wa[iw], isign);
1467 na = !na;
1468 break;
1469 case 3:
1470 ix2 = iw + idot;
1471 fftpack_passf3(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], isign);
1472 na = !na;
1473 break;
1474 case 5:
1475 ix2 = iw + idot;
1476 ix3 = ix2 + idot;
1477 ix4 = ix3 + idot;
1478 fftpack_passf5(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
1479 na = !na;
1480 break;
1481 default:
1482 fftpack_passf(&nac, idot, ip, l1, idl1, cinput, coutput, &wa[iw], isign);
1483 if (nac != 0) na = !na;
1484 }
1485 l1 = l2;
1486 iw += (ip - 1)*idot;
1487 }
1488 if (na == 0)
1489 return;
1490 for (i=0; i<2*n; i++)
1491 c[i] = ch[i];
1492 }
1493
1494
1495 void
1496 fftpack_cfftf(int n,
1497 real c[],
1498 real wsave[])
1499 {
1500 int iw1, iw2;
1501
1502 if (n == 1)
1503 return;
1504 iw1 = 2*n;
1505 iw2 = iw1 + 2*n;
1506 fftpack_cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), -1);
1507 }
1508
1509
1510 void
1511 fftpack_cfftb(int n,
1512 real c[],
1513 real wsave[])
1514 {
1515 int iw1, iw2;
1516
1517 if (n == 1)
1518 return;
1519 iw1 = 2*n;
1520 iw2 = iw1 + 2*n;
1521 fftpack_cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), +1);
1522 }
1523
1524
1525 static void
1526 fftpack_factorize(int n,
1527 int ifac[15])
1528 {
1529 static const int ntryh[4] = { 3,4,2,5 };
1530 int ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr;
1531
1532 startloop:
1533 if (j < 4)
1534 ntry = ntryh[j];
1535 else
1536 ntry+= 2;
1537 j++;
1538 do
1539 {
1540 nq = nl / ntry;
1541 nr = nl - ntry*nq;
1542 if (nr != 0) goto startloop;
1543 nf++;
1544 ifac[nf + 1] = ntry;
1545 nl = nq;
1546 if (ntry == 2 && nf != 1)
1547 {
1548 for (i=2; i<=nf; i++)
1549 {
1550 ib = nf - i + 2;
1551 ifac[ib + 1] = ifac[ib];
1552 }
1553 ifac[2] = 2;
1554 }
1555 }
1556 while (nl != 1);
1557 ifac[0] = n;
1558 ifac[1] = nf;
1559 }
1560
1561
1562 static void
1563 fftpack_cffti1(int n,
1564 real wa[],
1565 int ifac[15])
1566 {
1567 const real twopi = 6.28318530717959;
1568 real arg, argh, argld, fi;
1569 int idot, i, j;
1570 int i1, k1, l1, l2;
1571 int ld, ii, nf, ip;
1572 int ido, ipm;
1573
1574 fftpack_factorize(n,ifac);
1575 nf = ifac[1];
1576 argh = twopi/(real)n;
1577 i = 1;
1578 l1 = 1;
1579 for (k1=1; k1<=nf; k1++)
1580 {
1581 ip = ifac[k1+1];
1582 ld = 0;
1583 l2 = l1*ip;
1584 ido = n / l2;
1585 idot = ido + ido + 2;
1586 ipm = ip - 1;
1587 for (j=1; j<=ipm; j++)
1588 {
1589 i1 = i;
1590 wa[i-1] = 1;
1591 wa[i] = 0;
1592 ld += l1;
1593 fi = 0;
1594 argld = ld*argh;
1595 for (ii=4; ii<=idot; ii+=2)
1596 {
1597 i+= 2;
1598 fi+= 1;
1599 arg = fi*argld;
1600 wa[i-1] = cos(arg);
1601 wa[i] = sin(arg);
1602 }
1603 if (ip > 5)
1604 {
1605 wa[i1-1] = wa[i-1];
1606 wa[i1] = wa[i];
1607 }
1608 }
1609 l1 = l2;
1610 }
1611 }
1612
1613
1614
1615
1616 static void
1617 fftpack_rfftf1(int n,
1618 real c[],
1619 real ch[],
1620 real wa[],
1621 int ifac[15])
1622 {
1623 int i;
1624 int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
1625 real *cinput, *coutput;
1626 nf = ifac[1];
1627 na = 1;
1628 l2 = n;
1629 iw = n-1;
1630 for (k1 = 1; k1 <= nf; ++k1)
1631 {
1632 kh = nf - k1;
1633 ip = ifac[kh + 2];
1634 l1 = l2 / ip;
1635 ido = n / l2;
1636 idl1 = ido*l1;
1637 iw -= (ip - 1)*ido;
1638 na = !na;
1639 if (na)
1640 {
1641 cinput = ch;
1642 coutput = c;
1643 }
1644 else
1645 {
1646 cinput = c;
1647 coutput = ch;
1648 }
1649 switch (ip)
1650 {
1651 case 4:
1652 ix2 = iw + ido;
1653 ix3 = ix2 + ido;
1654 fftpack_radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
1655 break;
1656 case 2:
1657 fftpack_radf2(ido, l1, cinput, coutput, &wa[iw]);
1658 break;
1659 case 3:
1660 ix2 = iw + ido;
1661 fftpack_radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
1662 break;
1663 case 5:
1664 ix2 = iw + ido;
1665 ix3 = ix2 + ido;
1666 ix4 = ix3 + ido;
1667 fftpack_radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
1668 break;
1669 default:
1670 if (ido == 1)
1671 na = !na;
1672 if (na == 0)
1673 {
1674 fftpack_radfg(ido, ip, l1, idl1, c, ch, &wa[iw]);
1675 na = 1;
1676 }
1677 else
1678 {
1679 fftpack_radfg(ido, ip, l1, idl1, ch, c, &wa[iw]);
1680 na = 0;
1681 }
1682 }
1683 l2 = l1;
1684 }
1685 if (na == 1)
1686 return;
1687 for (i = 0; i < n; i++)
1688 c[i] = ch[i];
1689 }
1690
1691
1692 static void
1693 fftpack_rfftb1(int n,
1694 real c[],
1695 real ch[],
1696 real wa[],
1697 int ifac[15])
1698 {
1699 int i;
1700 int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
1701 real *cinput, *coutput;
1702 nf = ifac[1];
1703 na = 0;
1704 l1 = 1;
1705 iw = 0;
1706
1707 for (k1=1; k1<=nf; k1++)
1708 {
1709 ip = ifac[k1 + 1];
1710 l2 = ip*l1;
1711 ido = n / l2;
1712 idl1 = ido*l1;
1713 if (na)
1714 {
1715 cinput = ch;
1716 coutput = c;
1717 }
1718 else
1719 {
1720 cinput = c;
1721 coutput = ch;
1722 }
1723 switch (ip)
1724 {
1725 case 4:
1726 ix2 = iw + ido;
1727 ix3 = ix2 + ido;
1728 fftpack_radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
1729 na = !na;
1730 break;
1731 case 2:
1732 fftpack_radb2(ido, l1, cinput, coutput, &wa[iw]);
1733 na = !na;
1734 break;
1735 case 3:
1736 ix2 = iw + ido;
1737 fftpack_radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
1738 na = !na;
1739 break;
1740 case 5:
1741 ix2 = iw + ido;
1742 ix3 = ix2 + ido;
1743 ix4 = ix3 + ido;
1744 fftpack_radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
1745 na = !na;
1746 break;
1747 default:
1748 fftpack_radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]);
1749 if (ido == 1) na = !na;
1750 }
1751 l1 = l2;
1752 iw += (ip - 1)*ido;
1753 }
1754 if (na == 0)
1755 return;
1756 for (i=0; i<n; i++)
1757 c[i] = ch[i];
1758 }
1759
1760
1761
1762
1763 static void
1764 fftpack_rffti1(int n,
1765 real wa[],
1766 int ifac[15])
1767 {
1768 const real twopi = 6.28318530717959;
1769 real arg, argh, argld, fi;
1770 int i, j;
1771 int k1, l1, l2;
1772 int ld, ii, nf, ip, is;
1773 int ido, ipm, nfm1;
1774 fftpack_factorize(n,ifac);
1775 nf = ifac[1];
1776 argh = twopi / n;
1777 is = 0;
1778 nfm1 = nf - 1;
1779 l1 = 1;
1780 if (nfm1 == 0) return;
1781 for (k1 = 1; k1 <= nfm1; k1++)
1782 {
1783 ip = ifac[k1 + 1];
1784 ld = 0;
1785 l2 = l1*ip;
1786 ido = n / l2;
1787 ipm = ip - 1;
1788 for (j = 1; j <= ipm; ++j)
1789 {
1790 ld += l1;
1791 i = is;
1792 argld = (real) ld*argh;
1793 fi = 0;
1794 for (ii = 3; ii <= ido; ii += 2)
1795 {
1796 i += 2;
1797 fi += 1;
1798 arg = fi*argld;
1799 wa[i - 2] = cos(arg);
1800 wa[i - 1] = sin(arg);
1801 }
1802 is += ido;
1803 }
1804 l1 = l2;
1805 }
1806 }
1807
1808
1809
1810
1811 /* End of fftpack - begin GROMACS code */
1812
1813
1814 int
1815 gmx_fft_init_1d(gmx_fft_t * pfft,
1816 int nx,
1817 int flags)
1818 {
1819 gmx_fft_t fft;
1820
1821 if(pfft==NULL)
1822 {
1823 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
1824 return EINVAL;
1825 }
1826 *pfft = NULL;
1827
1828 if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
1829 {
1830 return ENOMEM;
1831 }
1832
1833 fft->next = NULL;
1834 fft->n = nx;
1835
1836 /* Need 4*n storage for 1D complex FFT */
1837 if( (fft->work = (real *)malloc(sizeof(real)*(4*nx))) == NULL)
1838 {
1839 free(fft);
1840 return ENOMEM;
1841 }
1842
1843 if(fft->n>1)
1844 fftpack_cffti1(nx,fft->work,fft->ifac);
1845
1846 *pfft = fft;
1847 return 0;
1848 };
1849
1850
1851
1852 int
1853 gmx_fft_init_1d_real(gmx_fft_t * pfft,
1854 int nx,
1855 int flags)
1856 {
1857 gmx_fft_t fft;
1858
1859 if(pfft==NULL)
1860 {
1861 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
1862 return EINVAL;
1863 }
1864 *pfft = NULL;
1865
1866 if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
1867 {
1868 return ENOMEM;
1869 }
1870
1871 fft->next = NULL;
1872 fft->n = nx;
1873
1874 /* Need 2*n storage for 1D real FFT */
1875 if((fft->work = (real *)malloc(sizeof(real)*(2*nx)))==NULL)
1876 {
1877 free(fft);
1878 return ENOMEM;
1879 }
1880
1881 if(fft->n>1)
1882 fftpack_rffti1(nx,fft->work,fft->ifac);
1883
1884 *pfft = fft;
1885 return 0;
1886 }
1887
1888
1889
1890 int
1891 gmx_fft_init_2d(gmx_fft_t * pfft,
1892 int nx,
1893 int ny,
1894 int flags)
1895 {
1896 gmx_fft_t fft;
1897 int rc;
1898
1899 if(pfft==NULL)
1900 {
1901 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
1902 return EINVAL;
1903 }
1904 *pfft = NULL;
1905
1906 /* Create the X transform */
1907 if( (rc = gmx_fft_init_1d(&fft,nx,flags)) != 0)
1908 {
1909 return rc;
1910 }
1911
1912 /* Create Y transform as a link from X */
1913 if( (rc=gmx_fft_init_1d(&(fft->next),ny,flags)) != 0)
1914 {
1915 free(fft);
1916 return rc;
1917 }
1918
1919 *pfft = fft;
1920 return 0;
1921 };
1922
1923
1924 int
1925 gmx_fft_init_2d_real(gmx_fft_t * pfft,
1926 int nx,
1927 int ny,
1928 int flags)
1929 {
1930 gmx_fft_t fft;
1931 int nyc = (ny/2 + 1);
1932 int rc;
1933
1934 if(pfft==NULL)
1935 {
1936 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
1937 return EINVAL;
1938 }
1939 *pfft = NULL;
1940
1941 /* Create the X transform */
1942 if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
1943 {
1944 return ENOMEM;
1945 }
1946
1947 fft->n = nx;
1948
1949 /* Need 4*nx storage for 1D complex FFT, and another
1950 * 2*nx*nyc elements for complex-to-real storage in our high-level routine.
1951 */
1952 if( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*nyc))) == NULL)
1953 {
1954 free(fft);
1955 return ENOMEM;
1956 }
1957 fftpack_cffti1(nx,fft->work,fft->ifac);
1958
1959 /* Create real Y transform as a link from X */
1960 if( (rc=gmx_fft_init_1d_real(&(fft->next),ny,flags)) != 0)
1961 {
1962 free(fft);
1963 return rc;
1964 }
1965
1966 *pfft = fft;
1967 return 0;
1968 }
1969
1970
1971 int
1972 gmx_fft_init_3d(gmx_fft_t * pfft,
1973 int nx,
1974 int ny,
1975 int nz,
1976 int flags)
1977 {
1978 gmx_fft_t fft;
1979 int rc;
1980
1981 if(pfft==NULL)
1982 {
1983 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
1984 return EINVAL;
1985 }
1986 *pfft = NULL;
1987
1988 /* Create the X transform */
1989
1990 if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
1991 {
1992 return ENOMEM;
1993 }
1994
1995 fft->n = nx;
1996
1997 /* Need 4*nx storage for 1D complex FFT, and another
1998 * 2*nz elements for gmx_fft_transpose_2d_nelem() storage.
1999 */
2000 if( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nz))) == NULL)
2001 {
2002 free(fft);
2003 return ENOMEM;
2004 }
2005
2006 fftpack_cffti1(nx,fft->work,fft->ifac);
2007
2008
2009 /* Create 2D Y/Z transforms as a link from X */
2010 if( (rc=gmx_fft_init_2d(&(fft->next),ny,nz,flags)) != 0)
2011 {
2012 free(fft);
2013 return rc;
2014 }
2015
2016 *pfft = fft;
2017 return 0;
2018 };
2019
2020
2021 int
2022 gmx_fft_init_3d_real(gmx_fft_t * pfft,
2023 int nx,
2024 int ny,
2025 int nz,
2026 int flags)
2027 {
2028 gmx_fft_t fft;
2029 int nzc = (nz/2 + 1);
2030 int rc;
2031
2032 if(pfft==NULL)
2033 {
2034 gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
2035 return EINVAL;
2036 }
2037 *pfft = NULL;
2038
2039 /* Create the X transform */
2040 if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
2041 {
2042 return ENOMEM;
2043 }
2044
2045 fft->n = nx;
2046
2047 /* Need 4*nx storage for 1D complex FFT, another
2048 * 2*nx*ny*nzc elements to copy the entire 3D matrix when
2049 * doing out-of-place complex-to-real FFTs, and finally
2050 * 2*nzc elements for transpose work space.
2051 */
2052 if( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*ny*nzc+2*nzc))) == NULL)
2053 {
2054 free(fft);
2055 return ENOMEM;
2056 }
2057 fftpack_cffti1(nx,fft->work,fft->ifac);
2058
2059 /* Create 2D real Y/Z transform as a link from X */
2060 if( (rc=gmx_fft_init_2d_real(&(fft->next),ny,nz,flags)) != 0)
2061 {
2062 free(fft);
2063 return rc;
2064 }
2065
2066 *pfft = fft;
2067 return 0;
2068 }
2069
2070
2071 int
2072 gmx_fft_1d (gmx_fft_t fft,
2073 enum gmx_fft_direction dir,
2074 void * in_data,
2075 void * out_data)
2076 {
2077 int i,n;
2078 real * p1;
2079 real * p2;
2080
2081 n=fft->n;
2082
2083 if(n==1)
2084 {
2085 p1 = (real *)in_data;
2086 p2 = (real *)out_data;
2087 p2[0] = p1[0];
2088 p2[1] = p1[1];
2089 }
2090
2091 /* FFTPACK only does in-place transforms, so emulate out-of-place
2092 * by copying data to the output array first.
2093 */
2094 if( in_data != out_data )
2095 {
2096 p1 = (real *)in_data;
2097 p2 = (real *)out_data;
2098
2099 /* n complex = 2*n real elements */
2100 for(i=0;i<2*n;i++)
2101 {
2102 p2[i] = p1[i];
2103 }
2104 }
2105
2106 /* Elements 0 .. 2*n-1 in work are used for ffac values,
2107 * Elements 2*n .. 4*n-1 are internal FFTPACK work space.
2108 */
2109
2110 if(dir == GMX_FFT_FORWARD)
2111 {
2112 fftpack_cfftf1(n,(real *)out_data,fft->work+2*n,fft->work,fft->ifac, -1);
2113 }
2114 else if(dir == GMX_FFT_BACKWARD)
2115 {
2116 fftpack_cfftf1(n,(real *)out_data,fft->work+2*n,fft->work,fft->ifac, 1);
2117 }
2118 else
2119 {
2120 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
2121 return EINVAL;
2122 }
2123
2124 return 0;
2125 }
2126
2127
2128
2129 int
2130 gmx_fft_1d_real (gmx_fft_t fft,
2131 enum gmx_fft_direction dir,
2132 void * in_data,
2133 void * out_data)
2134 {
2135 int i,n;
2136 real * p1;
2137 real * p2;
2138
2139 n = fft->n;
2140
2141 if(n==1)
2142 {
2143 p1 = (real *)in_data;
2144 p2 = (real *)out_data;
2145 p2[0] = p1[0];
2146 if(dir == GMX_FFT_REAL_TO_COMPLEX)
2147 p2[1] = 0.0;
2148 }
2149
2150 if(dir == GMX_FFT_REAL_TO_COMPLEX)
2151 {
2152 /* FFTPACK only does in-place transforms, so emulate out-of-place
2153 * by copying data to the output array first. This works fine, since
2154 * the complex array must be larger than the real.
2155 */
2156 if( in_data != out_data )
2157 {
2158 p1 = (real *)in_data;
2159 p2 = (real *)out_data;
2160
2161 for(i=0;i<2*(n/2+1);i++)
2162 {
2163 p2[i] = p1[i];
2164 }
2165 }
2166
2167 /* Elements 0 .. n-1 in work are used for ffac values,
2168 * Elements n .. 2*n-1 are internal FFTPACK work space.
2169 */
2170 fftpack_rfftf1(n,(real *)out_data,fft->work+n,fft->work,fft->ifac);
2171
2172 /*
2173 * FFTPACK has a slightly more compact storage than we, time to
2174 * convert it: ove most of the array one step up to make room for
2175 * zero imaginary parts.
2176 */
2177 p2 = (real *)out_data;
2178 for(i=n-1;i>0;i--)
2179 {
2180 p2[i+1] = p2[i];
2181 }
2182 /* imaginary zero freq. */
2183 p2[1] = 0;
2184
2185 /* Is n even? */
2186 if( (n & 0x1) == 0 )
2187 {
2188 p2[n+1] = 0;
2189 }
2190
2191 }
2192 else if(dir == GMX_FFT_COMPLEX_TO_REAL)
2193 {
2194 /* FFTPACK only does in-place transforms, and we cannot just copy
2195 * input to output first here since our real array is smaller than
2196 * the complex one. However, since the FFTPACK complex storage format
2197 * is more compact than ours (2 reals) it will fit, so compact it
2198 * and copy on-the-fly to the output array.
2199 */
2200 p1 = (real *) in_data;
2201 p2 = (real *)out_data;
2202
2203 p2[0] = p1[0];
2204 for(i=1;i<n;i++)
2205 {
2206 p2[i] = p1[i+1];
2207 }
2208 fftpack_rfftb1(n,(real *)out_data,fft->work+n,fft->work,fft->ifac);
2209 }
2210 else
2211 {
2212 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
2213 return EINVAL;
2214 }
2215
2216 return 0;
2217 }
2218
2219
2220 int
2221 gmx_fft_2d (gmx_fft_t fft,
2222 enum gmx_fft_direction dir,
2223 void * in_data,
2224 void * out_data)
2225 {
2226 int i,nx,ny;
2227 t_complex * data;
2228
2229 nx = fft->n;
2230 ny = fft->next->n;
2231
2232 /* FFTPACK only does in-place transforms, so emulate out-of-place
2233 * by copying data to the output array first.
2234 * For 2D there is likely enough data to benefit from memcpy().
2235 */
2236 if( in_data != out_data )
2237 {
2238 memcpy(out_data,in_data,sizeof(t_complex)*nx*ny);
2239 }
2240
2241 /* Much easier to do pointer arithmetic when base has the correct type */
2242 data = (t_complex *)out_data;
2243
2244 /* y transforms */
2245 for(i=0;i<nx;i++)
2246 {
2247 gmx_fft_1d(fft->next,dir,data+i*ny,data+i*ny);
2248 }
2249
2250 /* Transpose in-place to get data in place for x transform now */
2251 gmx_fft_transpose_2d(data,data,nx,ny);
2252
2253 /* x transforms */
2254 for(i=0;i<ny;i++)
2255 {
2256 gmx_fft_1d(fft,dir,data+i*nx,data+i*nx);
2257 }
2258
2259 /* Transpose in-place to get data back in original order */
2260 gmx_fft_transpose_2d(data,data,ny,nx);
2261
2262 return 0;
2263 }
2264
2265
2266
2267 int
2268 gmx_fft_2d_real (gmx_fft_t fft,
2269 enum gmx_fft_direction dir,
2270 void * in_data,
2271 void * out_data)
2272 {
2273 int i,j,nx,ny,nyc;
2274 t_complex * data;
2275 real * work;
2276 real * p1;
2277 real * p2;
2278
2279 nx=fft->n;
2280 ny=fft->next->n;
2281 /* Number of complex elements in y direction */
2282 nyc=(ny/2+1);
2283
2284 work = fft->work+4*nx;
2285
2286 if(dir==GMX_FFT_REAL_TO_COMPLEX)
2287 {
2288 /* If we are doing an in-place transform the 2D array is already
2289 * properly padded by the user, and we are all set.
2290 *
2291 * For out-of-place there is no array padding, but FFTPACK only
2292 * does in-place FFTs internally, so we need to start by copying
2293 * data from the input to the padded (larger) output array.
2294 */
2295 if( in_data != out_data )
2296 {
2297 p1 = (real *)in_data;
2298 p2 = (real *)out_data;
2299
2300 for(i=0;i<nx;i++)
2301 {
2302 for(j=0;j<ny;j++)
2303 {
2304 p2[i*nyc*2+j] = p1[i*ny+j];
2305 }
2306 }
2307 }
2308 data = (t_complex *)out_data;
2309
2310 /* y real-to-complex FFTs */
2311 for(i=0;i<nx;i++)
2312 {
2313 gmx_fft_1d_real(fft->next,GMX_FFT_REAL_TO_COMPLEX,data+i*nyc,data+i*nyc);
2314 }
2315
2316 /* Transform to get X data in place */
2317 gmx_fft_transpose_2d(data,data,nx,nyc);
2318
2319 /* Complex-to-complex X FFTs */
2320 for(i=0;i<nyc;i++)
2321 {
2322 gmx_fft_1d(fft,GMX_FFT_FORWARD,data+i*nx,data+i*nx);
2323 }
2324
2325 /* Transpose back */
2326 gmx_fft_transpose_2d(data,data,nyc,nx);
2327
2328 }
2329 else if(dir==GMX_FFT_COMPLEX_TO_REAL)
2330 {
2331 /* An in-place complex-to-real transform is straightforward,
2332 * since the output array must be large enough for the padding to fit.
2333 *
2334 * For out-of-place complex-to-real transforms we cannot just copy
2335 * data to the output array, since it is smaller than the input.
2336 * In this case there's nothing to do but employing temporary work data,
2337 * starting at work+4*nx and using nx*nyc*2 elements.
2338 */
2339 if(in_data != out_data)
2340 {
2341 memcpy(work,in_data,sizeof(t_complex)*nx*nyc);
2342 data = (t_complex *)work;
2343 }
2344 else
2345 {
2346 /* in-place */
2347 data = (t_complex *)out_data;
2348 }
2349
2350 /* Transpose to get X arrays */
2351 gmx_fft_transpose_2d(data,data,nx,nyc);
2352
2353 /* Do X iFFTs */
2354 for(i=0;i<nyc;i++)
2355 {
2356 gmx_fft_1d(fft,GMX_FFT_BACKWARD,data+i*nx,data+i*nx);
2357 }
2358
2359 /* Transpose to get Y arrays */
2360 gmx_fft_transpose_2d(data,data,nyc,nx);
2361
2362 /* Do Y iFFTs */
2363 for(i=0;i<nx;i++)
2364 {
2365 gmx_fft_1d_real(fft->next,GMX_FFT_COMPLEX_TO_REAL,data+i*nyc,data+i*nyc);
2366 }
2367
2368 if( in_data != out_data )
2369 {
2370 /* Output (pointed to by data) is now in padded format.
2371 * Pack it into out_data if we were doing an out-of-place transform.
2372 */
2373 p1 = (real *)data;
2374 p2 = (real *)out_data;
2375
2376 for(i=0;i<nx;i++)
2377 {
2378 for(j=0;j<ny;j++)
2379 {
2380 p2[i*ny+j] = p1[i*nyc*2+j];
2381 }
2382 }
2383 }
2384 }
2385 else
2386 {
2387 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
2388 return EINVAL;
2389 }
2390
2391 return 0;
2392 }
2393
2394
2395
2396 int
2397 gmx_fft_3d (gmx_fft_t fft,
2398 enum gmx_fft_direction dir,
2399 void * in_data,
2400 void * out_data)
2401 {
2402 int i,nx,ny,nz,rc;
2403 t_complex * data;
2404 t_complex * work;
2405 nx=fft->n;
2406 ny=fft->next->n;
2407 nz=fft->next->next->n;
2408
2409 /* First 4*nx positions are FFTPACK workspace, then ours starts */
2410 work = (t_complex *)(fft->work+4*nx);
2411
2412 /* FFTPACK only does in-place transforms, so emulate out-of-place
2413 * by copying data to the output array first.
2414 * For 3D there is likely enough data to benefit from memcpy().
2415 */
2416 if( in_data != out_data )
2417 {
2418 memcpy(out_data,in_data,sizeof(t_complex)*nx*ny*nz);
2419 }
2420
2421 /* Much easier to do pointer arithmetic when base has the correct type */
2422 data = (t_complex *)out_data;
2423
2424 /* Perform z transforms */
2425 for(i=0;i<nx*ny;i++)
2426 gmx_fft_1d(fft->next->next,dir,data+i*nz,data+i*nz);
2427
2428 /* For each X slice, transpose the y & z dimensions inside the slice */
2429 for(i=0;i<nx;i++)
2430 {
2431 gmx_fft_transpose_2d(data+i*ny*nz,data+i*ny*nz,ny,nz);
2432 }
2433
2434 /* Array is now (nx,nz,ny) - perform y transforms */
2435 for(i=0;i<nx*nz;i++)
2436 {
2437 gmx_fft_1d(fft->next,dir,data+i*ny,data+i*ny);
2438 }
2439
2440 /* Transpose back to (nx,ny,nz) */
2441 for(i=0;i<nx;i++)
2442 {
2443 gmx_fft_transpose_2d(data+i*ny*nz,data+i*ny*nz,nz,ny);
2444 }
2445
2446 /* Transpose entire x & y slices to go from
2447 * (nx,ny,nz) to (ny,nx,nz).
2448 * Use work data elements 4*n .. 4*n+2*nz-1.
2449 */
2450 rc=gmx_fft_transpose_2d_nelem(data,data,nx,ny,nz,work);
2451 if( rc != 0)
2452 {
2453 gmx_fatal(FARGS,"Cannot transpose X & Y/Z in gmx_fft_3d().");
2454 return rc;
2455 }
2456
2457 /* Then go from (ny,nx,nz) to (ny,nz,nx) */
2458 for(i=0;i<ny;i++)
2459 {
2460 gmx_fft_transpose_2d(data+i*nx*nz,data+i*nx*nz,nx,nz);
2461 }
2462
2463 /* Perform x transforms */
2464 for(i=0;i<ny*nz;i++)
2465 {
2466 gmx_fft_1d(fft,dir,data+i*nx,data+i*nx);
2467 }
2468
2469 /* Transpose back from (ny,nz,nx) to (ny,nx,nz) */
2470 for(i=0;i<ny;i++)
2471 {
2472 gmx_fft_transpose_2d(data+i*nz*nx,data+i*nz*nx,nz,nx);
2473 }
2474
2475 /* Transpose from (ny,nx,nz) to (nx,ny,nz)
2476 * Use work data elements 4*n .. 4*n+2*nz-1.
2477 */
2478 rc = gmx_fft_transpose_2d_nelem(data,data,ny,nx,nz,work);
2479 if( rc != 0)
2480 {
2481 gmx_fatal(FARGS,"Cannot transpose Y/Z & X in gmx_fft_3d().");
2482 return rc;
2483 }
2484
2485 return 0;
2486 }
2487
2488
2489 int
2490 gmx_fft_3d_real (gmx_fft_t fft,
2491 enum gmx_fft_direction dir,
2492 void * in_data,
2493 void * out_data)
2494 {
2495 int i,j,k;
2496 int nx,ny,nz,nzc,rc;
2497 t_complex * data;
2498 t_complex * work_transp;
2499 t_complex * work_c2r;
2500 real * p1;
2501 real * p2;
2502
2503 nx=fft->n;
2504 ny=fft->next->n;
2505 nz=fft->next->next->n;
2506 nzc=(nz/2+1);
2507
2508
2509 /* First 4*nx positions are FFTPACK workspace, then ours starts.
2510 * We have 2*nx*ny*nzc elements for temp complex-to-real storage when
2511 * doing out-of-place transforms, and another 2*nzc for transpose data.
2512 */
2513 work_c2r = (t_complex *)(fft->work+4*nx);
2514 work_transp = (t_complex *)(fft->work+4*nx+2*nx*ny*nzc);
2515
2516 /* Much easier to do pointer arithmetic when base has the correct type */
2517 data = (t_complex *)out_data;
2518
2519 if(dir==GMX_FFT_REAL_TO_COMPLEX)
2520 {
2521 /* FFTPACK only does in-place transforms, so emulate out-of-place
2522 * by copying data to the output array first. This is guaranteed to
2523 * work for real-to-complex since complex data is larger than the real.
2524 * For 3D there is likely enough data to benefit from memcpy().
2525 */
2526 if( in_data != out_data )
2527 {
2528 p1 = (real *)in_data;
2529 p2 = (real *)out_data;
2530
2531 for(i=0;i<nx;i++)
2532 {
2533 for(j=0;j<ny;j++)
2534 {
2535 for(k=0;k<nz;k++)
2536 {
2537 p2[(i*ny+j)*2*nzc+k] = p1[(i*ny+j)*nz+k];
2538 }
2539 }
2540 }
2541 }
2542 data = (t_complex *)out_data;
2543
2544 /* Transform the Y/Z slices real-to-complex */
2545 for(i=0;i<nx;i++)
2546 {
2547 gmx_fft_2d_real(fft->next,dir,data+i*ny*nzc,data+i*ny*nzc);
2548 }
2549
2550 /* Transpose x & y slices to go from
2551 * (nx,ny,nzc) to (ny,nx,nzc).
2552 */
2553 rc=gmx_fft_transpose_2d_nelem(data,data,nx,ny,nzc,work_transp);
2554 if( rc != 0)
2555 {
2556 gmx_fatal(FARGS,"Cannot transpose X & Y/Z gmx_fft_3d_real().");
2557 return rc;
2558 }
2559
2560 /* Then transpose from (ny,nx,nzc) to (ny,nzc,nx) */
2561 for(i=0;i<ny;i++)
2562 {
2563 gmx_fft_transpose_2d(data+i*nx*nzc,data+i*nx*nzc,nx,nzc);
2564 }
2565
2566 /* Perform x transforms */
2567 for(i=0;i<ny*nzc;i++)
2568 {
2569 gmx_fft_1d(fft,GMX_FFT_FORWARD,data+i*nx,data+i*nx);
2570 }
2571
2572 /* Transpose from (ny,nzc,nx) back to (ny,nx,nzc) */
2573 for(i=0;i<ny;i++)
2574 {
2575 gmx_fft_transpose_2d(data+i*nzc*nx,data+i*nzc*nx,nzc,nx);
2576 }
2577
2578 /* Transpose back from (ny,nx,nzc) to (nx,ny,nz) */
2579 rc=gmx_fft_transpose_2d_nelem(data,data,ny,nx,nzc,work_transp);
2580 if( rc != 0)
2581 {
2582 gmx_fatal(FARGS,"Cannot transpose Y/Z & X in gmx_fft_3d_real().");
2583 return rc;
2584 }
2585
2586 }
2587 else if(dir==GMX_FFT_COMPLEX_TO_REAL)
2588 {
2589 /* An in-place complex-to-real transform is straightforward,
2590 * since the output array must be large enough for the padding to fit.
2591 *
2592 * For out-of-place complex-to-real transforms we cannot just copy
2593 * data to the output array, since it is smaller than the input.
2594 * In this case there's nothing to do but employing temporary work data.
2595 */
2596 if(in_data != out_data)
2597 {
2598 memcpy(work_c2r,in_data,sizeof(t_complex)*nx*ny*nzc);
2599 data = (t_complex *)work_c2r;
2600 }
2601 else
2602 {
2603 /* in-place */
2604 data = (t_complex *)out_data;
2605 }
2606
2607 /* Transpose x & y slices to go from
2608 * (nx,ny,nz) to (ny,nx,nz).
2609 */
2610 gmx_fft_transpose_2d_nelem(data,data,nx,ny,nzc,work_transp);
2611
2612 /* Then go from (ny,nx,nzc) to (ny,nzc,nx) */
2613 for(i=0;i<ny;i++)
2614 {
2615 gmx_fft_transpose_2d(data+i*nx*nzc,data+i*nx*nzc,nx,nzc);
2616 }
2617
2618
2619 /* Perform x transforms */
2620 for(i=0;i<ny*nzc;i++)
2621 {
2622 gmx_fft_1d(fft,GMX_FFT_BACKWARD,data+i*nx,data+i*nx);
2623 }
2624
2625 /* Transpose back from (ny,nzc,nx) to (ny,nx,nzc) */
2626 for(i=0;i<ny;i++)
2627 {
2628 gmx_fft_transpose_2d(data+i*nzc*nx,data+i*nzc*nx,nzc,nx);
2629 }
2630
2631 /* Transpose back from (ny,nx,nzc) to (nx,ny,nz) */
2632 gmx_fft_transpose_2d_nelem(data,data,ny,nx,nzc,work_transp);
2633
2634
2635 /* Do 2D complex-to-real */
2636 for(i=0;i<nx;i++)
2637 {
2638 gmx_fft_2d_real(fft->next,dir,data+i*ny*nzc,data+i*ny*nzc);
2639 }
2640
2641 if( in_data != out_data )
2642 {
2643 /* Output (pointed to by data) is now in padded format.
2644 * Pack it into out_data if we were doing an out-of-place transform.
2645 */
2646 p1 = (real *)data;
2647 p2 = (real *)out_data;
2648
2649 for(i=0;i<nx;i++)
2650 {
2651 for(j=0;j<ny;j++)
2652 {
2653 for(k=0;k<nz;k++)
2654 {
2655 p2[(i*ny+j)*nz+k] = p1[(i*ny+j)*nzc*2+k];
2656 }
2657 }
2658 }
2659 }
2660
2661 }
2662 else
2663 {
2664 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
2665 return EINVAL;
2666 }
2667
2668 return 0;
2669 }
2670
2671
2672
2673
2674 void
2675 gmx_fft_destroy(gmx_fft_t fft)
2676 {
2677 if(fft != NULL)
2678 {
2679 free(fft->work);
2680 if(fft->next != NULL)
2681 gmx_fft_destroy(fft->next);
2682 free(fft);
2683 }
2684 }
2685 #else
2686 int
2687 gmx_fft_fftpack_empty;
2688 #endif /* GMX_FFT_FFTPACK */