NAMD
pub3dfft.C
Go to the documentation of this file.
1 
7 //
8 // XXX static and global variables are unsafe for shared memory builds.
9 // This FFT module cannot safely be used by multithreaded builds of NAMD.
10 //
11 
12 #include <math.h>
13 #include "pub3dfft.h"
14 #ifndef M_PI
15 #define M_PI 3.14159265358979323846
16 #endif
17 
18 /* Subroutine */ int passb(int *nac, int *ido, int *ip, int *
19  l1, int *idl1, double *cc, double *c1, double *c2,
20  double *ch, double *ch2, double *wa)
21 {
22  /* System generated locals */
23  int ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1,
24  c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset,
25  i_1, i_2, i_3;
26 
27  /* Local variables */
28  static int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj,
29  idl, inc, idp;
30  static double wai, war;
31  static int ipp2;
32 
33  /* Parameter adjustments */
34  cc_dim1 = *ido;
35  cc_dim2 = *ip;
36  cc_offset = cc_dim1 * (cc_dim2 + 1) + 1;
37  cc -= cc_offset;
38  c1_dim1 = *ido;
39  c1_dim2 = *l1;
40  c1_offset = c1_dim1 * (c1_dim2 + 1) + 1;
41  c1 -= c1_offset;
42  c2_dim1 = *idl1;
43  c2_offset = c2_dim1 + 1;
44  c2 -= c2_offset;
45  ch_dim1 = *ido;
46  ch_dim2 = *l1;
47  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
48  ch -= ch_offset;
49  ch2_dim1 = *idl1;
50  ch2_offset = ch2_dim1 + 1;
51  ch2 -= ch2_offset;
52  --wa;
53 
54  /* Function Body */
55  idot = *ido / 2;
56  nt = *ip * *idl1;
57  ipp2 = *ip + 2;
58  ipph = (*ip + 1) / 2;
59  idp = *ip * *ido;
60 
61  if (*ido < *l1) {
62  goto L106;
63  }
64  i_1 = ipph;
65  for (j = 2; j <= i_1; ++j) {
66  jc = ipp2 - j;
67  i_2 = *l1;
68  for (k = 1; k <= i_2; ++k) {
69  i_3 = *ido;
70  for (i = 1; i <= i_3; ++i) {
71  ch[i + (k + j * ch_dim2) * ch_dim1] = cc[i + (j + k * cc_dim2)
72  * cc_dim1] + cc[i + (jc + k * cc_dim2) * cc_dim1];
73  ch[i + (k + jc * ch_dim2) * ch_dim1] = cc[i + (j + k *
74  cc_dim2) * cc_dim1] - cc[i + (jc + k * cc_dim2) *
75  cc_dim1];
76 /* L101: */
77  }
78 /* L102: */
79  }
80 /* L103: */
81  }
82  i_1 = *l1;
83  for (k = 1; k <= i_1; ++k) {
84  i_2 = *ido;
85  for (i = 1; i <= i_2; ++i) {
86  ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) *
87  cc_dim1];
88 /* L104: */
89  }
90 /* L105: */
91  }
92  goto L112;
93 L106:
94  i_1 = ipph;
95  for (j = 2; j <= i_1; ++j) {
96  jc = ipp2 - j;
97  i_2 = *ido;
98  for (i = 1; i <= i_2; ++i) {
99  i_3 = *l1;
100  for (k = 1; k <= i_3; ++k) {
101  ch[i + (k + j * ch_dim2) * ch_dim1] = cc[i + (j + k * cc_dim2)
102  * cc_dim1] + cc[i + (jc + k * cc_dim2) * cc_dim1];
103  ch[i + (k + jc * ch_dim2) * ch_dim1] = cc[i + (j + k *
104  cc_dim2) * cc_dim1] - cc[i + (jc + k * cc_dim2) *
105  cc_dim1];
106 /* L107: */
107  }
108 /* L108: */
109  }
110 /* L109: */
111  }
112  i_1 = *ido;
113  for (i = 1; i <= i_1; ++i) {
114  i_2 = *l1;
115  for (k = 1; k <= i_2; ++k) {
116  ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) *
117  cc_dim1];
118 /* L110: */
119  }
120 /* L111: */
121  }
122 L112:
123  idl = 2 - *ido;
124  inc = 0;
125  i_1 = ipph;
126  for (l = 2; l <= i_1; ++l) {
127  lc = ipp2 - l;
128  idl += *ido;
129  i_2 = *idl1;
130  for (ik = 1; ik <= i_2; ++ik) {
131  c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1] * ch2[ik
132  + (ch2_dim1 << 1)];
133  c2[ik + lc * c2_dim1] = wa[idl] * ch2[ik + *ip * ch2_dim1];
134 /* L113: */
135  }
136  idlj = idl;
137  inc += *ido;
138  i_2 = ipph;
139  for (j = 3; j <= i_2; ++j) {
140  jc = ipp2 - j;
141  idlj += inc;
142  if (idlj > idp) {
143  idlj -= idp;
144  }
145  war = wa[idlj - 1];
146  wai = wa[idlj];
147  i_3 = *idl1;
148  for (ik = 1; ik <= i_3; ++ik) {
149  c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1];
150  c2[ik + lc * c2_dim1] += wai * ch2[ik + jc * ch2_dim1];
151 /* L114: */
152  }
153 /* L115: */
154  }
155 /* L116: */
156  }
157  i_1 = ipph;
158  for (j = 2; j <= i_1; ++j) {
159  i_2 = *idl1;
160  for (ik = 1; ik <= i_2; ++ik) {
161  ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
162 /* L117: */
163  }
164 /* L118: */
165  }
166  i_1 = ipph;
167  for (j = 2; j <= i_1; ++j) {
168  jc = ipp2 - j;
169  i_2 = *idl1;
170  for (ik = 2; ik <= i_2; ik += 2) {
171  ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - c2[ik +
172  jc * c2_dim1];
173  ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + c2[ik +
174  jc * c2_dim1];
175  ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + c2[ik - 1 + jc *
176  c2_dim1];
177  ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - c2[ik - 1 + jc *
178  c2_dim1];
179 /* L119: */
180  }
181 /* L120: */
182  }
183  *nac = 1;
184  if (*ido == 2) {
185  return 0;
186  }
187  *nac = 0;
188  i_1 = *idl1;
189  for (ik = 1; ik <= i_1; ++ik) {
190  c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
191 /* L121: */
192  }
193  i_1 = *ip;
194  for (j = 2; j <= i_1; ++j) {
195  i_2 = *l1;
196  for (k = 1; k <= i_2; ++k) {
197  c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) *
198  ch_dim1 + 1];
199  c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) *
200  ch_dim1 + 2];
201 /* L122: */
202  }
203 /* L123: */
204  }
205  if (idot > *l1) {
206  goto L127;
207  }
208  idij = 0;
209  i_1 = *ip;
210  for (j = 2; j <= i_1; ++j) {
211  idij += 2;
212  i_2 = *ido;
213  for (i = 4; i <= i_2; i += 2) {
214  idij += 2;
215  i_3 = *l1;
216  for (k = 1; k <= i_3; ++k) {
217  c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i
218  - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i
219  + (k + j * ch_dim2) * ch_dim1];
220  c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i + (
221  k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i - 1 + (
222  k + j * ch_dim2) * ch_dim1];
223 /* L124: */
224  }
225 /* L125: */
226  }
227 /* L126: */
228  }
229  return 0;
230 L127:
231  idj = 2 - *ido;
232  i_1 = *ip;
233  for (j = 2; j <= i_1; ++j) {
234  idj += *ido;
235  i_2 = *l1;
236  for (k = 1; k <= i_2; ++k) {
237  idij = idj;
238  i_3 = *ido;
239  for (i = 4; i <= i_3; i += 2) {
240  idij += 2;
241  c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i
242  - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i
243  + (k + j * ch_dim2) * ch_dim1];
244  c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i + (
245  k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i - 1 + (
246  k + j * ch_dim2) * ch_dim1];
247 /* L128: */
248  }
249 /* L129: */
250  }
251 /* L130: */
252  }
253  return 0;
254 } /* passb_ */
255 
256 
257 /* Subroutine */ int passb2(int *ido, int *l1, double *cc,
258  double *ch, double *wa1)
259 {
260  /* System generated locals */
261  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
262 
263  /* Local variables */
264  static int i, k;
265  static double ti2, tr2;
266 
267  /* Parameter adjustments */
268  cc_dim1 = *ido;
269  cc_offset = cc_dim1 * 3 + 1;
270  cc -= cc_offset;
271  ch_dim1 = *ido;
272  ch_dim2 = *l1;
273  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
274  ch -= ch_offset;
275  --wa1;
276 
277  /* Function Body */
278  if (*ido > 2) {
279  goto L102;
280  }
281  i_1 = *l1;
282  for (k = 1; k <= i_1; ++k) {
283  ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] +
284  cc[((k << 1) + 2) * cc_dim1 + 1];
285  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1
286  + 1] - cc[((k << 1) + 2) * cc_dim1 + 1];
287  ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2] +
288  cc[((k << 1) + 2) * cc_dim1 + 2];
289  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1
290  + 2] - cc[((k << 1) + 2) * cc_dim1 + 2];
291 /* L101: */
292  }
293  return 0;
294 L102:
295  i_1 = *l1;
296  for (k = 1; k <= i_1; ++k) {
297  i_2 = *ido;
298  for (i = 2; i <= i_2; i += 2) {
299  ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + ((k << 1) + 1) *
300  cc_dim1] + cc[i - 1 + ((k << 1) + 2) * cc_dim1];
301  tr2 = cc[i - 1 + ((k << 1) + 1) * cc_dim1] - cc[i - 1 + ((k << 1)
302  + 2) * cc_dim1];
303  ch[i + (k + ch_dim2) * ch_dim1] = cc[i + ((k << 1) + 1) * cc_dim1]
304  + cc[i + ((k << 1) + 2) * cc_dim1];
305  ti2 = cc[i + ((k << 1) + 1) * cc_dim1] - cc[i + ((k << 1) + 2) *
306  cc_dim1];
307  ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ti2 + wa1[i]
308  * tr2;
309  ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * tr2 -
310  wa1[i] * ti2;
311 /* L103: */
312  }
313 /* L104: */
314  }
315  return 0;
316 } /* passb2_ */
317 
318 
319 /* Subroutine */ int passb3(int *ido, int *l1, double *cc,
320  double *ch, double *wa1, double *wa2)
321 {
322  /* Initialized data */
323 
324  static double taur = -.5;
325  static double taui = .866025403784439;
326 
327  /* System generated locals */
328  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
329 
330  /* Local variables */
331  static int i, k;
332  static double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
333 
334  /* Parameter adjustments */
335  cc_dim1 = *ido;
336  cc_offset = (cc_dim1 << 2) + 1;
337  cc -= cc_offset;
338  ch_dim1 = *ido;
339  ch_dim2 = *l1;
340  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
341  ch -= ch_offset;
342  --wa1;
343  --wa2;
344 
345  /* Function Body */
346  if (*ido != 2) {
347  goto L102;
348  }
349  i_1 = *l1;
350  for (k = 1; k <= i_1; ++k) {
351  tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1];
352  cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
353  ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
354 
355  ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2];
356  ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2;
357  ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2;
358 
359  cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1] - cc[(k * 3 + 3) *
360  cc_dim1 + 1]);
361  ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2] - cc[(k * 3 + 3) *
362  cc_dim1 + 2]);
363  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
364  ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
365  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3;
366  ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3;
367 /* L101: */
368  }
369  return 0;
370 L102:
371  i_1 = *l1;
372  for (k = 1; k <= i_1; ++k) {
373  i_2 = *ido;
374  for (i = 2; i <= i_2; i += 2) {
375  tr2 = cc[i - 1 + (k * 3 + 2) * cc_dim1] + cc[i - 1 + (k * 3 + 3) *
376  cc_dim1];
377  cr2 = cc[i - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
378  ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 3 + 1) *
379  cc_dim1] + tr2;
380  ti2 = cc[i + (k * 3 + 2) * cc_dim1] + cc[i + (k * 3 + 3) *
381  cc_dim1];
382  ci2 = cc[i + (k * 3 + 1) * cc_dim1] + taur * ti2;
383  ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 3 + 1) * cc_dim1] +
384  ti2;
385  cr3 = taui * (cc[i - 1 + (k * 3 + 2) * cc_dim1] - cc[i - 1 + (k *
386  3 + 3) * cc_dim1]);
387  ci3 = taui * (cc[i + (k * 3 + 2) * cc_dim1] - cc[i + (k * 3 + 3) *
388  cc_dim1]);
389  dr2 = cr2 - ci3;
390  dr3 = cr2 + ci3;
391  di2 = ci2 + cr3;
392  di3 = ci2 - cr3;
393  ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 + wa1[i]
394  * dr2;
395  ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 -
396  wa1[i] * di2;
397  ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 + wa2[i] *
398  dr3;
399  ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 - wa2[
400  i] * di3;
401 /* L103: */
402  }
403 /* L104: */
404  }
405  return 0;
406 } /* passb3_ */
407 
408 /* Subroutine */ int passb4(int *ido, int *l1, double *cc,
409  double *ch, double *wa1, double *wa2, double *wa3)
410 {
411  /* System generated locals */
412  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
413 
414  /* Local variables */
415  static int i, k;
416  static double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1,
417  tr2, tr3, tr4;
418 
419  /* Parameter adjustments */
420  cc_dim1 = *ido;
421  cc_offset = cc_dim1 * 5 + 1;
422  cc -= cc_offset;
423  ch_dim1 = *ido;
424  ch_dim2 = *l1;
425  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
426  ch -= ch_offset;
427  --wa1;
428  --wa2;
429  --wa3;
430 
431  /* Function Body */
432  if (*ido != 2) {
433  goto L102;
434  }
435  i_1 = *l1;
436  for (k = 1; k <= i_1; ++k) {
437  ti1 = cc[((k << 2) + 1) * cc_dim1 + 2] - cc[((k << 2) + 3) * cc_dim1
438  + 2];
439  ti2 = cc[((k << 2) + 1) * cc_dim1 + 2] + cc[((k << 2) + 3) * cc_dim1
440  + 2];
441  tr4 = cc[((k << 2) + 4) * cc_dim1 + 2] - cc[((k << 2) + 2) * cc_dim1
442  + 2];
443  ti3 = cc[((k << 2) + 2) * cc_dim1 + 2] + cc[((k << 2) + 4) * cc_dim1
444  + 2];
445  tr1 = cc[((k << 2) + 1) * cc_dim1 + 1] - cc[((k << 2) + 3) * cc_dim1
446  + 1];
447  tr2 = cc[((k << 2) + 1) * cc_dim1 + 1] + cc[((k << 2) + 3) * cc_dim1
448  + 1];
449  ti4 = cc[((k << 2) + 2) * cc_dim1 + 1] - cc[((k << 2) + 4) * cc_dim1
450  + 1];
451  tr3 = cc[((k << 2) + 2) * cc_dim1 + 1] + cc[((k << 2) + 4) * cc_dim1
452  + 1];
453  ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
454  ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
455  ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3;
456  ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3;
457  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4;
458  ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4;
459  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4;
460  ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4;
461 /* L101: */
462  }
463  return 0;
464 L102:
465  i_1 = *l1;
466  for (k = 1; k <= i_1; ++k) {
467  i_2 = *ido;
468  for (i = 2; i <= i_2; i += 2) {
469  ti1 = cc[i + ((k << 2) + 1) * cc_dim1] - cc[i + ((k << 2) + 3) *
470  cc_dim1];
471  ti2 = cc[i + ((k << 2) + 1) * cc_dim1] + cc[i + ((k << 2) + 3) *
472  cc_dim1];
473  ti3 = cc[i + ((k << 2) + 2) * cc_dim1] + cc[i + ((k << 2) + 4) *
474  cc_dim1];
475  tr4 = cc[i + ((k << 2) + 4) * cc_dim1] - cc[i + ((k << 2) + 2) *
476  cc_dim1];
477  tr1 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] - cc[i - 1 + ((k << 2)
478  + 3) * cc_dim1];
479  tr2 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] + cc[i - 1 + ((k << 2)
480  + 3) * cc_dim1];
481  ti4 = cc[i - 1 + ((k << 2) + 2) * cc_dim1] - cc[i - 1 + ((k << 2)
482  + 4) * cc_dim1];
483  tr3 = cc[i - 1 + ((k << 2) + 2) * cc_dim1] + cc[i - 1 + ((k << 2)
484  + 4) * cc_dim1];
485  ch[i - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
486  cr3 = tr2 - tr3;
487  ch[i + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
488  ci3 = ti2 - ti3;
489  cr2 = tr1 + tr4;
490  cr4 = tr1 - tr4;
491  ci2 = ti1 + ti4;
492  ci4 = ti1 - ti4;
493  ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * cr2 -
494  wa1[i] * ci2;
495  ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ci2 + wa1[i]
496  * cr2;
497  ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * cr3 - wa2[
498  i] * ci3;
499  ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * ci3 + wa2[i] *
500  cr3;
501  ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * cr4 -
502  wa3[i] * ci4;
503  ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * ci4 + wa3[i]
504  * cr4;
505 /* L103: */
506  }
507 /* L104: */
508  }
509  return 0;
510 } /* passb4_ */
511 
512 /* Subroutine */ int passb5(int *ido, int *l1, double *cc,
513  double *ch, double *wa1, double *wa2, double *wa3,
514  double *wa4)
515 {
516  /* Initialized data */
517 
518  static double tr11 = .309016994374947;
519  static double ti11 = .951056516295154;
520  static double tr12 = -.809016994374947;
521  static double ti12 = .587785252292473;
522 
523  /* System generated locals */
524  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
525 
526  /* Local variables */
527  static int i, k;
528  static double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5,
529  cr4, ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
530 
531  /* Parameter adjustments */
532  cc_dim1 = *ido;
533  cc_offset = cc_dim1 * 6 + 1;
534  cc -= cc_offset;
535  ch_dim1 = *ido;
536  ch_dim2 = *l1;
537  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
538  ch -= ch_offset;
539  --wa1;
540  --wa2;
541  --wa3;
542  --wa4;
543 
544  /* Function Body */
545  if (*ido != 2) {
546  goto L102;
547  }
548  i_1 = *l1;
549  for (k = 1; k <= i_1; ++k) {
550  ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2];
551  ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2];
552  ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2];
553  ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2];
554  tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1];
555  tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
556  tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1];
557  tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1];
558  ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1] + tr2
559  + tr3;
560  ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2] + ti2
561  + ti3;
562  cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
563  ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3;
564  cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
565  ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3;
566  cr5 = ti11 * tr5 + ti12 * tr4;
567  ci5 = ti11 * ti5 + ti12 * ti4;
568  cr4 = ti12 * tr5 - ti11 * tr4;
569  ci4 = ti12 * ti5 - ti11 * ti4;
570  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
571  ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
572  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5;
573  ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4;
574  ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
575  ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
576  ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4;
577  ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5;
578 /* L101: */
579  }
580  return 0;
581 L102:
582  i_1 = *l1;
583  for (k = 1; k <= i_1; ++k) {
584  i_2 = *ido;
585  for (i = 2; i <= i_2; i += 2) {
586  ti5 = cc[i + (k * 5 + 2) * cc_dim1] - cc[i + (k * 5 + 5) *
587  cc_dim1];
588  ti2 = cc[i + (k * 5 + 2) * cc_dim1] + cc[i + (k * 5 + 5) *
589  cc_dim1];
590  ti4 = cc[i + (k * 5 + 3) * cc_dim1] - cc[i + (k * 5 + 4) *
591  cc_dim1];
592  ti3 = cc[i + (k * 5 + 3) * cc_dim1] + cc[i + (k * 5 + 4) *
593  cc_dim1];
594  tr5 = cc[i - 1 + (k * 5 + 2) * cc_dim1] - cc[i - 1 + (k * 5 + 5) *
595  cc_dim1];
596  tr2 = cc[i - 1 + (k * 5 + 2) * cc_dim1] + cc[i - 1 + (k * 5 + 5) *
597  cc_dim1];
598  tr4 = cc[i - 1 + (k * 5 + 3) * cc_dim1] - cc[i - 1 + (k * 5 + 4) *
599  cc_dim1];
600  tr3 = cc[i - 1 + (k * 5 + 3) * cc_dim1] + cc[i - 1 + (k * 5 + 4) *
601  cc_dim1];
602  ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 5 + 1) *
603  cc_dim1] + tr2 + tr3;
604  ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 5 + 1) * cc_dim1] +
605  ti2 + ti3;
606  cr2 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 * tr3;
607 
608  ci2 = cc[i + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;
609  cr3 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 * tr3;
610 
611  ci3 = cc[i + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;
612  cr5 = ti11 * tr5 + ti12 * tr4;
613  ci5 = ti11 * ti5 + ti12 * ti4;
614  cr4 = ti12 * tr5 - ti11 * tr4;
615  ci4 = ti12 * ti5 - ti11 * ti4;
616  dr3 = cr3 - ci4;
617  dr4 = cr3 + ci4;
618  di3 = ci3 + cr4;
619  di4 = ci3 - cr4;
620  dr5 = cr2 + ci5;
621  dr2 = cr2 - ci5;
622  di5 = ci2 - cr5;
623  di2 = ci2 + cr5;
624  ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 -
625  wa1[i] * di2;
626  ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 + wa1[i]
627  * dr2;
628  ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 - wa2[
629  i] * di3;
630  ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 + wa2[i] *
631  dr3;
632  ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * dr4 -
633  wa3[i] * di4;
634  ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * di4 + wa3[i]
635  * dr4;
636  ch[i - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * dr5 - wa4[
637  i] * di5;
638  ch[i + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * di5 + wa4[i] *
639  dr5;
640 /* L103: */
641  }
642 /* L104: */
643  }
644  return 0;
645 } /* passb5_ */
646 
647 /* Subroutine */ int passf(int *nac, int *ido, int *ip, int *
648  l1, int *idl1, double *cc, double *c1, double *c2,
649  double *ch, double *ch2, double *wa)
650 {
651  /* System generated locals */
652  int ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1,
653  c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset,
654  i_1, i_2, i_3;
655 
656  /* Local variables */
657  static int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj,
658  idl, inc, idp;
659  static double wai, war;
660  static int ipp2;
661 
662  /* Parameter adjustments */
663  cc_dim1 = *ido;
664  cc_dim2 = *ip;
665  cc_offset = cc_dim1 * (cc_dim2 + 1) + 1;
666  cc -= cc_offset;
667  c1_dim1 = *ido;
668  c1_dim2 = *l1;
669  c1_offset = c1_dim1 * (c1_dim2 + 1) + 1;
670  c1 -= c1_offset;
671  c2_dim1 = *idl1;
672  c2_offset = c2_dim1 + 1;
673  c2 -= c2_offset;
674  ch_dim1 = *ido;
675  ch_dim2 = *l1;
676  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
677  ch -= ch_offset;
678  ch2_dim1 = *idl1;
679  ch2_offset = ch2_dim1 + 1;
680  ch2 -= ch2_offset;
681  --wa;
682 
683  /* Function Body */
684  idot = *ido / 2;
685  nt = *ip * *idl1;
686  ipp2 = *ip + 2;
687  ipph = (*ip + 1) / 2;
688  idp = *ip * *ido;
689 
690  if (*ido < *l1) {
691  goto L106;
692  }
693  i_1 = ipph;
694  for (j = 2; j <= i_1; ++j) {
695  jc = ipp2 - j;
696  i_2 = *l1;
697  for (k = 1; k <= i_2; ++k) {
698  i_3 = *ido;
699  for (i = 1; i <= i_3; ++i) {
700  ch[i + (k + j * ch_dim2) * ch_dim1] = cc[i + (j + k * cc_dim2)
701  * cc_dim1] + cc[i + (jc + k * cc_dim2) * cc_dim1];
702  ch[i + (k + jc * ch_dim2) * ch_dim1] = cc[i + (j + k *
703  cc_dim2) * cc_dim1] - cc[i + (jc + k * cc_dim2) *
704  cc_dim1];
705 /* L101: */
706  }
707 /* L102: */
708  }
709 /* L103: */
710  }
711  i_1 = *l1;
712  for (k = 1; k <= i_1; ++k) {
713  i_2 = *ido;
714  for (i = 1; i <= i_2; ++i) {
715  ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) *
716  cc_dim1];
717 /* L104: */
718  }
719 /* L105: */
720  }
721  goto L112;
722 L106:
723  i_1 = ipph;
724  for (j = 2; j <= i_1; ++j) {
725  jc = ipp2 - j;
726  i_2 = *ido;
727  for (i = 1; i <= i_2; ++i) {
728  i_3 = *l1;
729  for (k = 1; k <= i_3; ++k) {
730  ch[i + (k + j * ch_dim2) * ch_dim1] = cc[i + (j + k * cc_dim2)
731  * cc_dim1] + cc[i + (jc + k * cc_dim2) * cc_dim1];
732  ch[i + (k + jc * ch_dim2) * ch_dim1] = cc[i + (j + k *
733  cc_dim2) * cc_dim1] - cc[i + (jc + k * cc_dim2) *
734  cc_dim1];
735 /* L107: */
736  }
737 /* L108: */
738  }
739 /* L109: */
740  }
741  i_1 = *ido;
742  for (i = 1; i <= i_1; ++i) {
743  i_2 = *l1;
744  for (k = 1; k <= i_2; ++k) {
745  ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) *
746  cc_dim1];
747 /* L110: */
748  }
749 /* L111: */
750  }
751 L112:
752  idl = 2 - *ido;
753  inc = 0;
754  i_1 = ipph;
755  for (l = 2; l <= i_1; ++l) {
756  lc = ipp2 - l;
757  idl += *ido;
758  i_2 = *idl1;
759  for (ik = 1; ik <= i_2; ++ik) {
760  c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1] * ch2[ik
761  + (ch2_dim1 << 1)];
762  c2[ik + lc * c2_dim1] = -wa[idl] * ch2[ik + *ip * ch2_dim1];
763 /* L113: */
764  }
765  idlj = idl;
766  inc += *ido;
767  i_2 = ipph;
768  for (j = 3; j <= i_2; ++j) {
769  jc = ipp2 - j;
770  idlj += inc;
771  if (idlj > idp) {
772  idlj -= idp;
773  }
774  war = wa[idlj - 1];
775  wai = wa[idlj];
776  i_3 = *idl1;
777  for (ik = 1; ik <= i_3; ++ik) {
778  c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1];
779  c2[ik + lc * c2_dim1] -= wai * ch2[ik + jc * ch2_dim1];
780 /* L114: */
781  }
782 /* L115: */
783  }
784 /* L116: */
785  }
786  i_1 = ipph;
787  for (j = 2; j <= i_1; ++j) {
788  i_2 = *idl1;
789  for (ik = 1; ik <= i_2; ++ik) {
790  ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
791 /* L117: */
792  }
793 /* L118: */
794  }
795  i_1 = ipph;
796  for (j = 2; j <= i_1; ++j) {
797  jc = ipp2 - j;
798  i_2 = *idl1;
799  for (ik = 2; ik <= i_2; ik += 2) {
800  ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - c2[ik +
801  jc * c2_dim1];
802  ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + c2[ik +
803  jc * c2_dim1];
804  ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + c2[ik - 1 + jc *
805  c2_dim1];
806  ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - c2[ik - 1 + jc *
807  c2_dim1];
808 /* L119: */
809  }
810 /* L120: */
811  }
812  *nac = 1;
813  if (*ido == 2) {
814  return 0;
815  }
816  *nac = 0;
817  i_1 = *idl1;
818  for (ik = 1; ik <= i_1; ++ik) {
819  c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
820 /* L121: */
821  }
822  i_1 = *ip;
823  for (j = 2; j <= i_1; ++j) {
824  i_2 = *l1;
825  for (k = 1; k <= i_2; ++k) {
826  c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) *
827  ch_dim1 + 1];
828  c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) *
829  ch_dim1 + 2];
830 /* L122: */
831  }
832 /* L123: */
833  }
834  if (idot > *l1) {
835  goto L127;
836  }
837  idij = 0;
838  i_1 = *ip;
839  for (j = 2; j <= i_1; ++j) {
840  idij += 2;
841  i_2 = *ido;
842  for (i = 4; i <= i_2; i += 2) {
843  idij += 2;
844  i_3 = *l1;
845  for (k = 1; k <= i_3; ++k) {
846  c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i
847  - 1 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i
848  + (k + j * ch_dim2) * ch_dim1];
849  c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i + (
850  k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i - 1 + (
851  k + j * ch_dim2) * ch_dim1];
852 /* L124: */
853  }
854 /* L125: */
855  }
856 /* L126: */
857  }
858  return 0;
859 L127:
860  idj = 2 - *ido;
861  i_1 = *ip;
862  for (j = 2; j <= i_1; ++j) {
863  idj += *ido;
864  i_2 = *l1;
865  for (k = 1; k <= i_2; ++k) {
866  idij = idj;
867  i_3 = *ido;
868  for (i = 4; i <= i_3; i += 2) {
869  idij += 2;
870  c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i
871  - 1 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i
872  + (k + j * ch_dim2) * ch_dim1];
873  c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i + (
874  k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i - 1 + (
875  k + j * ch_dim2) * ch_dim1];
876 /* L128: */
877  }
878 /* L129: */
879  }
880 /* L130: */
881  }
882  return 0;
883 } /* passf_ */
884 
885 /* Subroutine */ int passf2(int *ido, int *l1, double *cc,
886  double *ch, double *wa1)
887 {
888  /* System generated locals */
889  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
890 
891  /* Local variables */
892  static int i, k;
893  static double ti2, tr2;
894 
895  /* Parameter adjustments */
896  cc_dim1 = *ido;
897  cc_offset = cc_dim1 * 3 + 1;
898  cc -= cc_offset;
899  ch_dim1 = *ido;
900  ch_dim2 = *l1;
901  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
902  ch -= ch_offset;
903  --wa1;
904 
905  /* Function Body */
906  if (*ido > 2) {
907  goto L102;
908  }
909  i_1 = *l1;
910  for (k = 1; k <= i_1; ++k) {
911  ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] +
912  cc[((k << 1) + 2) * cc_dim1 + 1];
913  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1
914  + 1] - cc[((k << 1) + 2) * cc_dim1 + 1];
915  ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2] +
916  cc[((k << 1) + 2) * cc_dim1 + 2];
917  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1
918  + 2] - cc[((k << 1) + 2) * cc_dim1 + 2];
919 /* L101: */
920  }
921  return 0;
922 L102:
923  i_1 = *l1;
924  for (k = 1; k <= i_1; ++k) {
925  i_2 = *ido;
926  for (i = 2; i <= i_2; i += 2) {
927  ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + ((k << 1) + 1) *
928  cc_dim1] + cc[i - 1 + ((k << 1) + 2) * cc_dim1];
929  tr2 = cc[i - 1 + ((k << 1) + 1) * cc_dim1] - cc[i - 1 + ((k << 1)
930  + 2) * cc_dim1];
931  ch[i + (k + ch_dim2) * ch_dim1] = cc[i + ((k << 1) + 1) * cc_dim1]
932  + cc[i + ((k << 1) + 2) * cc_dim1];
933  ti2 = cc[i + ((k << 1) + 1) * cc_dim1] - cc[i + ((k << 1) + 2) *
934  cc_dim1];
935  ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ti2 - wa1[i]
936  * tr2;
937  ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * tr2 +
938  wa1[i] * ti2;
939 /* L103: */
940  }
941 /* L104: */
942  }
943  return 0;
944 } /* passf2_ */
945 
946 /* Subroutine */ int passf3(int *ido, int *l1, double *cc,
947  double *ch, double *wa1, double *wa2)
948 {
949  /* Initialized data */
950 
951  static double taur = -.5;
952  static double taui = -.866025403784439;
953 
954  /* System generated locals */
955  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
956 
957  /* Local variables */
958  static int i, k;
959  static double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
960 
961  /* Parameter adjustments */
962  cc_dim1 = *ido;
963  cc_offset = (cc_dim1 << 2) + 1;
964  cc -= cc_offset;
965  ch_dim1 = *ido;
966  ch_dim2 = *l1;
967  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
968  ch -= ch_offset;
969  --wa1;
970  --wa2;
971 
972  /* Function Body */
973  if (*ido != 2) {
974  goto L102;
975  }
976  i_1 = *l1;
977  for (k = 1; k <= i_1; ++k) {
978  tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1];
979  cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
980  ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
981 
982  ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2];
983  ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2;
984  ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2;
985 
986  cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1] - cc[(k * 3 + 3) *
987  cc_dim1 + 1]);
988  ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2] - cc[(k * 3 + 3) *
989  cc_dim1 + 2]);
990  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
991  ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
992  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3;
993  ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3;
994 /* L101: */
995  }
996  return 0;
997 L102:
998  i_1 = *l1;
999  for (k = 1; k <= i_1; ++k) {
1000  i_2 = *ido;
1001  for (i = 2; i <= i_2; i += 2) {
1002  tr2 = cc[i - 1 + (k * 3 + 2) * cc_dim1] + cc[i - 1 + (k * 3 + 3) *
1003  cc_dim1];
1004  cr2 = cc[i - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
1005  ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 3 + 1) *
1006  cc_dim1] + tr2;
1007  ti2 = cc[i + (k * 3 + 2) * cc_dim1] + cc[i + (k * 3 + 3) *
1008  cc_dim1];
1009  ci2 = cc[i + (k * 3 + 1) * cc_dim1] + taur * ti2;
1010  ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 3 + 1) * cc_dim1] +
1011  ti2;
1012  cr3 = taui * (cc[i - 1 + (k * 3 + 2) * cc_dim1] - cc[i - 1 + (k *
1013  3 + 3) * cc_dim1]);
1014  ci3 = taui * (cc[i + (k * 3 + 2) * cc_dim1] - cc[i + (k * 3 + 3) *
1015  cc_dim1]);
1016  dr2 = cr2 - ci3;
1017  dr3 = cr2 + ci3;
1018  di2 = ci2 + cr3;
1019  di3 = ci2 - cr3;
1020  ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 - wa1[i]
1021  * dr2;
1022  ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 +
1023  wa1[i] * di2;
1024  ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 - wa2[i] *
1025  dr3;
1026  ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 + wa2[
1027  i] * di3;
1028 /* L103: */
1029  }
1030 /* L104: */
1031  }
1032  return 0;
1033 } /* passf3_ */
1034 
1035 /* Subroutine */ int passf4(int *ido, int *l1, double *cc,
1036  double *ch, double *wa1, double *wa2, double *wa3)
1037 {
1038  /* System generated locals */
1039  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
1040 
1041  /* Local variables */
1042  static int i, k;
1043  static double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1,
1044  tr2, tr3, tr4;
1045 
1046  /* Parameter adjustments */
1047  cc_dim1 = *ido;
1048  cc_offset = cc_dim1 * 5 + 1;
1049  cc -= cc_offset;
1050  ch_dim1 = *ido;
1051  ch_dim2 = *l1;
1052  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
1053  ch -= ch_offset;
1054  --wa1;
1055  --wa2;
1056  --wa3;
1057 
1058  /* Function Body */
1059  if (*ido != 2) {
1060  goto L102;
1061  }
1062  i_1 = *l1;
1063  for (k = 1; k <= i_1; ++k) {
1064  ti1 = cc[((k << 2) + 1) * cc_dim1 + 2] - cc[((k << 2) + 3) * cc_dim1
1065  + 2];
1066  ti2 = cc[((k << 2) + 1) * cc_dim1 + 2] + cc[((k << 2) + 3) * cc_dim1
1067  + 2];
1068  tr4 = cc[((k << 2) + 2) * cc_dim1 + 2] - cc[((k << 2) + 4) * cc_dim1
1069  + 2];
1070  ti3 = cc[((k << 2) + 2) * cc_dim1 + 2] + cc[((k << 2) + 4) * cc_dim1
1071  + 2];
1072  tr1 = cc[((k << 2) + 1) * cc_dim1 + 1] - cc[((k << 2) + 3) * cc_dim1
1073  + 1];
1074  tr2 = cc[((k << 2) + 1) * cc_dim1 + 1] + cc[((k << 2) + 3) * cc_dim1
1075  + 1];
1076  ti4 = cc[((k << 2) + 4) * cc_dim1 + 1] - cc[((k << 2) + 2) * cc_dim1
1077  + 1];
1078  tr3 = cc[((k << 2) + 2) * cc_dim1 + 1] + cc[((k << 2) + 4) * cc_dim1
1079  + 1];
1080  ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
1081  ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
1082  ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3;
1083  ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3;
1084  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4;
1085  ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4;
1086  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4;
1087  ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4;
1088 /* L101: */
1089  }
1090  return 0;
1091 L102:
1092  i_1 = *l1;
1093  for (k = 1; k <= i_1; ++k) {
1094  i_2 = *ido;
1095  for (i = 2; i <= i_2; i += 2) {
1096  ti1 = cc[i + ((k << 2) + 1) * cc_dim1] - cc[i + ((k << 2) + 3) *
1097  cc_dim1];
1098  ti2 = cc[i + ((k << 2) + 1) * cc_dim1] + cc[i + ((k << 2) + 3) *
1099  cc_dim1];
1100  ti3 = cc[i + ((k << 2) + 2) * cc_dim1] + cc[i + ((k << 2) + 4) *
1101  cc_dim1];
1102  tr4 = cc[i + ((k << 2) + 2) * cc_dim1] - cc[i + ((k << 2) + 4) *
1103  cc_dim1];
1104  tr1 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] - cc[i - 1 + ((k << 2)
1105  + 3) * cc_dim1];
1106  tr2 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] + cc[i - 1 + ((k << 2)
1107  + 3) * cc_dim1];
1108  ti4 = cc[i - 1 + ((k << 2) + 4) * cc_dim1] - cc[i - 1 + ((k << 2)
1109  + 2) * cc_dim1];
1110  tr3 = cc[i - 1 + ((k << 2) + 2) * cc_dim1] + cc[i - 1 + ((k << 2)
1111  + 4) * cc_dim1];
1112  ch[i - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
1113  cr3 = tr2 - tr3;
1114  ch[i + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
1115  ci3 = ti2 - ti3;
1116  cr2 = tr1 + tr4;
1117  cr4 = tr1 - tr4;
1118  ci2 = ti1 + ti4;
1119  ci4 = ti1 - ti4;
1120  ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * cr2 +
1121  wa1[i] * ci2;
1122  ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ci2 - wa1[i]
1123  * cr2;
1124  ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * cr3 + wa2[
1125  i] * ci3;
1126  ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * ci3 - wa2[i] *
1127  cr3;
1128  ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * cr4 +
1129  wa3[i] * ci4;
1130  ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * ci4 - wa3[i]
1131  * cr4;
1132 /* L103: */
1133  }
1134 /* L104: */
1135  }
1136  return 0;
1137 } /* passf4_ */
1138 
1139 /* Subroutine */ int passf5(int *ido, int *l1, double *cc,
1140  double *ch, double *wa1, double *wa2, double *wa3,
1141  double *wa4)
1142 {
1143  /* Initialized data */
1144 
1145  static double tr11 = .309016994374947;
1146  static double ti11 = -.951056516295154;
1147  static double tr12 = -.809016994374947;
1148  static double ti12 = -.587785252292473;
1149 
1150  /* System generated locals */
1151  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
1152 
1153  /* Local variables */
1154  static int i, k;
1155  static double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5,
1156  cr4, ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
1157 
1158  /* Parameter adjustments */
1159  cc_dim1 = *ido;
1160  cc_offset = cc_dim1 * 6 + 1;
1161  cc -= cc_offset;
1162  ch_dim1 = *ido;
1163  ch_dim2 = *l1;
1164  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
1165  ch -= ch_offset;
1166  --wa1;
1167  --wa2;
1168  --wa3;
1169  --wa4;
1170 
1171  /* Function Body */
1172  if (*ido != 2) {
1173  goto L102;
1174  }
1175  i_1 = *l1;
1176  for (k = 1; k <= i_1; ++k) {
1177  ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2];
1178  ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2];
1179  ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2];
1180  ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2];
1181  tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1];
1182  tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
1183  tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1];
1184  tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1];
1185  ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1] + tr2
1186  + tr3;
1187  ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2] + ti2
1188  + ti3;
1189  cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
1190  ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3;
1191  cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
1192  ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3;
1193  cr5 = ti11 * tr5 + ti12 * tr4;
1194  ci5 = ti11 * ti5 + ti12 * ti4;
1195  cr4 = ti12 * tr5 - ti11 * tr4;
1196  ci4 = ti12 * ti5 - ti11 * ti4;
1197  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
1198  ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
1199  ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5;
1200  ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4;
1201  ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
1202  ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
1203  ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4;
1204  ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5;
1205 /* L101: */
1206  }
1207  return 0;
1208 L102:
1209  i_1 = *l1;
1210  for (k = 1; k <= i_1; ++k) {
1211  i_2 = *ido;
1212  for (i = 2; i <= i_2; i += 2) {
1213  ti5 = cc[i + (k * 5 + 2) * cc_dim1] - cc[i + (k * 5 + 5) *
1214  cc_dim1];
1215  ti2 = cc[i + (k * 5 + 2) * cc_dim1] + cc[i + (k * 5 + 5) *
1216  cc_dim1];
1217  ti4 = cc[i + (k * 5 + 3) * cc_dim1] - cc[i + (k * 5 + 4) *
1218  cc_dim1];
1219  ti3 = cc[i + (k * 5 + 3) * cc_dim1] + cc[i + (k * 5 + 4) *
1220  cc_dim1];
1221  tr5 = cc[i - 1 + (k * 5 + 2) * cc_dim1] - cc[i - 1 + (k * 5 + 5) *
1222  cc_dim1];
1223  tr2 = cc[i - 1 + (k * 5 + 2) * cc_dim1] + cc[i - 1 + (k * 5 + 5) *
1224  cc_dim1];
1225  tr4 = cc[i - 1 + (k * 5 + 3) * cc_dim1] - cc[i - 1 + (k * 5 + 4) *
1226  cc_dim1];
1227  tr3 = cc[i - 1 + (k * 5 + 3) * cc_dim1] + cc[i - 1 + (k * 5 + 4) *
1228  cc_dim1];
1229  ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 5 + 1) *
1230  cc_dim1] + tr2 + tr3;
1231  ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 5 + 1) * cc_dim1] +
1232  ti2 + ti3;
1233  cr2 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 * tr3;
1234 
1235  ci2 = cc[i + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;
1236  cr3 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 * tr3;
1237 
1238  ci3 = cc[i + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;
1239  cr5 = ti11 * tr5 + ti12 * tr4;
1240  ci5 = ti11 * ti5 + ti12 * ti4;
1241  cr4 = ti12 * tr5 - ti11 * tr4;
1242  ci4 = ti12 * ti5 - ti11 * ti4;
1243  dr3 = cr3 - ci4;
1244  dr4 = cr3 + ci4;
1245  di3 = ci3 + cr4;
1246  di4 = ci3 - cr4;
1247  dr5 = cr2 + ci5;
1248  dr2 = cr2 - ci5;
1249  di5 = ci2 - cr5;
1250  di2 = ci2 + cr5;
1251  ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 +
1252  wa1[i] * di2;
1253  ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 - wa1[i]
1254  * dr2;
1255  ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 + wa2[
1256  i] * di3;
1257  ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 - wa2[i] *
1258  dr3;
1259  ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * dr4 +
1260  wa3[i] * di4;
1261  ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * di4 - wa3[i]
1262  * dr4;
1263  ch[i - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * dr5 + wa4[
1264  i] * di5;
1265  ch[i + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * di5 - wa4[i] *
1266  dr5;
1267 /* L103: */
1268  }
1269 /* L104: */
1270  }
1271  return 0;
1272 } /* passf5_ */
1273 
1274 /* Subroutine */ int passb4(int *ido, int *l1, double *cc,
1275  double *ch, double *wa1, double *wa2, double *wa3);
1276 
1277 /* Subroutine */ int cfftb1(int *n, double *c, double *ch,
1278  double *wa, int *ifac)
1279 {
1280  /* System generated locals */
1281  int i_1;
1282 
1283  /* Local variables */
1284  static int idot, i;
1285  static int k1, l1, l2, n2;
1286  static int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
1287 
1288  /* Parameter adjustments */
1289  --c;
1290  --ch;
1291  --wa;
1292  --ifac;
1293 
1294  /* Function Body */
1295  nf = ifac[2];
1296  na = 0;
1297  l1 = 1;
1298  iw = 1;
1299  i_1 = nf;
1300  for (k1 = 1; k1 <= i_1; ++k1) {
1301  ip = ifac[k1 + 2];
1302  l2 = ip * l1;
1303  ido = *n / l2;
1304  idot = ido + ido;
1305  idl1 = idot * l1;
1306  if (ip != 4) {
1307  goto L103;
1308  }
1309  ix2 = iw + idot;
1310  ix3 = ix2 + idot;
1311  if (na != 0) {
1312  goto L101;
1313  }
1314  passb4(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);
1315  goto L102;
1316 L101:
1317  passb4(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3]);
1318 L102:
1319  na = 1 - na;
1320  goto L115;
1321 L103:
1322  if (ip != 2) {
1323  goto L106;
1324  }
1325  if (na != 0) {
1326  goto L104;
1327  }
1328  passb2(&idot, &l1, &c[1], &ch[1], &wa[iw]);
1329  goto L105;
1330 L104:
1331  passb2(&idot, &l1, &ch[1], &c[1], &wa[iw]);
1332 L105:
1333  na = 1 - na;
1334  goto L115;
1335 L106:
1336  if (ip != 3) {
1337  goto L109;
1338  }
1339  ix2 = iw + idot;
1340  if (na != 0) {
1341  goto L107;
1342  }
1343  passb3(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2]);
1344  goto L108;
1345 L107:
1346  passb3(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2]);
1347 L108:
1348  na = 1 - na;
1349  goto L115;
1350 L109:
1351  if (ip != 5) {
1352  goto L112;
1353  }
1354  ix2 = iw + idot;
1355  ix3 = ix2 + idot;
1356  ix4 = ix3 + idot;
1357  if (na != 0) {
1358  goto L110;
1359  }
1360  passb5(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
1361  ix4]);
1362  goto L111;
1363 L110:
1364  passb5(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
1365  ix4]);
1366 L111:
1367  na = 1 - na;
1368  goto L115;
1369 L112:
1370  if (na != 0) {
1371  goto L113;
1372  }
1373  passb(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1], &ch[
1374  1], &wa[iw]);
1375  goto L114;
1376 L113:
1377  passb(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1], &
1378  c[1], &wa[iw]);
1379 L114:
1380  if (nac != 0) {
1381  na = 1 - na;
1382  }
1383 L115:
1384  l1 = l2;
1385  iw += (ip - 1) * idot;
1386 /* L116: */
1387  }
1388  if (na == 0) {
1389  return 0;
1390  }
1391  n2 = *n + *n;
1392  i_1 = n2;
1393  for (i = 1; i <= i_1; ++i) {
1394  c[i] = ch[i];
1395 /* L117: */
1396  }
1397  return 0;
1398 } /* cfftb1_ */
1399 
1400 
1401 /* Subroutine */ int cfftb(int *n, double *c, double *wsave)
1402 {
1403  static int iw1, iw2;
1404 
1405  /* Parameter adjustments */
1406  --c;
1407  --wsave;
1408 
1409  /* Function Body */
1410  if (*n == 1) {
1411  return 0;
1412  }
1413  iw1 = *n + *n + 1;
1414  iw2 = iw1 + *n + *n;
1415  cfftb1(n, &c[1], &wsave[1], &wsave[iw1], (int*) &wsave[iw2]);
1416  return 0;
1417 } /* cfftb_ */
1418 
1419 
1420 /* Subroutine */ int cfftf1(int *n, double *c, double *ch,
1421  double *wa, int *ifac)
1422 {
1423  /* System generated locals */
1424  int i_1;
1425 
1426  /* Local variables */
1427  static int idot, i;
1428  static int k1, l1, l2, n2;
1429  static int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
1430 
1431  /* Parameter adjustments */
1432  --c;
1433  --ch;
1434  --wa;
1435  --ifac;
1436 
1437  /* Function Body */
1438  nf = ifac[2];
1439  na = 0;
1440  l1 = 1;
1441  iw = 1;
1442  i_1 = nf;
1443  for (k1 = 1; k1 <= i_1; ++k1) {
1444  ip = ifac[k1 + 2];
1445  l2 = ip * l1;
1446  ido = *n / l2;
1447  idot = ido + ido;
1448  idl1 = idot * l1;
1449  if (ip != 4) {
1450  goto L103;
1451  }
1452  ix2 = iw + idot;
1453  ix3 = ix2 + idot;
1454  if (na != 0) {
1455  goto L101;
1456  }
1457  passf4(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);
1458  goto L102;
1459 L101:
1460  passf4(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3]);
1461 L102:
1462  na = 1 - na;
1463  goto L115;
1464 L103:
1465  if (ip != 2) {
1466  goto L106;
1467  }
1468  if (na != 0) {
1469  goto L104;
1470  }
1471  passf2(&idot, &l1, &c[1], &ch[1], &wa[iw]);
1472  goto L105;
1473 L104:
1474  passf2(&idot, &l1, &ch[1], &c[1], &wa[iw]);
1475 L105:
1476  na = 1 - na;
1477  goto L115;
1478 L106:
1479  if (ip != 3) {
1480  goto L109;
1481  }
1482  ix2 = iw + idot;
1483  if (na != 0) {
1484  goto L107;
1485  }
1486  passf3(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2]);
1487  goto L108;
1488 L107:
1489  passf3(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2]);
1490 L108:
1491  na = 1 - na;
1492  goto L115;
1493 L109:
1494  if (ip != 5) {
1495  goto L112;
1496  }
1497  ix2 = iw + idot;
1498  ix3 = ix2 + idot;
1499  ix4 = ix3 + idot;
1500  if (na != 0) {
1501  goto L110;
1502  }
1503  passf5(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
1504  ix4]);
1505  goto L111;
1506 L110:
1507  passf5(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
1508  ix4]);
1509 L111:
1510  na = 1 - na;
1511  goto L115;
1512 L112:
1513  if (na != 0) {
1514  goto L113;
1515  }
1516  passf(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1], &ch[
1517  1], &wa[iw]);
1518  goto L114;
1519 L113:
1520  passf(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1], &
1521  c[1], &wa[iw]);
1522 L114:
1523  if (nac != 0) {
1524  na = 1 - na;
1525  }
1526 L115:
1527  l1 = l2;
1528  iw += (ip - 1) * idot;
1529 /* L116: */
1530  }
1531  if (na == 0) {
1532  return 0;
1533  }
1534  n2 = *n + *n;
1535  i_1 = n2;
1536  for (i = 1; i <= i_1; ++i) {
1537  c[i] = ch[i];
1538 /* L117: */
1539  }
1540  return 0;
1541 } /* cfftf1_ */
1542 
1543 
1544 /* Subroutine */ int cfftf(int *n, double *c, double *wsave)
1545 {
1546  static int iw1, iw2;
1547 
1548  /* Parameter adjustments */
1549  --c;
1550  --wsave;
1551 
1552  /* Function Body */
1553  if (*n == 1) {
1554  return 0;
1555  }
1556  iw1 = *n + *n + 1;
1557  iw2 = iw1 + *n + *n;
1558  cfftf1(n, &c[1], &wsave[1], &wsave[iw1], (int*) &wsave[iw2]);
1559  return 0;
1560 } /* cfftf_ */
1561 
1562 
1563 /* Subroutine */ int cffti1(int *n, double *wa, int *ifac)
1564 {
1565  /* Initialized data */
1566 
1567  static int ntryh[4] = { 3,4,2,5 };
1568 
1569  /* System generated locals */
1570  int i_1, i_2, i_3;
1571 
1572  /* Local variables */
1573  static double argh;
1574  static int idot, ntry, i, j;
1575  static double argld;
1576  static int i1, k1, l1, l2, ib;
1577  static double fi;
1578  static int ld, ii, nf, ip, nl, nq, nr;
1579  static double arg;
1580  static int ido, ipm;
1581  static double tpi;
1582 
1583  /* Parameter adjustments */
1584  --wa;
1585  --ifac;
1586 
1587  /* Function Body */
1588  nl = *n;
1589  nf = 0;
1590  j = 0;
1591 L101:
1592  ++j;
1593  if (j - 4 <= 0) {
1594  goto L102;
1595  } else {
1596  goto L103;
1597  }
1598 L102:
1599  ntry = ntryh[j - 1];
1600  goto L104;
1601 L103:
1602  ntry += 2;
1603 L104:
1604  nq = nl / ntry;
1605  nr = nl - ntry * nq;
1606  if (nr != 0) {
1607  goto L101;
1608  } else {
1609  goto L105;
1610  }
1611 L105:
1612  ++nf;
1613  ifac[nf + 2] = ntry;
1614  nl = nq;
1615  if (ntry != 2) {
1616  goto L107;
1617  }
1618  if (nf == 1) {
1619  goto L107;
1620  }
1621  i_1 = nf;
1622  for (i = 2; i <= i_1; ++i) {
1623  ib = nf - i + 2;
1624  ifac[ib + 2] = ifac[ib + 1];
1625 /* L106: */
1626  }
1627  ifac[3] = 2;
1628 L107:
1629  if (nl != 1) {
1630  goto L104;
1631  }
1632  ifac[1] = *n;
1633  ifac[2] = nf;
1634  tpi = 6.28318530717959;
1635  argh = tpi / (double) (*n);
1636  i = 2;
1637  l1 = 1;
1638  i_1 = nf;
1639  for (k1 = 1; k1 <= i_1; ++k1) {
1640  ip = ifac[k1 + 2];
1641  ld = 0;
1642  l2 = l1 * ip;
1643  ido = *n / l2;
1644  idot = ido + ido + 2;
1645  ipm = ip - 1;
1646  i_2 = ipm;
1647  for (j = 1; j <= i_2; ++j) {
1648  i1 = i;
1649  wa[i - 1] = 1.;
1650  wa[i] = 0.;
1651  ld += l1;
1652  fi = 0.;
1653  argld = (double) ld * argh;
1654  i_3 = idot;
1655  for (ii = 4; ii <= i_3; ii += 2) {
1656  i += 2;
1657  fi += 1.;
1658  arg = fi * argld;
1659  wa[i - 1] = cos(arg);
1660  wa[i] = sin(arg);
1661 /* L108: */
1662  }
1663  if (ip <= 5) {
1664  goto L109;
1665  }
1666  wa[i1 - 1] = wa[i - 1];
1667  wa[i1] = wa[i];
1668 L109:
1669  ;}
1670  l1 = l2;
1671 /* L110: */
1672  }
1673  return 0;
1674 } /* cffti1_ */
1675 
1676 /* Subroutine */ int cffti(int *n, double *wsave)
1677 {
1678  static int iw1, iw2;
1679 
1680  /* Parameter adjustments */
1681  --wsave;
1682 
1683  /* Function Body */
1684  if (*n == 1) {
1685  return 0;
1686  }
1687  iw1 = *n + *n + 1;
1688  iw2 = iw1 + *n + *n;
1689  cffti1(n, &wsave[iw1], (int*) &wsave[iw2]);
1690  return 0;
1691 } /* cffti_ */
1692 
1693 
1694 /* typedef struct { double r, i; } doublecomplex; */
1695 
1696 /****************************************************************************
1697 **/
1698 
1699 /* 3D (slow) Fourier Transform */
1700 /* this 1d->3d code is brute force approach */
1701 /* the 1d code is a double precision version of fftpack from netlib */
1702 /* due to Paul N Swartztrauber at NCAR Boulder Coloraso */
1703 
1704 /****************************************************************************
1705 **/
1706 /* Subroutine */ int pubz3di(int *n1, int *n2, int *n3,
1707  double *table, int *ntable)
1708 {
1709  /* System generated locals */
1710  int table_dim1, table_offset;
1711 
1712  /* Local variables */
1713 
1714  /* Parameter adjustments */
1715  table_dim1 = *ntable;
1716  table_offset = table_dim1 + 1;
1717  table -= table_offset;
1718 
1719  /* Function Body */
1720 /* ntable should be 4*max(n1,n2,n3) +15 */
1721  cffti(n1, &table[table_dim1 + 1]);
1722  cffti(n2, &table[(table_dim1 << 1) + 1]);
1723  cffti(n3, &table[table_dim1 * 3 + 1]);
1724  return 0;
1725 } /* pubz3di_ */
1726 
1727 /****************************************************************************
1728 **/
1729 /* Subroutine */ int pubz3d(int *isign, int *n1, int *n2,
1730  int *n3, doublecomplex *w, int *ld1, int *ld2, double
1731  *table, int *ntable, doublecomplex *work /*, int * DUMMY1 nwork */)
1732 {
1733  /* System generated locals */
1734  int w_dim1, w_dim2, w_offset, table_dim1, table_offset, i_1, i_2, i_3,
1735  i_4, i_5;
1736 
1737  /* Local variables */
1738  static int i, j, k;
1739 
1740  /* Parameter adjustments */
1741  w_dim1 = *ld1;
1742  w_dim2 = *ld2;
1743  w_offset = w_dim1 * (w_dim2 + 1) + 1;
1744  w -= w_offset;
1745  table_dim1 = *ntable;
1746  table_offset = table_dim1 + 1;
1747  table -= table_offset;
1748  --work;
1749 
1750  /* Function Body */
1751 /* ntable should be 4*max(n1,n2,n3) +15 */
1752 /* nwork should be max(n1,n2,n3) */
1753 
1754 /* transform along X first ... */
1755 
1756  i_1 = *n3;
1757  for (k = 1; k <= i_1; ++k) {
1758  i_2 = *n2;
1759  for (j = 1; j <= i_2; ++j) {
1760  i_3 = *n1;
1761  for (i = 1; i <= i_3; ++i) {
1762  i_4 = i;
1763  i_5 = i + (j + k * w_dim2) * w_dim1;
1764  work[i_4].r = w[i_5].r, work[i_4].i = w[i_5].i;
1765 /* L70: */
1766  }
1767  if (*isign == -1) {
1768  cfftf(n1, (double*) &work[1], &table[table_dim1 + 1]);
1769  }
1770  if (*isign == 1) {
1771  cfftb(n1, (double*) &work[1], &table[table_dim1 + 1]);
1772  }
1773  i_3 = *n1;
1774  for (i = 1; i <= i_3; ++i) {
1775  i_4 = i + (j + k * w_dim2) * w_dim1;
1776  i_5 = i;
1777  w[i_4].r = work[i_5].r, w[i_4].i = work[i_5].i;
1778 /* L80: */
1779  }
1780 /* L90: */
1781  }
1782 /* L100: */
1783  }
1784 
1785 /* transform along Y then ... */
1786 
1787  i_1 = *n3;
1788  for (k = 1; k <= i_1; ++k) {
1789  i_2 = *n1;
1790  for (i = 1; i <= i_2; ++i) {
1791  i_3 = *n2;
1792  for (j = 1; j <= i_3; ++j) {
1793  i_4 = j;
1794  i_5 = i + (j + k * w_dim2) * w_dim1;
1795  work[i_4].r = w[i_5].r, work[i_4].i = w[i_5].i;
1796 /* L170: */
1797  }
1798  if (*isign == -1) {
1799  cfftf(n2, (double*) &work[1], &table[(table_dim1 << 1) + 1]);
1800  }
1801  if (*isign == 1) {
1802  cfftb(n2, (double*) &work[1], &table[(table_dim1 << 1) + 1]);
1803  }
1804  i_3 = *n2;
1805  for (j = 1; j <= i_3; ++j) {
1806  i_4 = i + (j + k * w_dim2) * w_dim1;
1807  i_5 = j;
1808  w[i_4].r = work[i_5].r, w[i_4].i = work[i_5].i;
1809 /* L180: */
1810  }
1811 /* L190: */
1812  }
1813 /* L200: */
1814  }
1815 
1816 /* transform along Z finally ... */
1817 
1818  i_1 = *n1;
1819  for (i = 1; i <= i_1; ++i) {
1820  i_2 = *n2;
1821  for (j = 1; j <= i_2; ++j) {
1822  i_3 = *n3;
1823  for (k = 1; k <= i_3; ++k) {
1824  i_4 = k;
1825  i_5 = i + (j + k * w_dim2) * w_dim1;
1826  work[i_4].r = w[i_5].r, work[i_4].i = w[i_5].i;
1827 /* L270: */
1828  }
1829  if (*isign == -1) {
1830  cfftf(n3, (double*) &work[1], &table[table_dim1 * 3 + 1]);
1831  }
1832  if (*isign == 1) {
1833  cfftb(n3, (double*) &work[1], &table[table_dim1 * 3 + 1]);
1834  }
1835  i_3 = *n3;
1836  for (k = 1; k <= i_3; ++k) {
1837  i_4 = i + (j + k * w_dim2) * w_dim1;
1838  i_5 = k;
1839  w[i_4].r = work[i_5].r, w[i_4].i = work[i_5].i;
1840 /* L280: */
1841  }
1842 /* L290: */
1843  }
1844 /* L300: */
1845  }
1846  return 0;
1847 } /* pubz3d_ */
1848 
1849 /****************************************************************************/
1850 
1851 int pubd3di(int n1, int n2, int n3, double *table, int ntable) {
1852 
1853  int n1over2;
1854 
1855  n1over2 = n1 / 2;
1856  return pubz3di(&n1over2,&n2,&n3,table,&ntable);
1857 
1858 } /* pubd3di */
1859 
1860 /****************************************************************************/
1861 /* real to complex fft */
1862 
1863 int pubdz3d(int isign, int n1, int n2,
1864  int n3, double *w, int ld1, int ld2, double
1865  *table, int ntable, double *work) {
1866 
1867  int n1over2, ld1over2, rval;
1868  int i, j, j2, k, k2, i1, i2, imax;
1869  double *data, *data2;
1870  double TwoPiOverN, tmp1r, tmp1i, tmp2r, tmp2i;
1871 
1872  /* complex transform */
1873  n1over2 = n1 / 2;
1874  ld1over2 = ld1 / 2;
1875  rval = pubz3d(&isign, &n1over2, &n2, &n3, (doublecomplex*) w,
1876  &ld1over2, &ld2, table, &ntable, (doublecomplex*) work);
1877 
1878  /* rearrange data */
1879  TwoPiOverN = isign * 2.0 * M_PI / n1;
1880  imax = n1/4+1;
1881  for ( i=0; i<imax; ++i ) {
1882  work[2*i] = cos(i * TwoPiOverN);
1883  work[2*i+1] = sin(i * TwoPiOverN);
1884  }
1885  for ( k=0; k<n3; ++k ) {
1886  for ( j=0; j<n2; ++j ) {
1887  data = w + ld1*(ld2*k + j);
1888  data[n1] = data[0];
1889  data[n1+1] = data[1];
1890  }
1891  }
1892  for ( k=0; k<n3; ++k ) {
1893  k2 = k?(n3-k):0;
1894  for ( j=0; j<n2; ++j ) {
1895  j2 = j?(n2-j):0;
1896  data = w + ld1*(ld2*k + j);
1897  data2 = w + ld1*(ld2*k2 + j2);
1898  imax = n1/4;
1899  if ( (n1/2) & 1 ) imax += 1;
1900  else {
1901  if ( (2*j<n2) || (2*j==n2 && 2*k<=n3) ) imax +=1;
1902  if ( j==0 && 2*k>n3 ) imax -=1;
1903  }
1904  for ( i=0; i<imax; ++i ) {
1905  i1 = 2*i; i2 = n1-i1;
1906  tmp1r = data[i1] - data2[i2];
1907  tmp1i = data[i1+1] + data2[i2+1];
1908  tmp2r = tmp1r * work[i1+1] + tmp1i * work[i1];
1909  tmp2i = tmp1i * work[i1+1] - tmp1r * work[i1];
1910  tmp1r = data[i1] + data2[i2];
1911  tmp1i = data[i1+1] - data2[i2+1];
1912  data[i1] = 0.5 * ( tmp1r + tmp2r );
1913  data[i1+1] = 0.5 * ( tmp1i + tmp2i );
1914  data2[i2] = 0.5 * ( tmp1r - tmp2r );
1915  data2[i2+1] = 0.5 * ( tmp2i - tmp1i );
1916  }
1917  }
1918  }
1919 
1920  return rval;
1921 
1922 } /* pubdz3d */
1923 
1924 /****************************************************************************/
1925 /* complex to real fft */
1926 
1927 int pubzd3d(int isign, int n1, int n2,
1928  int n3, double *w, int ld1, int ld2, double
1929  *table, int ntable, double *work) {
1930 
1931  int n1over2, ld1over2;
1932  int i, j, j2, k, k2, i1, i2, imax;
1933  double *data, *data2;
1934  double TwoPiOverN, tmp1r, tmp1i, tmp2r, tmp2i;
1935 
1936  /* rearrange data */
1937  TwoPiOverN = isign * 2.0 * M_PI / n1;
1938  imax = n1/4+1;
1939  for ( i=0; i<imax; ++i ) {
1940  work[2*i] = -cos(i * TwoPiOverN);
1941  work[2*i+1] = -sin(i * TwoPiOverN);
1942  }
1943  for ( k=0; k<n3; ++k ) {
1944  k2 = k?(n3-k):0;
1945  for ( j=0; j<n2; ++j ) {
1946  j2 = j?(n2-j):0;
1947  data = w + ld1*(ld2*k + j);
1948  data2 = w + ld1*(ld2*k2 + j2);
1949  imax = n1/4;
1950  if ( (n1/2) & 1 ) imax += 1;
1951  else {
1952  if ( (2*j<n2) || (2*j==n2 && 2*k<=n3) ) imax +=1;
1953  if ( j==0 && 2*k>n3 ) imax -=1;
1954  }
1955  for ( i=0; i<imax; ++i ) {
1956  i1 = 2*i; i2 = n1-i1;
1957  tmp1r = data[i1] - data2[i2];
1958  tmp1i = data[i1+1] + data2[i2+1];
1959  tmp2r = tmp1r * work[i1+1] + tmp1i * work[i1];
1960  tmp2i = tmp1i * work[i1+1] - tmp1r * work[i1];
1961  tmp1r = data[i1] + data2[i2];
1962  tmp1i = data[i1+1] - data2[i2+1];
1963  data[i1] = tmp1r + tmp2r;
1964  data[i1+1] = tmp1i + tmp2i;
1965  data2[i2] = tmp1r - tmp2r;
1966  data2[i2+1] = tmp2i - tmp1i;
1967  }
1968  }
1969  }
1970 
1971  /* complex transform */
1972  n1over2 = n1 / 2;
1973  ld1over2 = ld1 / 2;
1974  return pubz3d(&isign, &n1over2, &n2, &n3, (doublecomplex*) w,
1975  &ld1over2, &ld2, table, &ntable, (doublecomplex*) work);
1976 
1977 } /* pubzd3d */
1978 
1979 /****************************************************************************/
1980 
cfftb
int cfftb(int *n, double *c, double *wsave)
Definition: pub3dfft.C:1401
pubdz3d
int pubdz3d(int isign, int n1, int n2, int n3, double *w, int ld1, int ld2, double *table, int ntable, double *work)
Definition: pub3dfft.C:1863
passf5
int passf5(int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2, double *wa3, double *wa4)
Definition: pub3dfft.C:1139
cfftf1
int cfftf1(int *n, double *c, double *ch, double *wa, int *ifac)
Definition: pub3dfft.C:1420
passf3
int passf3(int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2)
Definition: pub3dfft.C:946
doublecomplex::r
double r
Definition: pub3dfft.h:10
passf
int passf(int *nac, int *ido, int *ip, int *l1, int *idl1, double *cc, double *c1, double *c2, double *ch, double *ch2, double *wa)
Definition: pub3dfft.C:647
cffti1
int cffti1(int *n, double *wa, int *ifac)
Definition: pub3dfft.C:1563
passf4
int passf4(int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2, double *wa3)
Definition: pub3dfft.C:1035
pubz3di
int pubz3di(int *n1, int *n2, int *n3, double *table, int *ntable)
Definition: pub3dfft.C:1706
M_PI
#define M_PI
Definition: pub3dfft.C:15
passb4
int passb4(int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2, double *wa3)
Definition: pub3dfft.C:408
cfftf
int cfftf(int *n, double *c, double *wsave)
Definition: pub3dfft.C:1544
doublecomplex::i
double i
Definition: pub3dfft.h:10
pubzd3d
int pubzd3d(int isign, int n1, int n2, int n3, double *w, int ld1, int ld2, double *table, int ntable, double *work)
Definition: pub3dfft.C:1927
cfftb1
int cfftb1(int *n, double *c, double *ch, double *wa, int *ifac)
Definition: pub3dfft.C:1277
passf2
int passf2(int *ido, int *l1, double *cc, double *ch, double *wa1)
Definition: pub3dfft.C:885
passb2
int passb2(int *ido, int *l1, double *cc, double *ch, double *wa1)
Definition: pub3dfft.C:257
pubd3di
int pubd3di(int n1, int n2, int n3, double *table, int ntable)
Definition: pub3dfft.C:1851
cffti
int cffti(int *n, double *wsave)
Definition: pub3dfft.C:1676
pubz3d
int pubz3d(int *isign, int *n1, int *n2, int *n3, doublecomplex *w, int *ld1, int *ld2, double *table, int *ntable, doublecomplex *work)
Definition: pub3dfft.C:1729
passb3
int passb3(int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2)
Definition: pub3dfft.C:319
passb5
int passb5(int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2, double *wa3, double *wa4)
Definition: pub3dfft.C:512
passb
int passb(int *nac, int *ido, int *ip, int *l1, int *idl1, double *cc, double *c1, double *c2, double *ch, double *ch2, double *wa)
Definition: pub3dfft.C:18
Generated by   doxygen 1.8.14