casacore
SCSL.h
Go to the documentation of this file.
1 //# extern_fft.h: C++ wrapper functions for FORTRAN FFT code
2 //# Copyright (C) 1993,1994,1995,1997,1999,2000
3 //# Associated Universities, Inc. Washington DC, USA.
4 //#
5 //# This library is free software; you can redistribute it and/or modify it
6 //# under the terms of the GNU Library General Public License as published by
7 //# the Free Software Foundation; either version 2 of the License, or (at your
8 //# option) any later version.
9 //#
10 //# This library is distributed in the hope that it will be useful, but WITHOUT
11 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 //# License for more details.
14 //#
15 //# You should have received a copy of the GNU Library General Public License
16 //# along with this library; if not, write to the Free Software Foundation,
17 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 //#
19 //# Correspondence concerning AIPS++ should be addressed as follows:
20 //# Internet email: aips2-request@nrao.edu.
21 //# Postal address: AIPS++ Project Office
22 //# National Radio Astronomy Observatory
23 //# 520 Edgemont Road
24 //# Charlottesville, VA 22903-2475 USA
25 //#
26 //# $Id$
27 
28 #ifndef SCIMATH_SCSL_H
29 #define SCIMATH_SCSL_H
30 
31 #include <casacore/casa/aips.h>
32 #include <casacore/casa/BasicSL/Complex.h>
33 
34 
35 namespace casacore { //# NAMESPACE CASACORE - BEGIN
36 
37 // <summary>C++ Interface to the Sgi/Cray Scientific Library (SCSL)</summary>
38 // <synopsis>
39 // These are C++ wrapper functions for the transform routines in the SGI/Cray
40 // Scientific Library (SCSL). The purpose of these definitions is to overload
41 // the functions so that C++ users can access the functions in SCSL with
42 // identical function names.
43 //
44 // <note role=warning>
45 // Currently, the SCSL is available only on SGI machines.
46 // </note>
47 // </synopsis>
48 
49 class SCSL
50 {
51 public:
52 // These routines compute the Fast Fourier Transform (FFT) of the complex
53 // vector x, and store the result in vector y. <src>ccfft</src> does the
54 // complex-to-complex transform and <src>zzfft</src> does the same for double
55 // precision arrays.
56 //
57 // In FFT applications, it is customary to use zero-based subscripts; the
58 // formulas are simpler that way. Suppose that the arrays are
59 // dimensioned as follows:
60 //
61 // <srcblock>
62 // COMPLEX X(0:N-1), Y(0:N-1)
63 // </srcblock>
64 //
65 // The output array is the FFT of the input array, using the following
66 // formula for the FFT:
67 //
68 // <srcblock>
69 // n-1
70 // Y(k) = scale * Sum [ X(j)*w**(isign*j*k) ] for k = 0, ..., n-1
71 // j=0
72 //
73 // where:
74 // w = exp(2*pi*i/n),
75 // i = + sqrt(-1),
76 // pi = 3.14159...,
77 // isign = +1 or -1
78 // </srcblock>
79 //
80 // Different authors use different conventions for which of the
81 // transforms, isign = +1 or isign = -1, is the forward or inverse
82 // transform, and what the scale factor should be in either case. You
83 // can make this routine compute any of the various possible definitions,
84 // however, by choosing the appropriate values for isign and scale.
85 //
86 // The relevant fact from FFT theory is this: If you take the FFT with
87 // any particular values of isign and scale, the mathematical inverse
88 // function is computed by taking the FFT with -isign and 1/(n*scale).
89 // In particular, if you use isign = +1 and scale = 1.0 you can compute
90 // the inverse FFT by using isign = -1 and scale = 1.0/n.
91 //
92 // The output array may be the same as the input array.
93 //
94 // <h3>Initialization</h3>
95 // The table array stores the trigonometric tables used in calculation of
96 // the FFT. You must initialize table by calling the routine with isign
97 // = 0 prior to doing the transforms. If the value of the problem size,
98 // n, does not change, table does not have to be reinitialized.
99 //
100 // <h3>Dimensions</h3>
101 // In the preceding description, it is assumed that array subscripts were
102 // zero-based, as is customary in FFT applications. Thus, the input and
103 // output arrays are declared as follows:
104 //
105 // <srcblock>
106 // COMPLEX X(0:N-1)
107 // COMPLEX Y(0:N-1)
108 // </srcblock>
109 //
110 // However, if you prefer to use the more customary FORTRAN style with
111 // subscripts starting at 1 you do not have to change the calling
112 // sequence, as in the following (assuming N > 0):
113 //
114 // <srcblock>
115 // COMPLEX X(N)
116 // COMPLEX Y(N)
117 // </srcblock>
118 //
119 // <example>
120 // These examples use the table and workspace sizes appropriate to the
121 // Origin series.
122 //
123 // Example 1: Initialize the complex array table in preparation for
124 // doing an FFT of size 1024. Only the isign, n, and table arguments are
125 // used in this case. You can use dummy arguments or zeros for the other
126 // arguments in the subroutine call.
127 //
128 // <srcblock>
129 // REAL TABLE(30 + 2048)
130 // CALL CCFFT(0, 1024, 0.0, DUMMY, DUMMY, TABLE, DUMMY, 0)
131 // </srcblock>
132 //
133 // Example 2: x and y are complex arrays of dimension (0:1023). Take
134 // the FFT of x and store the results in y. Before taking the FFT,
135 // initialize the table array, as in example 1.
136 //
137 // <srcblock>
138 // COMPLEX X(0:1023), Y(0:1023)
139 // REAL TABLE(30 + 2048)
140 // REAL WORK(2048)
141 // CALL CCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
142 // CALL CCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)
143 // </srcblock>
144 //
145 // Example 3: Using the same x and y as in example 2, take the inverse
146 // FFT of y and store it back in x. The scale factor 1/1024 is used.
147 // Assume that the table array is already initialized.
148 //
149 // <srcblock>
150 // CALL CCFFT(-1, 1024, 1.0/1024.0, Y, X, TABLE, WORK, 0)
151 // </srcblock>
152 //
153 // Example 4: Perform the same computation as in example 2, but assume
154 // that the lower bound of each array is 1, rather than 0. No change was
155 // needed in the subroutine calls.
156 //
157 // <srcblock>
158 // COMPLEX X(1024), Y(1024)
159 // CALL CCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
160 // CALL CCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)
161 // </srcblock>
162 //
163 // Example 5: Do the same computation as in example 4, but put the
164 // output back in array x to save storage space. Assume that table is
165 // already initialized.
166 //
167 // <srcblock>
168 // COMPLEX X(1024)
169 // CALL CCFFT(1, 1024, 1.0, X, X, TABLE, WORK, 0)
170 // </srcblock>
171 // </example>
172 //
173 // Input parameters:
174 // <dl compact>
175 // <dt><b>isign</b>
176 // <dd> Integer.
177 // Specifies whether to initialize the table array or to do the
178 // forward or inverse Fourier transform, as follows:
179 //
180 // If isign = 0, the routine initializes the table array and
181 // returns. In this case, the only arguments used or checked
182 // are isign, n, and table.
183 //
184 // If isign = +1 or -1, the value of isign is the sign of the
185 // exponent used in the FFT formula.
186 // <dt><b>n</b>
187 // <dd> Integer. Size of the transform (the number of values in
188 // the input array). n >= 1.
189 // <dt><b>scale</b>
190 // <dd> Scale factor.
191 // <src>ccfft</src>: real.
192 // <src>zzfft</src>: double precision.
193 // Each element of the output array is multiplied by scale
194 // after taking the Fourier transform, as defined by the previous
195 // formula.
196 // <dt><b>x</b>
197 // <dd> Array of dimension (0:n-1).
198 // <src>ccfft</src>: complex array.
199 // <src>zzfft</src>: double complex array.
200 //
201 // Input array of values to be transformed.
202 // <dt><b>isys</b>
203 // <dd> Integer.
204 // Algorithm used; value dependent on hardware system. Currently, no
205 // special options are supported; therefore, you must always specify
206 // an isys argument as constant 0.
207 // </dl>
208 // Output parameters:
209 // <dl compact>
210 // <dt><b>y</b>
211 // <dd> Array of dimension (0:n-1).
212 // <src>ccfft</src>: complex array.
213 // <src>zzfft</src>: double complex array.
214 // Output array of transformed values. The output array may be
215 // the same as the input array. In that case, the transform is
216 // done in place and the input array is overwritten with the
217 // transformed values.
218 // <dt><b>table</b>
219 // <dd> Real array; dimension 2*n+30.
220 //
221 // Table of factors and trigonometric functions.
222 //
223 // If isign = 0, the routine initializes table (table is output
224 // only).
225 //
226 // If isign = +1 or -1, the values in table are assumed to be
227 // initialized already by a prior call with isign = 0 (table is
228 // input only).
229 // <dt><b>work</b>
230 // <dd> Real array; dimension 2*n.
231 //
232 // Work array. This is a scratch array used for intermediate
233 // calculations. Its address space must be different address
234 // space from that of the input and output arrays.
235 // </dl>
236 // <group>
237 static void ccfft(Int isign, Int n, Float scale, Complex* x,
238  Complex* y, Float* table, Float* work, Int isys);
239 static void ccfft(Int isign, Int n, Double scale, DComplex* x,
240  DComplex* y, Double* table, Double* work, Int isys);
241 static void zzfft(Int isign, Int n, Double scale, DComplex* x,
242  DComplex* y, Double* table, Double* work, Int isys);
243 // </group>
244 
245 // <src>scfft/dzfft</src> computes the FFT of the real array x, and it stores
246 // the results in the complex array y. <src>csfft/zdfft</src> computes the
247 // corresponding inverse complex-to-real transform.
248 //
249 // It is customary in FFT applications to use zero-based subscripts; the
250 // formulas are simpler that way. For these routines, suppose that the
251 // arrays are dimensioned as follows:
252 //
253 // <srcblock>
254 // REAL X(0:n-1)
255 // COMPLEX Y(0:n/2)
256 // </srcblock>
257 //
258 // Then the output array is the FFT of the input array, using the
259 // following formula for the FFT:
260 //
261 // <srcblock>
262 // n-1
263 // Y(k) = scale * Sum [ X(j)*w**(isign*j*k) ] for k = 0, ..., n/2
264 // j=0
265 //
266 // where:
267 // w = exp(2*pi*i/n),
268 // i = + sqrt(-1),
269 // pi = 3.14159...,
270 // isign = +1 or -1.
271 // </srcblock>
272 //
273 // Different authors use different conventions for which of the
274 // transforms, isign = +1 or isign = -1, is the forward or inverse
275 // transform, and what the scale factor should be in either case. You
276 // can make these routines compute any of the various possible
277 // definitions, however, by choosing the appropriate values for isign and
278 // scale.
279 //
280 // The relevant fact from FFT theory is this: If you call <src>scfft</src>
281 // with any particular values of isign and scale, the mathematical inverse
282 // function is computed by calling <src>csfft</src> with -isign and
283 // 1/(n*scale). In particular, if you use isign = +1 and scale = 1.0 in
284 // <src>scfft</src> for the forward FFT, you can compute the inverse FFT by
285 // using <src>ccfft</src> with isign = -1 and scale = 1.0/n.
286 //
287 // <h3>Real-to-complex FFTs</h3>
288 // Notice in the preceding formula that there are n real input values,
289 // and n/2 + 1 complex output values. This property is characteristic of
290 // real-to-complex FFTs.
291 //
292 // The mathematical definition of the Fourier transform takes a sequence
293 // of n complex values and transforms it to another sequence of n complex
294 // values. A complex-to-complex FFT routine, such as <src>ccfft</src>, will
295 // take n complex input values, and produce n complex output values. In
296 // fact, one easy way to compute a real-to-complex FFT is to store the
297 // input data in a complex array, then call routine <src>ccfft</src> to
298 // compute the FFT. You get the same answer when using the <src>scfft</src>
299 // routine.
300 //
301 // The reason for having a separate real-to-complex FFT routine is
302 // efficiency. Because the input data is real, you can make use of this
303 // fact to save almost half of the computational work. The theory of
304 // Fourier transforms tells us that for real input data, you have to
305 // compute only the first n/2 + 1 complex output values, because the
306 // remaining values can be computed from the first half of the values by
307 // the simple formula:
308 //
309 // <srcblock>
310 // Y(k) = conjg(Y(n-k)) for n/2 <= k <= n-1
311 // </srcblock>
312 //
313 // where the notation conjgY represents the complex conjugate of y.
314 //
315 // In fact, in many applications, the second half of the complex output
316 // data is never explicitly computed or stored. Likewise, as explained
317 // later, only the first half of the complex data has to be supplied for
318 // the complex-to-real FFT.
319 //
320 // Another implication of FFT theory is that, for real input data, the
321 // first output value, Y(0), will always be a real number; therefore, the
322 // imaginary part will always be 0. If n is an even number, Y(n/2) will
323 // also be real and thus, have zero imaginary parts.
324 //
325 // <h3>Complex-to-real FFTs</h3>
326 // Consider the complex-to-real case. The effect of the computation is
327 // given by the preceding formula, but with X complex and Y real.
328 //
329 // Generally, the FFT transforms a complex sequence into a complex
330 // sequence. However, in a certain application we may know the output
331 // sequence is real. Often, this is the case because the complex input
332 // sequence was the transform of a real sequence. In this case, you can
333 // save about half of the computational work.
334 //
335 // According to the theory of Fourier transforms, for the output
336 // sequence, Y, to be a real sequence, the following identity on the
337 // input sequence, X, must be true:
338 //
339 // <srcblock>
340 // X(k) = conjg(X(n-k)) for n/2 <= k <= n-1
341 // </srcblock>
342 //
343 // And, in fact, the input values X(k) for k > n/2 need not be supplied;
344 // they can be inferred from the first half of the input.
345 //
346 // Thus, in the complex-to-real routine, <src>csfft</src>, the arrays can be
347 // dimensioned as follows:
348 //
349 // <srcblock>
350 // COMPLEX X(0:n/2)
351 // REAL Y(0:n-1)
352 // </srcblock>
353 //
354 // There are n/2 + 1 complex input values and n real output values. Even
355 // though only n/2 + 1 input values are supplied, the size of the
356 // transform is still n in this case, because implicitly you are using
357 // the FFT formula for a sequence of length n.
358 //
359 // Another implication of the theory is that X(0) must be a real number
360 // (that is, it must have zero imaginary part). Also, if n is even,
361 // X(n/2) must also be real. Routine <src>CSFFT</src> assumes that these
362 // values are real; if you specify a nonzero imaginary part, it is ignored.
363 //
364 // <h3>Table Initialization</h3>
365 // The table array stores the trigonometric tables used in calculation of
366 // the FFT. This table must be initialized by calling the routine with
367 // isign = 0 prior to doing the transforms. The table does not have to
368 // be reinitialized if the value of the problem size, n, does not change.
369 // Because <src>scfft</src> and <src>csfft</src> use the same format for
370 // table, either can be used to initialize it (note that CCFFT uses a
371 // different table format).
372 //
373 // <h3>Dimensions</h3>
374 // In the preceding description, it is assumed that array subscripts were
375 // zero-based, as is customary in FFT applications. Thus, the input and
376 // output arrays are declared (assuming n > 0):
377 //
378 // <srcblock>
379 // REAL X(0:n-1)
380 // COMPLEX Y(0:n/2)
381 // </srcblock>
382 //
383 // No change is needed in the calling sequence; however, if you prefer
384 // you can use the more customary Fortran style with subscripts starting
385 // at 1, as in the following:
386 //
387 // <srcblock>
388 // REAL X(n)
389 // COMPLEX Y(n/2 + 1)
390 // </srcblock>
391 //
392 // <example>
393 // These examples use the table and workspace sizes appropriate to Origin
394 // series.
395 //
396 // Example 1: Initialize the complex array TABLE in preparation for
397 // doing an FFT of size 1024. In this case only the arguments isign, n,
398 // and table are used. You can use dummy arguments or zeros for the other
399 // arguments in the subroutine call.
400 //
401 // <srcblock>
402 // REAL TABLE(15 + 1024)
403 // CALL SCFFT(0, 1024, 0.0, DUMMY, DUMMY, TABLE, DUMMY, 0)
404 // </srcblock>
405 //
406 // Example 2: X is a real array of dimension (0:1023), and Y is a
407 // complex array of dimension (0:512). Take the FFT of X and store the
408 // results in Y. Before taking the FFT, initialize the TABLE array, as
409 // in example 1.
410 //
411 // <srcblock>
412 // REAL X(0:1023)
413 // COMPLEX Y(0:512)
414 // REAL TABLE(15 + 1024)
415 // REAL WORK(1024)
416 // CALL SCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
417 // CALL SCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)
418 // </srcblock>
419 //
420 // Example 3: With X and Y as in example 2, take the inverse FFT of Y
421 // and store it back in X. The scale factor 1/1024 is used. Assume that
422 // the TABLE array is initialized already.
423 //
424 // <srcblock>
425 // CALL CSFFT(-1, 1024, 1.0/1024.0, Y, X, TABLE, WORK, 0)
426 // </srcblock>
427 //
428 // Example 4: Perform the same computation as in example 2, but assume
429 // that the lower bound of each array is 1, rather than 0. The
430 // subroutine calls are not changed.
431 //
432 // <srcblock>
433 // REAL X(1024)
434 // COMPLEX Y(513)
435 // CALL SCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
436 // CALL SCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)
437 // </srcblock>
438 //
439 // Example 5: Perform the same computation as in example 4, but
440 // equivalence the input and output arrays to save storage space. Assume
441 // that the TABLE array is initialized already.
442 //
443 // <srcblock>
444 // REAL X(1024)
445 // COMPLEX Y(513)
446 // EQUIVALENCE ( X(1), Y(1) )
447 // CALL SCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)
448 // </srcblock>
449 // </example>
450 //
451 // Input parameters:
452 // <dl compact>
453 // <dt><b>isign</b>
454 // <dd> Integer.
455 // Specifies whether to initialize the table array or to do the
456 // forward or inverse Fourier transform, as follows:
457 //
458 // If isign = 0, the routine initializes the table array and
459 // returns. In this case, the only arguments used or checked
460 // are isign, n, and table.
461 //
462 // If isign = +1 or -1, the value of isign is the sign of the
463 // exponent used in the FFT formula.
464 // <dt><b>n</b>
465 // <dd> Integer.
466 // Size of transform. If n <= 2, <src>scfft/dzfft</src>
467 // returns without calculating the transform.
468 // <dt><b>scale</b>
469 // <dd> Scale factor.
470 // <src>scfft</src>: real.
471 // <src>dzfft</src>: double precision.
472 // <src>csfft</src>: real.
473 // <src>zdfft</src>: double precision.
474 // Each element of the output array is multiplied by scale
475 // after taking the Fourier transform, as defined by the previous
476 // formula.
477 // <dt><b>x</b>
478 // <dd> Input array of values to be transformed.
479 // <src>scfft</src>: real array of dimension (0:n-1).
480 // <src>dzfft</src>: double precision array of dimension (0:n-1).
481 // <src>csfft</src>: complex array of dimension (0:n/2).
482 // <src>zdfft</src>: double complex array of dimension (0:n/2).
483 // <dt><b>isys</b>
484 // <dd> Integer array of dimension (0:isys(0)).
485 // Use isys to specify certain processor-specific parameters or
486 // options. The first element of the array specifies how many
487 // more elements are in the array.
488 //
489 // If isys(0) = 0, the default values of such parameters are
490 // used. In this case, you can specify the argument value as
491 // the scalar integer constant 0. If isys(0) > 0, isys(0)
492 // gives the upper bound of the isys array; that is, if
493 // il = isys(0), user-specified parameters are expected in
494 // isys(1) through isys(il).
495 // </dl>
496 // Output parameters:
497 // <dl compact>
498 // <dt><b>y</b>
499 // <dd> Output array of transformed values.
500 // <src>scfft</src>: complex array of dimension (0:n/2).
501 // <src>dzfft</src>: double complex array of dimension (0:n/2).
502 // <src>csfft</src>: real array of dimension (0:n-1).
503 // <src>zdfft</src>: double precision array of dimension (0:n-1).
504 //
505 // The output array, y, is the FFT of the the input array, x,
506 // computed according to the preceding formula. The output
507 // array may be equivalenced to the input array in the calling
508 // program. Be careful when dimensioning the arrays, in this
509 // case, to allow for the fact that the complex array contains
510 // two (real) words more than the real array.
511 // <dt><b>table</b>
512 // <dd> Real array; dimension n+15.
513 //
514 // Table of factors and trigonometric functions.
515 //
516 // If isign = 0, the table array is initialized to contain
517 // trigonometric tables needed to compute an FFT of size n.
518 //
519 // If isign = +1 or -1, the values in table are assumed to be
520 // initialized already by a prior call with isign = 0.
521 // <dt><b>work</b>
522 // <dd> Real array; dimension n.
523 //
524 // Work array used for intermediate calculations. Its address
525 // space must be different from that of the input and output
526 // arrays.
527 // </dl>
528 // <group>
529 static void scfft(Int isign, Int n, Float scale, Float* x,
530  Complex* y, Float* table, Float* work, Int isys);
531 static void scfft(Int isign, Int n, Double scale, Double* x,
532  DComplex* y, Double* table, Double* work, Int isys);
533 static void dzfft(Int isign, Int n, Double scale, Double* x,
534  DComplex* y, Double* table, Double* work, Int isys);
535 static void csfft(Int isign, Int n, Float scale, Complex* x,
536  Float* y, Float* table, Float* work, Int isys);
537 static void csfft(Int isign, Int n, Double scale, DComplex* x,
538  Double* y, Double* table, Double* work, Int isys);
539 static void zdfft(Int isign, Int n, Double scale, DComplex* x,
540  Double* y, Double* table, Double* work, Int isys);
541 // </group>
542 
543 // <src>ccfftm/zzfftm</src> computes the FFT of each column of the
544 // complex matrix x, and stores the results in the columns of complex
545 // matrix y.
546 //
547 // Suppose the arrays are dimensioned as follows:
548 //
549 // <srcblock>
550 // COMPLEX X(0:ldx-1, 0:lot-1)
551 // COMPLEX Y(0:ldy-1, 0:lot-1)
552 //
553 // where ldx >= n, ldy >= n.
554 // </srcblock>
555 //
556 // Then column L of the output array is the FFT of column L of the
557 // input array, using the following formula for the FFT:
558 //
559 // <srcblock>
560 // n-1
561 // Y(k, L) = scale * Sum [ X(j)*w**(isign*j*k) ]
562 // j=0
563 // for k = 0, ..., n-1
564 // L = 0, ..., lot-1
565 // where:
566 // w = exp(2*pi*i/n),
567 // i = + sqrt(-1),
568 // pi = 3.14159...,
569 // isign = +1 or -1
570 // lot = the number of columns to transform
571 // </srcblock>
572 //
573 // Different authors use different conventions for which of the
574 // transforms, isign = +1 or isign = -1, is the forward or inverse
575 // transform, and what the scale factor should be in either case. You
576 // can make this routine compute any of the various possible definitions,
577 // however, by choosing the appropriate values for isign and scale.
578 //
579 // The relevant fact from FFT theory is this: If you take the FFT with
580 // any particular values of isign and scale, the mathematical inverse
581 // function is computed by taking the FFT with -isign and 1/(n * scale).
582 // In particular, if you use isign = +1 and scale = 1.0 for the forward
583 // FFT, you can compute the inverse FFT by using the following: isign =
584 // -1 and scale = 1.0/n.
585 //
586 // This section contains information about the algorithm for these
587 // routines, the initialization of the table array, the declaration of
588 // dimensions for x and y arrays, some performance tips, and some
589 // implementation-dependent details.
590 //
591 // <h3>Algorithm</h3>
592 // These routines use decimation-in-frequency type FFT. It takes the FFT
593 // of the columns and vectorizes the operations along the rows of the
594 // matrix. Thus, the vector length in the calculations depends on the
595 // row size, and the strides for vector loads and stores are the leading
596 // dimensions, ldx and ldy.
597 //
598 // <h3>Initialization</h3>
599 // The table array stores the trigonometric tables used in calculation of
600 // the FFT. You must initialize the table array by calling the routine
601 // with isign = 0 prior to doing the transforms. If the value of the
602 // problem size, n, does not change, table does not have to be
603 // reinitialized.
604 //
605 // <h3>Dimensions</h3>
606 // In the preceding description, it is assumed that array subscripts were
607 // zero-based, as is customary in FFT applications. Thus, the input and
608 // output arrays are declared as follows:
609 //
610 // <srcblock>
611 // COMPLEX X(0:ldx-1, 0:lot-1)
612 // COMPLEX Y(0:ldy-1, 0:lot-1)
613 // </srcblock>
614 //
615 // The calling sequence does not have to change, however, if you prefer
616 // to use the more customary Fortran style with subscripts starting at 1.
617 // The same values of ldx and ldy would be passed to the subroutine even
618 // if the input and output arrays were dimensioned as follows:
619 //
620 // <srcblock>
621 // COMPLEX X(ldx, lot)
622 // COMPLEX Y(ldy, lot)
623 // </srcblock>
624 //
625 // <example>
626 // Example 1: Initialize the TABLE array in preparation for doing an FFT
627 // of size 128. Only the isign, n, and table arguments are used in this
628 // case. You can use dummy arguments or zeros for the other arguments in
629 // the subroutine call.
630 //
631 // <srcblock>
632 // REAL TABLE(30 + 256)
633 // CALL CCFFTM(0, 128, 0, 0., DUMMY, 1, DUMMY, 1, TABLE, DUMMY, 0)
634 // </srcblock>
635 //
636 // Example 2: X and Y are complex arrays of dimension (0:128) by (0:55).
637 // The first 128 elements of each column contain data. For performance
638 // reasons, the extra element forces the leading dimension to be an odd
639 // number. Take the FFT of the first 50 columns of X and store the
640 // results in the first 50 columns of Y. Before taking the FFT,
641 // initialize the TABLE array, as in example 1.
642 //
643 // <srcblock>
644 // COMPLEX X(0:128, 0:55)
645 // COMPLEX Y(0:128, 0:55)
646 // REAL TABLE(30 + 256)
647 // REAL WORK(256)
648 // ...
649 // CALL CCFFTM(0, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
650 // CALL CCFFTM(1, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
651 // </srcblock>
652 //
653 // Example 3: With X and Y as in example 2, take the inverse FFT of Y
654 // and store it back in X. The scale factor 1/128 is used. Assume that
655 // the TABLE array is already initialized.
656 //
657 // <srcblock>
658 // CALL CCFFTM(-1, 128, 50, 1./128., Y, 129, X, 129, TABLE,WORK,0)
659 // </srcblock>
660 //
661 // Example 4: Perform the same computation as in example 2, but assume
662 // that the lower bound of each array is 1, rather than 0. The
663 // subroutine calls are not changed.
664 //
665 // <srcblock>
666 // COMPLEX X(129, 55)
667 // COMPLEX Y(129, 55)
668 // ...
669 // CALL CCFFTM(0, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
670 // CALL CCFFTM(1, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
671 // </srcblock>
672 //
673 // Example 5: Perform the same computation as in example 4, but put the
674 // output back in array X to save storage space. Assume that the TABLE
675 // array is already initialized.
676 //
677 // <srcblock>
678 // COMPLEX X(129, 55)
679 // ...
680 // CALL CCFFTM(1, 128, 50, 1.0, X, 129, X, 129, TABLE, WORK, 0)
681 // </srcblock>
682 //
683 // </example>
684 //
685 // Input parameters:
686 // <dl compact>
687 // <dt><b>isign</b>
688 // <dd> Integer.
689 // Specifies whether to initialize the table array or to do the
690 // forward or inverse Fourier transform, as follows:
691 //
692 // If isign = 0, the routine initializes the table array and
693 // returns. In this case, the only arguments used or checked
694 // are isign, n, and table.
695 //
696 // If isign = +1 or -1, the value of isign is the sign of the
697 // exponent used in the FFT formula.
698 // <dt><b>n</b>
699 // <dd> Integer.
700 // Size of each transform (the number of elements in each
701 // column of the input and output matrix to be transformed).
702 // Performance depends on the value of n, as explained above.
703 // n >= 0; if n = 0, the routine returns.
704 // <dt><b>lot</b>
705 // <dd> Integer.
706 // The number of transforms to be computed (lot size). This is
707 // the number of elements in each row of the input and output
708 // matrix. lot >= 0. If lot = 0, the routine returns.
709 // <dt><b>scale</b>
710 // <dd> Scale factor.
711 // <src>ccfftm</src>: real.
712 // <src>zzfftm</src>: double precision.
713 // Each element of the output array is multiplied by scale
714 // factor after taking the Fourier transform, as defined
715 // previously.
716 // <dt><b>x</b>
717 // <dd> Array of dimension (0:ldx-1, 0:n2-1).
718 // <src>ccfftm</src>: real array.
719 // <src>zzfftm</src>: double precision array.
720 // Input array of values to be transformed.
721 // <dt><b>ldx</b>
722 // <dd> The number of rows in the x array, as it was declared in the
723 // calling program (the leading dimension of X). ldx >= MAX(n, 1).
724 // <dt><b>ldy</b>
725 // <dd> Integer.
726 // The number of rows in the y array, as it was declared in the
727 // calling program (the leading dimension of y). ldy >= MAX(n,
728 // 1).
729 // <dt><b>isys</b>
730 // <dd> Integer array of dimension (0:isys(0)).
731 // The first element of the array specifies how many more
732 // elements are in the array. Use isys to specify certain
733 // processor-specific parameters or options.
734 //
735 // If isys(0) = 0, the default values of such parameters are
736 // used. In this case, you can specify the argument value as
737 // the scalar integer constant 0.
738 //
739 // If isys(0) > 0, isys(0) gives the upper bound of the isys
740 // array. Therefore, if il = isys(0), user-specified
741 // parameters are expected in isys(1) through isys(il).
742 // </dl>
743 // Output parameters:
744 // <dl compact>
745 // <dt><b>y</b>
746 // <dd> Array of dimension (0:ldy-1, 0:lot-1).
747 // <src>ccfftm</src>: complex array.
748 // <src>zzfftm</src>: double complex array.
749 // Output array of transformed values. Each column of the
750 // output array, y, is the FFT of the corresponding column of
751 // the input array, x, computed according to the preceding
752 // formula.
753 //
754 // The output array may be the same as the input array. In that
755 // case, the transform is done in place. The input array is
756 // overwritten with the transformed values. In this case, it
757 // is necessary that ldx = ldy.
758 // <dt><b>table</b>
759 // <dd> Real array; dimension (30 + 2n).
760 // Table of factors and trigonometric functions.
761 //
762 // If isign = 0, the routine initializes table (table is output
763 // only).
764 //
765 // If isign = +1 or -1, the values in table are assumed to be
766 // initialized already by a prior call with isign = 0 (table is
767 // input only).
768 // <dt><b>work</b>
769 // <dd> Real array; dimension 2n.
770 // Work array. This is a scratch array used for intermediate
771 // calculations. Its address space must be different from that
772 // of the input and output arrays.
773 // </dl>
774 // <group>
775 static void ccfftm(Int isign, Int n, Int lot, Float scale, Complex*
776  x, Int ldx, Complex* y, Int ldy, Float* table,
777  Float* work, Int isys);
778 static void zzfftm(Int isign, Int n, Int lot, Double scale, DComplex*
779  x, Int ldx, DComplex* y, Int ldy, Double* table,
780  Double* work, Int isys);
781 // </group>
782 
783 // <src>scfftm/dzfftm</src> computes the FFT of each column of the real matrix
784 // X, and it stores the results in the corresponding column of the complex
785 // matrix Y. <src>csfftm/zdfftm</src> computes the corresponding inverse
786 // transforms.
787 //
788 // In FFT applications, it is customary to use zero-based subscripts; the
789 // formulas are simpler that way. First, the function of <src>scfftm</src> is
790 // described. Suppose that the arrays are dimensioned as follows:
791 //
792 // <srcblock>
793 // REAL X(0:ldx-1, 0:lot-1)
794 // COMPLEX Y(0:ldy-1, 0:lot-1)
795 //
796 // where ldx >= n, ldy >= n/2 + 1.
797 // </srcblock>
798 //
799 // Then column L of the output array is the FFT of column L of the input
800 // array, using the following formula for the FFT:
801 //
802 // <srcblock>
803 // n-1
804 // Y(k, L) = scale * Sum [ X(j, L)*w**(isign*j*k) ]
805 // j=0
806 //
807 // for k = 0, ..., n/2
808 // L = 0, ..., lot-1 where:
809 // w = exp(2*pi*i/n),
810 // i = + sqrt(-1)
811 // pi = 3.14159...,
812 // isign = +1 or -1,
813 // lot = the number of columns to transform
814 // </srcblock>
815 //
816 // Different authors use different conventions for which transform
817 // (isign = +1 or isign = -1) is used in the real-to-complex case, and
818 // what the scale factor should be. Some adopt the convention that isign
819 // = 1 for the real-to-complex transform, and isign = -1 for the
820 // complex-to-real inverse. Others use the opposite convention. You can
821 // make these routines compute any of the various possible definitions,
822 // however, by choosing the appropriate values for isign and scale.
823 //
824 // The relevant fact from FFT theory is this: If you use <src>scfftm</src> to
825 // take the real-to-complex FFT, using any particular values of isign and
826 // scale, the mathematical inverse function is computed by using
827 // <src>csfftm</src> with -isign and 1/ (n*scale). In particular, if you call
828 // <src>scfftm</src> with isign = +1 and scale = 1.0, you can use
829 // <src>csfftm</src> to compute the inverse complex-to-real FFT by using isign
830 // = -1 and scale = 1.0/n.
831 //
832 // <h3>Real-to-complex FFTs</h3>
833 // Notice in the preceding formula that there are n real input values and
834 // (n/2) + 1 complex output values for each column. This property is
835 // characteristic of real-to-complex FFTs.
836 //
837 // The mathematical definition of the Fourier transform takes a sequence
838 // of n complex values and transforms it to another sequence of n complex
839 // values. A complex-to-complex FFT routine, such as <src>ccfftm</src>, will
840 // take n complex input values and produce n complex output values. In fact,
841 // one easy way to compute a real-to-complex FFT is to store the input
842 // data x in a complex array, then call routine <src>ccfftm</src> to compute
843 // the FFT. You get the same answer when using the <src>scfftm</src> routine.
844 //
845 // A separate real-to-complex FFT routine is more efficient than the
846 // equivalent complex-to-complex routine. Because the input data is
847 // real, you can make use of this fact to save almost half of the
848 // computational work. According to the theory of Fourier transforms,
849 // for real input data, you have to compute only the first n/2 + 1
850 // complex output values in each column, because the second half of the
851 // FFT values in each column can be computed from the first half of the
852 // values by the simple formula:
853 //
854 // <srcblock>
855 // Y = conjgY for n/2 <= k <= n-1
856 // k,L n-k, L
857 //
858 // where the notation conjg(z) represents the complex conjugate of z.
859 // </srcblock>
860 //
861 // In fact, in many applications, the second half of the complex output
862 // data is never explicitly computed or stored. Likewise, you must
863 // supply only the first half of the complex data in each column has to
864 // be supplied for the complex-to-real FFT.
865 //
866 // Another implication of FFT theory is that for real input data, the
867 // first output value in each column, Y(0, L), will always be a real
868 // number; therefore, the imaginary part will always be 0. If n is an
869 // even number, Y(n/2, L) will also be real and have 0 imaginary parts.
870 //
871 // <h3>Complex-to-real FFTs</h3>
872 // Consider the complex-to-real case. The effect of the computation is
873 // given by the preceding formula, but with X complex and Y real.
874 //
875 // In general, the FFT transforms a complex sequence into a complex
876 // sequence; however, in a certain application you may know the output
877 // sequence is real, perhaps because the complex input sequence was the
878 // transform of a real sequence. In this case, you can save about half
879 // of the computational work.
880 //
881 // According to the theory of Fourier transforms, for the output
882 // sequence, Y, to be a real sequence, the following identity on the
883 // input sequence, X, must be true:
884 //
885 // <srcblock>
886 // X = conjgX for n/2 <= k <= n-1
887 // k,L n-k,L
888 // And, in fact, the following input values
889 //
890 // X for k > n/2
891 // k,L
892 // do not have to be supplied, because they can be inferred from the
893 // first half of the input.
894 // </srcblock>
895 //
896 // Thus, in the complex-to-real routine, CSFFTM, the arrays can be
897 // dimensioned as follows:
898 //
899 // <srcblock>
900 // COMPLEX X(0:ldx-1, 0:lot-1)
901 // REAL Y(0:ldy-1, 0:lot-1)
902 //
903 // where ldx >= n/2 + 1, ldy >= n.
904 // </srcblock>
905 //
906 // In each column, there are (n/2) + 1 complex input values and n real
907 // output values. Even though only (n/2) + 1 input values are supplied,
908 // the size of the transform is still n in this case, because implicitly
909 // the FFT formula for a sequence of length n is used.
910 //
911 // Another implication of the theory is that X(0, L) must be a real
912 // number (that is, must have zero imaginary part). If n is an even
913 // number, X(n/2, L) must also be real. Routine CSFFTM assumes that each
914 // of these values is real; if a nonzero imaginary part is given, it is
915 // ignored.
916 //
917 // <h3>Table Initialization</h3>
918 // The table array contains the trigonometric tables used in calculation
919 // of the FFT. You must initialize this table by calling the routine
920 // with isign = 0 prior to doing the transforms. table does not have to
921 // be reinitialized if the value of the problem size, n, does not change.
922 //
923 // <h3>Dimensions</h3>
924 // In the preceding description, it is assumed that array subscripts were
925 // zero-based, as is customary in FFT applications. Thus, the input and
926 // output arrays are declared (for SCFFTM):
927 //
928 // <srcblock>
929 // REAL X(0:ldx-1, 0:lot-1)
930 // COMPLEX Y(0:ldy-1, 0:lot-1)
931 // </srcblock>
932 //
933 // No change is made in the calling sequence, however, if you prefer to
934 // use the more customary Fortran style with subscripts starting at 1.
935 // The same values of ldx and ldy would be passed to the subroutine even
936 // if the input and output arrays were dimensioned as follows:
937 //
938 // <srcblock>
939 // REAL X(ldx, lot)
940 // COMPLEX Y(ldy, lot)
941 // </srcblock>
942 
943 // </example>
944 // Example 1: Initialize the complex array TABLE in preparation for
945 // doing an FFT of size 128. In this case only the isign, n, and table
946 // arguments are used; you may use dummy arguments or zeros for the other
947 // arguments in the subroutine call.
948 //
949 // <srcblock>
950 // REAL TABLE(15 + 128)
951 // CALL SCFFTM(0, 128, 1, 0.0, DUMMY, 1, DUMMY, 1,
952 // & TABLE, DUMMY, 0)
953 // </srcblock>
954 //
955 // Example 2: X is a real array of dimension (0:128, 0:55), and Y is a
956 // complex array of dimension (0:64, 0:55). The first 128 elements in
957 // each column of X contain data; the extra element forces an odd leading
958 // dimension. Take the FFT of the first 50 columns of X and store the
959 // results in the first 50 columns of Y. Before taking the FFT,
960 // initialize the TABLE array, as in example 1.
961 //
962 // <srcblock>
963 // REAL X(0:128, 0:55)
964 // COMPLEX Y(0:64, 0:55)
965 // REAL TABLE(15 + 128)
966 // REAL WORK((128)
967 // ...
968 // CALL SCFFTM(0, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
969 // CALL SCFFTM(1, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
970 // </srcblock>
971 //
972 // Example 3: With X and Y as in example 2, take the inverse FFT of Y
973 // and store it back in X. The scale factor 1/128 is used. Assume that
974 // the TABLE array is initialized already.
975 //
976 // <srcblock>
977 // CALL CSFFTM(-1, 128, 50, 1.0/128.0, Y, 65, X, 129,
978 // & TABLE, WORK, 0)
979 // </srcblock>
980 //
981 // Example 4: Perform the same computation as in example 2, but assume
982 // that the lower bound of each array is 1, rather than 0. No change is
983 // made in the subroutine calls.
984 //
985 // <srcblock>
986 // REAL X(129, 56)
987 // COMPLEX Y(65, 56)
988 // ...
989 // CALL SCFFTM(0, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
990 // CALL SCFFTM(1, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
991 // </srcblock>
992 //
993 // Example 5: Perform the same computation as in example 4, but
994 // equivalence the input and output arrays to save storage space. In
995 // this case, a row must be added to X, because it is equivalenced to a
996 // complex array. The leading dimension of X is two times an odd number;
997 // therefore, memory bank conflicts are minimal. Assume that TABLE is
998 // initialized already.
999 //
1000 // <srcblock>
1001 // REAL X(130, 56)
1002 // COMPLEX Y(65, 56)
1003 // EQUIVALENCE ( X(1, 1), Y(1, 1) )
1004 // ...
1005 // CALL SCFFTM(1, 128, 50, 1.0, X, 130, Y, 65, TABLE, WORK, 0)
1006 // </srcblock>
1007 // </example>
1008 //
1009 // Input parameters:
1010 // <dl compact>
1011 // <dt><b>isign</b>
1012 // <dd> Integer.
1013 // Specifies whether to initialize the table array or to do the
1014 // forward or inverse Fourier transform, as follows:
1015 //
1016 // If isign = 0, the routine initializes the table array and
1017 // returns. In this case, the only arguments used or checked
1018 // are isign, n, and table.
1019 //
1020 // If isign = +1 or -1, the value of isign is the sign of the
1021 // exponent used in the FFT formula.
1022 // <dt><b>n</b>
1023 // <dd> Integer.
1024 // Size of the transforms (the number of elements in each
1025 // column of the input and output matrix to be transformed).
1026 // If n is not positive, <src>scfftm</src> or <src>csfftm</src> returns
1027 // without computing a transforms.
1028 // <dt><b>lot</b>
1029 // <dd> Integer.
1030 // The number of transforms to be computed (or "lot size").
1031 // This is the number of elements in each row of the input and
1032 // output matrix. If lot is not positive, <src>csfftm</src> or
1033 // <src>scfftm</src> returns without computing a transforms.
1034 // <dt><b>scale</b>
1035 // <dd> Scale factor.
1036 // <src>scfftm</src>: real.
1037 // <src>dzfftm</src>: double precision.
1038 // <src>csfftm</src>: real.
1039 // <src>zdfftm</src>: double precision.
1040 // Each element of the output array is multiplied by scale
1041 // after taking the transform, as defined in the preceding
1042 // formula.
1043 // <dt><b>x</b>
1044 // <dd> Input array of values to be transformed. Dimension (0:ldx-1,
1045 // 0:lot-1).
1046 // <src>scfftm</src>: real array.
1047 // <src>dzfftm</src>: double precision array.
1048 // <src>csfftm</src>: complex array.
1049 // <src>zdfftm</src>: double complex array.
1050 // <dt><b>ldx</b>
1051 // <dd> Integer.
1052 // The number of rows in the x array, as it was declared in the
1053 // calling program. That is, the leading dimension of x.
1054 // <src>scfftm, dzfftm</src>: ldx >= MAX(n, 1).
1055 // <src>csfftm, zdfftm</src>: ldx >= MAX(n/2 + 1, 1).
1056 // <dt><b>ldy</b>
1057 // <dd> Integer.
1058 // The number of rows in the y array, as it was declared in the
1059 // calling program (the leading dimension of y).
1060 // <src>scfftm, dzfftm</src>: ldy >= MAX(n/2 + 1, 1).
1061 // <src>csfftm, zdfftm</src>: ldy >= MAX(n, 1).
1062 // <dt><b>isys</b>
1063 // <dd> Integer array of dimension (0:isys(0)).
1064 // The first element of the array specifies how many more
1065 // elements are in the array. Use isys to specify certain
1066 // processor-specific parameters or options.
1067 //
1068 // If isys(0) = 0, the default values of such parameters are
1069 // used. In this case, you can specify the argument value as
1070 // the scalar integer constant 0.
1071 //
1072 // If isys(0) > 0, isys(0) gives the upper bound of the isys
1073 // array. Therefore, if il = isys(0), user-specified
1074 // parameters are expected in isys(1) through isys(il).
1075 // </dl>
1076 // Output parameters:
1077 // <dl compact>
1078 // <dt><b>y</b>
1079 // <dd> Output array of transformed values. Dimension (0:ldy-1,
1080 // 0:lot-1).
1081 // <src>scfftm</src>: complex array.
1082 // <src>dzfftm</src>: double complex array.
1083 // <src>csfftm</src>: real array.
1084 // <src>zdfftm</src>: double precision array.
1085 //
1086 // Each column of the output array, y, is the FFT of the
1087 // corresponding column of the input array, x, computed
1088 // according to the preceding formula. The output array may be
1089 // equivalenced to the input array. In that case, the transform
1090 // is done in place and the input array is overwritten with the
1091 // transformed values. In this case, the following conditions
1092 // on the leading dimensions must hold:
1093 //
1094 // <src>scfftm, dzfftm</src>: ldx = 2ldy.
1095 // <src>csfftm, zdfftm</src>: ldy = 2ldx.
1096 // <dt><b>table</b>
1097 // <dd> Real array; dimension (15 + n).
1098 // Table of factors and trigonometric functions.
1099 // This array must be initialized by a call to <src>scfftm</src> (or
1100 // <src>csfftm</src>) with isign = 0.
1101 //
1102 // If isign = 0, table is initialized to contain trigonometric
1103 // tables needed to compute an FFT of length n.
1104 // <dt><b>work</b>
1105 // <dd> Real array; dimension n.
1106 // Work array used for intermediate calculations. Its address
1107 // space must be different from that of the input and output
1108 // arrays.
1109 // </dl>
1110 // <group>
1111 static void scfftm(Int isign, Int n, Int lot, Float scale, Float*
1112  x, Int ldx, Complex* y, Int ldy, Float* table,
1113  Float* work, Int isys);
1114 static void dzfftm(Int isign, Int n, Int lot, Double scale, Double*
1115  x, Int ldx, DComplex* y, Int ldy, Double* table,
1116  Double* work, Int isys);
1117 static void csfftm(Int isign, Int n, Int lot, Float scale, Complex*
1118  x, Int ldx, Float* y, Int ldy, Float* table,
1119  Float* work, Int isys);
1120 static void zdfftm(Int isign, Int n, Int lot, Double scale, DComplex*
1121  x, Int ldx, Double* y, Int ldy, Double* table,
1122  Double* work, Int isys);
1123 // </group>
1124 
1125 // These routines compute the two-dimensional complex Fast Fourier
1126 // Transform (FFT) of the complex matrix x, and store the results in the
1127 // complex matrix y. <src>ccfft2d</src> does the complex-to-complex
1128 // transform and <src>zzfft</src> does the same for double
1129 // precision arrays.
1130 //
1131 // In FFT applications, it is customary to use zero-based subscripts; the
1132 // formulas are simpler that way. Suppose that the arrays are
1133 // dimensioned as follows:
1134 //
1135 // <srcblock>
1136 // COMPLEX X(0:n1-1, 0:n2-1)
1137 // COMPLEX Y(0:n1-1, 0:n2-1)
1138 // </srcblock>
1139 //
1140 // These routines compute the formula:
1141 //
1142 // <srcblock>
1143 // n2-1 n1-1
1144 // Y(k1, k2) = scale * Sum Sum [ X(j1, j2)*w1**(j1*k1)*w2**(j2*k2) ]
1145 // j2=0 j1=0
1146 //
1147 // for k1 = 0, ..., n1-1
1148 // k2 = 0, ..., n2-1
1149 //
1150 // where:
1151 // w1 = exp(isign*2*pi*i/n1)
1152 // w2 = exp(isign*2*pi*i/n2)
1153 // i = + sqrt(-1)
1154 // pi = 3.14159...,
1155 // isign = +1 or -1
1156 // </srcblock>
1157 //
1158 // Different authors use different conventions for which of the
1159 // transforms, isign = +1 or isign = -1, is the forward or inverse
1160 // transform, and what the scale factor should be in either case. You
1161 // can make this routine compute any of the various possible definitions,
1162 // however, by choosing the appropriate values for isign and scale.
1163 //
1164 // The relevant fact from FFT theory is this: If you take the FFT with
1165 // any particular values of isign and scale, the mathematical inverse
1166 // function is computed by taking the FFT with -isign and
1167 // 1/(n1*n2*scale). In particular, if you use isign = +1 and scale = 1.0
1168 // for the forward FFT, you can compute the inverse FFT by using isign =
1169 // -1 and scale = 1.0/(n1*n2).
1170 //
1171 // <h3>Algorithm</h3>
1172 // These routines use a routine very much like <src>ccfftm/zzfftm</src> to do
1173 // multiple FFTs first on all columns in an input matrix and then on all
1174 // of the rows.
1175 //
1176 // <h3>Initialization</h3>
1177 // The table array stores factors of n1 and n2 and also trigonometric
1178 // tables that are used in calculation of the FFT. This table must be
1179 // initialized by calling the routine with isign = 0. If the values of
1180 // the problem sizes, n1 and n2, do not change, the table does not have
1181 // to be reinitialized.
1182 //
1183 // <h3>Dimensions</h3>
1184 // In the preceding description, it is assumed that array subscripts were
1185 // zero-based, as is customary in FFT applications. Thus, the input and
1186 // output arrays are declared as follows:
1187 //
1188 // <srcblock>
1189 // COMPLEX X(0:ldx-1, 0:n2-1)
1190 // COMPLEX Y(0:ldy-1, 0:n2-1)
1191 // </srcblock>
1192 //
1193 // However, the calling sequence does not change if you prefer to use the
1194 // more customary Fortran style with subscripts starting at 1. The same
1195 // values of ldx and ldy would be passed to the subroutine even if the
1196 // input and output arrays were dimensioned as follows:
1197 //
1198 // <srcblock>
1199 // COMPLEX X(ldx, n2)
1200 // COMPLEX Y(ldy, n2)
1201 // </srcblock>
1202 //
1203 // <example>
1204 // All examples here are for Origin series only.
1205 //
1206 // Example 1: Initialize the TABLE array in preparation for doing a
1207 // two-dimensional FFT of size 128 by 256. In this case only the isign,
1208 // n1, n2, and table arguments are used; you can use dummy arguments or
1209 // zeros for other arguments.
1210 //
1211 // <srcblock>
1212 // REAL TABLE ((30 + 256) + (30 + 512))
1213 // CALL CCFFT2D (0, 128, 256, 0.0, DUMMY, 1, DUMMY, 1,
1214 // & TABLE, DUMMY, 0)
1215 // </srcblock>
1216 //
1217 // Example 2: X and Y are complex arrays of dimension (0:128, 0:255).
1218 // The first 128 elements of each column contain data. For performance
1219 // reasons, the extra element forces the leading dimension to be an odd
1220 // number. Take the two-dimensional FFT of X and store it in Y.
1221 // Initialize the TABLE array, as in example 1.
1222 //
1223 // <srcblock>
1224 // COMPLEX X(0:128, 0:255)
1225 // COMPLEX Y(0:128, 0:255)
1226 // REAL TABLE((30 + 256) + (30 + 512))
1227 // REAL WORK(2*128*256)
1228 // ...
1229 // CALL CCFFT2D(0, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
1230 // CALL CCFFT2D(1, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
1231 // </srcblock>
1232 //
1233 // Example 3: With X and Y as in example 2, take the inverse FFT of Y
1234 // and store it back in X. The scale factor 1/(128*256) is used. Assume
1235 // that the TABLE array is already initialized.
1236 //
1237 // <srcblock>
1238 // CALL CCFFT2D(-1, 128, 256, 1.0/(128.0*256.0), Y, 129,
1239 // & X, 129, TABLE, WORK, 0)
1240 // </srcblock>
1241 //
1242 // Example 4: Perform the same computation as in example 2, but assume
1243 // that the lower bound of each array is 1, rather than 0. The
1244 // subroutine calls are not changed.
1245 //
1246 // <srcblock>
1247 // COMPLEX X(129, 256)
1248 // COMPLEX Y(129, 256)
1249 // ...
1250 // CALL CCFFT2D(0, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
1251 // CALL CCFFT2D(1, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
1252 // </srcblock>
1253 //
1254 // Example 5: Perform the same computation as in example 4, but put the
1255 // output back in array X to save storage space. Assume that the TABLE
1256 // array is already initialized.
1257 //
1258 // <srcblock>
1259 // COMPLEX X(129, 256)
1260 // ...
1261 // CALL CCFFT2D(1, 128, 256, 1.0, X, 129, X, 129, TABLE, WORK, 0)
1262 // </srcblock>
1263 // </example>
1264 //
1265 // Input parameters:
1266 // <dl compact>
1267 // <dt><b>isign</b>
1268 // <dd> Integer.
1269 // Specifies whether to initialize the table array or to do the
1270 // forward or inverse transform as follows:
1271 //
1272 // If isign = 0, the routine initializes the table array and
1273 // returns. In this case, the only arguments used or checked
1274 // are isign, n1, n2, table.
1275 //
1276 // If isign = +1 or -1, the value of isign is the sign of the
1277 // exponent used in the FFT formula.
1278 // <dt><b>n1</b>
1279 // <dd> Integer.
1280 // Transform size in the first dimension. If n1 is not
1281 // positive, the routine returns without performing a
1282 // transform.
1283 // <dt><b>n2</b>
1284 // <dd> Integer.
1285 // Transform size in the second dimension. If n2 is not
1286 // positive, the routine returns without performing a
1287 // transform.
1288 // <dt><b>scale</b>
1289 // <dd> Scale factor.
1290 // ccfft2d: real.
1291 // zzfft2d: double precision.
1292 // Each element of the output array is multiplied by scale
1293 // factor after taking the Fourier transform, as defined
1294 // previously.
1295 // <dt><b>x</b>
1296 // <dd> Array of dimension (0:ldx-1, 0:n2-1).
1297 // ccfft2d: complex array.
1298 // zzfft2d: double complex array.
1299 // Input array of values to be transformed.
1300 // <dt><b>ldx</b>
1301 // <dd> Integer.
1302 // The number of rows in the x array, as it was declared in the
1303 // calling program (the leading dimension of x). ldx >=
1304 // MAX(n1, 1).
1305 // <dt><b>ldy</b>
1306 // <dd> Integer.
1307 //
1308 // The number of rows in the y array, as it was declared in the
1309 // calling program (the leading dimension of y). ldy >=
1310 // MAX(n1, 1).
1311 // <dt><b>isys</b>
1312 // <dd> Algorithm used; value dependent on hardware system. Currently, no
1313 // special options are supported; therefore, you must always specify
1314 // an isys argument as constant 0.
1315 // </dl>
1316 // Output parameters:
1317 // <dl compact>
1318 // <dt><b>y</b>
1319 // <dd> Array of dimension (0:ldy-1, 0:n2-1).
1320 // ccfft2d: complex array.
1321 // zzfft2d: double complex array.
1322 // Output array of transformed values. The output array may be
1323 // the same as the input array, in which case, the transform is
1324 // done in place (the input array is overwritten with the
1325 // transformed values). In this case, it is necessary that
1326 // ldx = ldy.
1327 // <dt><b>table</b>
1328 // <dd> Real array; dimension (30+ 2 * n1) + (30 + 2 * n2).
1329 //
1330 // Table of factors and trigonometric functions.
1331 //
1332 // If isign = 0, the routine initializes table (table is output
1333 // only).
1334 //
1335 // If isign = +1 or -1, the values in table are assumed to be
1336 // initialized already by a prior call with isign = 0 (table is
1337 // input only).
1338 // <dt><b>work</b>
1339 // <dd> Real array; dimension 2 * (n1*n2).
1340 //
1341 // Work array. This is a scratch array used for intermediate
1342 // calculations. Its address space must be different from that
1343 // of the input and output arrays.
1344 // </dl>
1345 // <group>
1346 static void ccfft2d(Int isign, Int n1, Int n2, Float scale, Complex*
1347  x, Int ldx, Complex* y, Int ldy, Float* table,
1348  Float* work, Int isys);
1349 static void zzfft2d(Int isign, Int n1, Int n2, Double scale, DComplex*
1350  x, Int ldx, DComplex* y, Int ldy, Double* table,
1351  Double* work, Int isys);
1352 // </group>
1353 
1354 // <src>scfft2d/dzfft2d</src> computes the two-dimensional Fast Fourier
1355 // Transform (FFT) of the real matrix x, and it stores the results in the
1356 // complex matrix y. <src>csfft2d/zdfft2d</src> computes the corresponding
1357 // inverse transform.
1358 //
1359 // In FFT applications, it is customary to use zero-based subscripts; the
1360 // formulas are simpler that way. First the function of <src>scfft2d</src> is
1361 // described. Suppose the arrays are dimensioned as follows:
1362 //
1363 // <srcblock>
1364 // REAL X(0:ldx-1, 0:n2-1)
1365 // COMPLEX Y(0:ldy-1, 0:n2-1)
1366 //
1367 // where ldx >= n1 ldy >= (n1/2) + 1.
1368 // </srcblock>
1369 //
1370 // <src>scfft2d</src> computes the formula:
1371 //
1372 // <srcblock>
1373 // n2-1 n1-1
1374 // Y(k1, k2) = scale * Sum Sum [ X(j1, j2)*w1**(j1*k1)*w2**(j2*k2) ]
1375 // j2=0 j1=0
1376 //
1377 // for k1 = 0, ..., n1/2 + 1
1378 // k2 = 0, ..., n2-1
1379 //
1380 // where:
1381 // w1 = exp(isign*2*pi*i/n1)
1382 // w2 = exp(isign*2*pi*i/n2)
1383 // i = + sqrt(-1)
1384 // pi = 3.14159...,
1385 // isign = +1 or -1
1386 // </srcblock>
1387 //
1388 // Different authors use different conventions for which of the
1389 // transforms, isign = +1 or isign = -1, is the forward or inverse
1390 // transform, and what the scale factor should be in either case. You
1391 // can make these routines compute any of the various possible
1392 // definitions, however, by choosing the appropriate values for isign and
1393 // scale.
1394 //
1395 // The relevant fact from FFT theory is this: If you take the FFT with
1396 // any particular values of isign and scale, the mathematical inverse
1397 // function is computed by taking the FFT with -isign and 1/(n1 * n2 *
1398 // scale). In particular, if you use isign = +1 and scale = 1.0 for the
1399 // forward FFT, you can compute the inverse FFT by using isign = -1 and
1400 // scale = 1.0/(n1 . n2).
1401 //
1402 // <src>scfft2d</src> is very similar in function to <src>ccfft2d</src>, but
1403 // it takes the real-to-complex transform in the first dimension, followed by
1404 // the complex-to-complex transform in the second dimension.
1405 //
1406 // <src>csfft2d</src> does the reverse. It takes the complex-to-complex FFT
1407 // in the second dimension, followed by the complex-to-real FFT in the first
1408 // dimension.
1409 //
1410 // See the <src>scfft</src> man page for more information about real-to-complex
1411 // and complex-to-real FFTs. The two-dimensional analog of the conjugate
1412 // formula is as follows:
1413 //
1414 // <srcblock>
1415 // Y = conjg Y
1416 // k , k n1 - k , n2 - k
1417 // 1 2 1 2
1418 //
1419 // for n1/2 < k <= n1 - 1
1420 // 1
1421 //
1422 // 0 <= k <= n2 - 1
1423 // 2
1424 // where the notation conjg(z) represents the complex conjugate of z.
1425 // </srcblock>
1426 //
1427 // Thus, you have to compute only (slightly more than) half of the output
1428 // values, namely:
1429 //
1430 // <srcblock>
1431 // Y for 0 <= k <= n1/2 0 <= k <= n2 - 1
1432 // k , k 1 2
1433 // 1 2
1434 // </srcblock>
1435 //
1436 // <h3>Algorithm</h3>
1437 // <src>scfft2d</src> uses a routine similar to <src>scfftm</src> to do a
1438 // real-to-complex FFT on the columns, then uses a routine similar to
1439 // <src>ccfftm</src> to do a complex-to-complex FFT on the rows.
1440 //
1441 // <src>csfft2d</src> uses a routine similar to <src>ccfftm</src> to do a
1442 // complex-to-complex FFT on the rows, then uses a routine similar to
1443 // <src>csfftm</src> to do a complex-to-real FFT on the columns.
1444 //
1445 // <h3>Table Initialization</h3>
1446 // The table array stores factors of n1 and n2, and trigonometric tables
1447 // that are used in calculation of the FFT. table must be initialized by
1448 // calling the routine with isign = 0. table does not have to be
1449 // reinitialized if the values of the problem sizes, n1 and n2, do not
1450 // change.
1451 //
1452 // <h3>Dimensions</h3>
1453 // In the preceding description, it is assumed that array subscripts were
1454 // zero-based, as is customary in FFT applications. Thus, the input and
1455 // output arrays are declared:
1456 //
1457 // <srcblock>
1458 // REAL X(0:ldx-1, 0:n2-1)
1459 // COMPLEX Y(0:ldy-1, 0:n2-1)
1460 // </srcblock>
1461 //
1462 // No change is made in the calling sequence, however, if you prefer to
1463 // use the more customary Fortran style with subscripts starting at 1.
1464 // The same values of ldx and ldy would be passed to the subroutine even
1465 // if the input and output arrays were dimensioned as follows:
1466 //
1467 // <srcblock>
1468 // REAL X(ldx, n2)
1469 // COMPLEX Y(ldy, n2)
1470 // </srcblock>
1471 //
1472 // <example>
1473 // The following examples are for Origin series only.
1474 //
1475 // Example 1: Initialize the TABLE array in preparation for doing a
1476 // two-dimensional FFT of size 128 by 256. In this case, only the isign,
1477 // n1, n2, and table arguments are used; you can use dummy arguments or
1478 // zeros for other arguments.
1479 //
1480 // <srcblock>
1481 // REAL TABLE ((15 + 128) + 2(15 + 256))
1482 // CALL SCFFT2D (0, 128, 256, 0.0, DUMMY, 1, DUMMY, 1,
1483 // & TABLE, DUMMY, 0)
1484 // </srcblock>
1485 //
1486 // Example 2: X is a real array of size (0:128, 0: 255), and Y is a
1487 // complex array of dimension (0:64, 0:255). The first 128 elements of
1488 // each column of X contain data; for performance reasons, the extra
1489 // element forces the leading dimension to be an odd number. Take the
1490 // two-dimensional FFT of X and store it in Y. Initialize the TABLE
1491 // array, as in example 1.
1492 //
1493 // <srcblock>
1494 // REAL X(0:128, 0:255)
1495 // COMPLEX Y(0:64, 0:255)
1496 // REAL TABLE ((15 + 128) + 2(15 + 256))
1497 // REAL WORK(128*256)
1498 // ...
1499 // CALL SCFFT2D(0, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
1500 // CALL SCFFT2D(1, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
1501 // </srcblock>
1502 //
1503 // Example 3: With X and Y as in example 2, take the inverse FFT of Y
1504 // and store it back in X. The scale factor 1/(128*256) is used. Assume
1505 // that the TABLE array is initialized already.
1506 //
1507 // <srcblock>
1508 // CALL CSFFT2D(-1, 128, 256, 1.0/(128.0*256.0), Y, 65,
1509 // & X, 130, TABLE, WORK, 0)
1510 // </srcblock>
1511 //
1512 // Example 4: Perform the same computation as in example 2, but assume
1513 // that the lower bound of each array is 1, rather than 0. No change is
1514 // needed in the subroutine calls.
1515 //
1516 // <srcblock>
1517 // REAL X(129, 256)
1518 // COMPLEX Y(65, 256)
1519 // ...
1520 // CALL SCFFT2D(0, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
1521 // CALL SCFFT2D(1, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
1522 // </srcblock>
1523 //
1524 // Example 5: Perform the same computation as in example 4, but
1525 // equivalence the input and output arrays to save storage space. In
1526 // this case, a row must be added to X, because it is equivalenced to a
1527 // complex array. Assume that TABLE is already initialized.
1528 //
1529 // <srcblock>
1530 // REAL X(130, 256)
1531 // COMPLEX Y(65, 256)
1532 // EQUIVALENCE ( X(1, 1), Y(1, 1) )
1533 // ...
1534 // CALL SCFFT2D(1, 128, 256, 1.0, X, 130, Y, 65, TABLE, WORK, 0)
1535 // </srcblock>
1536 // </example>
1537 //
1538 // Input parameters:
1539 // <dl compact>
1540 // <dt><b>isign</b>
1541 // <dd> Integer.
1542 // Specifies whether to initialize the table array or to do the
1543 // forward or inverse Fourier transform, as follows:
1544 //
1545 // If isign = 0, the routine initializes the table array and
1546 // returns. In this case, the only arguments used or checked
1547 // are isign, n, and table.
1548 //
1549 // If isign = +1 or -1, the value of isign is the sign of the
1550 // exponent used in the FFT formula.
1551 // <dt><b>n1</b>
1552 // <dd> Integer.
1553 // Transform size in the first dimension. If n1 is not
1554 // positive, <src>scfft2d</src> returns without calculating a
1555 // transform.
1556 // <dt><b>n2</b>
1557 // <dd> Integer.
1558 // Transform size in the second dimension. If n2 is not
1559 // positive, <src>scfft2d</src> returns without calculating a
1560 // transform.
1561 // <dt><b>scale</b>
1562 // <dd> Scale factor.
1563 // <src>scfft2d</src>: real.
1564 // <src>dzfft2d</src>: double precision.
1565 // <src>csfft2d</src>: real.
1566 // <src>zdfft2d</src>: double precision.
1567 // Each element of the output array is multiplied by scale
1568 // factor after taking the Fourier transform, as defined
1569 // previously.
1570 // <dt><b>x</b>
1571 // <dd> Array of dimension (0:ldx-1, 0:n2-1).
1572 // <src>scfft2d</src>: real array.
1573 // <src>dzfft2d</src>: double precision array.
1574 // <src>csfft2d</src>: complex array.
1575 // <src>zdfft2d</src>: double complex array.
1576 //
1577 // Array of values to be transformed.
1578 // <dt><b>ldx</b>
1579 // <dd> Integer.
1580 // The number of rows in the x array, as it was declared in the
1581 // calling program. That is, the leading dimension of x.
1582 // <src>scfft2d, dzfft2d</src>: ldx >= MAX(n1, 1).
1583 // <src>csfft2d, zdfft2d</src>: ldx >= MAX(n1/2 + 1, 1).
1584 // <dt><b>ldy</b>
1585 // <dd> Integer.
1586 //
1587 // The number of rows in the y array, as it was declared in the
1588 // calling program (the leading dimension of y).
1589 //
1590 // <src>scfft2d, dzfft2d</src>: ldy >= MAX(n1/2 + 1, 1).
1591 // <src>csfft2d, zdfft2d</src>: ldy >= MAX(n1 + 2, 1).
1592 //
1593 // In the complex-to-real routine, two extra elements are in
1594 // the first dimension (ldy >= n1 + 2, rather than just ldy >=
1595 // n1). These elements are needed for intermediate storage
1596 // during the computation. On exit, their value is undefined.
1597 // <dt><b>isys</b>
1598 // <dd> Integer array of dimension (0:isys(0)).
1599 // The first element of the array specifies how many more
1600 // elements are in the array. Use isys to specify certain
1601 // processor-specific parameters or options.
1602 //
1603 // If isys(0) = 0, the default values of such parameters are
1604 // used. In this case, you can specify the argument value as
1605 // the scalar integer constant 0.
1606 //
1607 // If isys(0) > 0, isys(0) gives the upper bound of the isys
1608 // array. Therefore, if il = isys(0), user-specified
1609 // parameters are expected in isys(1) through isys(il).
1610 // <dt><b>isys</b>
1611 // <dd> Algorithm used; value dependent on hardware system. Currently, no
1612 // special options are supported; therefore, you must always specify
1613 // an isys argument as constant 0.
1614 // </dl>
1615 // Output parameters:
1616 // <dl compact>
1617 // <dt><b>y</b>
1618 // <dd> <src>scfft2d</src>: complex array.
1619 // <src>dzfft2d</src>: double complex array.
1620 // <src>csfft2d</src>: real array.
1621 // <src>zdfft2d</src>: double precision array.
1622 //
1623 // Output array of transformed values. The output array can be
1624 // the same as the input array, in which case, the transform is
1625 // done in place and the input array is overwritten with the
1626 // transformed values. In this case, it is necessary that the
1627 // following equalities hold:
1628 //
1629 // <src>scfft2d, dzfft2d</src>: ldx = 2 * ldy.
1630 // <src>csfft2d, zdfft2d</src>: ldy = 2 * ldx.
1631 // <dt><b>table</b>
1632 // <dd> Real array; dimension (15 + n1) + 2(15 + n2).
1633 //
1634 // Table of factors and trigonometric functions.
1635 //
1636 // If isign = 0, the routine initializes table (table is output
1637 // only).
1638 //
1639 // If isign = +1 or -1, the values in table are assumed to be
1640 // initialized already by a prior call with isign = 0 (table is
1641 // input only).
1642 // <dt><b>work</b>
1643 // <dd> Real array; dimension (n1 * n2).
1644 //
1645 // Work array. This is a scratch array used for intermediate
1646 // calculations. Its address space must be different from that
1647 // of the input and output arrays.
1648 // </dl>
1649 // <group>
1650 static void scfft2d(Int isign, Int n1, Int n2, Float scale, Float*
1651  x, Int ldx, Complex* y, Int ldy, Float* table,
1652  Float* work, Int isys);
1653 static void dzfft2d(Int isign, Int n1, Int n2, Double scale, Double*
1654  x, Int ldx, DComplex* y, Int ldy, Double* table,
1655  Double* work, Int isys);
1656 static void csfft2d(Int isign, Int n1, Int n2, Float scale, Complex*
1657  x, Int ldx, Float* y, Int ldy, Float* table,
1658  Float* work, Int isys);
1659 static void zdfft2d(Int isign, Int n1, Int n2, Double scale, DComplex*
1660  x, Int ldx, Double* y, Int ldy, Double* table,
1661  Double* work, Int isys);
1662 // </group>
1663 
1664 // These routines compute the three-dimensional complex FFT of the
1665 // complex matrix X, and store the results in the complex matrix Y.
1666 //
1667 // In FFT applications, it is customary to use zero-based subscripts; the
1668 // formulas are simpler that way. So suppose the arrays are dimensioned
1669 // as follows:
1670 //
1671 // <srcblock>
1672 // COMPLEX X(0:n1-1, 0:n2-1, 0:n3-1)
1673 // COMPLEX Y(0:n1-1, 0:n2-1, 0:n3-1)
1674 // </srcblock>
1675 //
1676 // These routines compute the formula:
1677 //
1678 // <srcblock>
1679 // Y(k1,k2,k3) =
1680 // n1-1 n2-1 n3-1
1681 // scale * Sum Sum Sum [X(j1,j2,j3)*w1**(j1*k1)*w2**(j2*k2)*w3**(j3*k3)]
1682 // j1=0 j2=0 j3=0
1683 //
1684 // for k1 = 0, ..., n1 - 1,
1685 // k2 = 0, ..., n2 - 1,
1686 // k3 = 0, ..., n3 - 1,
1687 //
1688 // where:
1689 // w1 = exp(isign*2*pi*i/n1),
1690 // w2 = exp(isign*2*pi*i/n2),
1691 // w3 = exp(isign*2*pi*i/n3),
1692 // i = + sqrt(-1)
1693 // pi = 3.14159...
1694 // isign = +1 or -1
1695 // </srcblock>
1696 //
1697 // Different authors use different conventions for which of the
1698 // transforms, isign = +1 or isign = -1, is the forward or inverse
1699 // transform, and what the scale factor should be in either case. You
1700 // can make this routine compute any of the various possible definitions,
1701 // however, by choosing the appropriate values for isign and scale.
1702 //
1703 // The relevant fact from FFT theory is this: If you take the FFT with
1704 // any particular values of isign and scale, the mathematical inverse
1705 // function is computed by taking the FFT with -isign and 1/(n1 * n2 * n3
1706 // * scale). In particular, if you use isign = +1 and scale = 1.0 for
1707 // the forward FFT, you can compute the inverse FFT by using isign = -1
1708 // and scale = 1/(n1 . n2 . n3).
1709 //
1710 // <example>
1711 // The following examples are for Origin series only.
1712 //
1713 // Example 1: Initialize the TABLE array in preparation for doing a
1714 // three-dimensional FFT of size 128 by 128 by 128. In this case, only
1715 // the isign, n1, n2, n3, and table arguments are used; you can use dummy
1716 // arguments or zeros for other arguments.
1717 //
1718 // <srcblock>
1719 // REAL TABLE ((30 + 256) + (30 + 256) + (30 + 256))
1720 // CALL CCFFT3D (0, 128, 128, 128, 0.0, DUMMY, 1, 1, DUMMY, 1, 1,
1721 // & TABLE, DUMMY, 0)
1722 // </srcblock>
1723 //
1724 // Example 2: X and Y are complex arrays of dimension (0:128, 0:128,
1725 // 0:128). The first 128 elements of each dimension contain data; for
1726 // performance reasons, the extra element forces the leading dimensions
1727 // to be odd numbers. Take the three-dimensional FFT of X and store it
1728 // in Y. Initialize the TABLE array, as in example 1.
1729 //
1730 // <srcblock>
1731 // COMPLEX X(0:128, 0:128, 0:128)
1732 // COMPLEX Y(0:128, 0:128, 0:128)
1733 // REAL TABLE ((30+256) + (30 + 256) + (30 + 256))
1734 // REAL WORK 2(128*128*128)
1735 // ...
1736 // CALL CCFFT3D(0, 128, 128, 128, 1.0, DUMMY, 1, 1,
1737 // & DUMMY, 1, 1, TABLE, WORK, 0)
1738 // CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
1739 // & Y, 129, 129, TABLE, WORK, 0)
1740 // </srcblock>
1741 //
1742 // Example 3: With X and Y as in example 2, take the inverse FFT of Y
1743 // and store it back in X. The scale factor 1.0/(128.0**3) is used.
1744 // Assume that the TABLE array is already initialized.
1745 //
1746 // <srcblock>
1747 // CALL CCFFT3D(-1, 128, 128, 128, 1.0/(128.0**3), Y, 129, 129,
1748 // & X, 129, 129, TABLE, WORK, 0)
1749 // </srcblock>
1750 //
1751 // Example 4: Perform the same computation as in example 2, but assume
1752 // that the lower bound of each array is 1, rather than 0. The
1753 // subroutine calls do not change.
1754 //
1755 // <srcblock>
1756 // COMPLEX X(129, 129, 129)
1757 // COMPLEX Y(129, 129, 129)
1758 // ...
1759 // CALL CCFFT3D(0, 128, 128, 128, 1.0, DUMMY, 1, 1,
1760 // & DUMMY, 1, 1, TABLE, WORK, 0)
1761 // CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
1762 // & Y, 129, 129, TABLE, WORK, 0)
1763 // </srcblock>
1764 //
1765 // Example 5: Perform the same computation as in example 4, but put the
1766 // output back in the array X to save storage space. Assume that the
1767 // TABLE array is already initialized.
1768 //
1769 // <srcblock>
1770 // COMPLEX X(129, 129, 129)
1771 // ...
1772 // CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
1773 // & X, 129, 129, TABLE, WORK, 0)
1774 // </srcblock>
1775 // </example>
1776 //
1777 // Input parameters:
1778 // <dl compact>
1779 // <dt><b>isign</b>
1780 // <dd> Integer.
1781 // Specifies whether to initialize the table array or to do the
1782 // forward or inverse Fourier transform, as follows:
1783 //
1784 // If isign = 0, the routine initializes the table array and
1785 // returns. In this case, the only arguments used or checked
1786 // are isign, n1, n2, n3, and table.
1787 //
1788 // If isign = +1 or -1, the value of isign is the sign of the
1789 // exponent used in the FFT formula.
1790 //
1791 // <dt><b>n1</b>
1792 // <dd> Integer.
1793 // Transform size in the first dimension. If n1 is not
1794 // positive, the routine returns without computing a transform.
1795 //
1796 // <dt><b>n2</b>
1797 // <dd> Integer.
1798 // Transform size in the second dimension. If n2 is not
1799 // positive, the routine returns without computing a transform.
1800 //
1801 // <dt><b>n3</b>
1802 // <dd> Integer.
1803 // Transform size in the third dimension. If n3 is not
1804 // positive, the routine returns without computing a transform.
1805 //
1806 // <dt><b>scale</b>
1807 // <dd> Scale factor.
1808 // <src>ccfft3d</src>: real.
1809 // <src>zzfft3d</src>: double precision.
1810 //
1811 // Each element of the output array is multiplied by scale
1812 // after taking the Fourier transform, as defined previously.
1813 //
1814 // <dt><b>x</b>
1815 // <dd> Array of dimension (0:ldx-1, 0:ldx2-1, 0:n3-1).
1816 // <src>ccfft3d</src>: complex array.
1817 // <src>zzfft3d</src>: double complex array.
1818 //
1819 // Input array of values to be transformed.
1820 //
1821 // <dt><b>ldx</b>
1822 // <dd> Integer.
1823 // The first dimension of x, as it was declared in the calling
1824 // program (the leading dimension of x). ldx >= MAX(n1, 1).
1825 //
1826 // <dt><b>ldx2</b>
1827 // <dd> Integer.
1828 // The second dimension of x, as it was declared in the calling
1829 // program. ldx2 >= MAX(n2, 1).
1830 //
1831 // <dt><b>ldy</b>
1832 // <dd> Integer.
1833 // The first dimension of y, as it was declared in the calling
1834 // program (the leading dimension of y). ldy >= MAX(n1, 1).
1835 //
1836 // <dt><b>ldy2</b>
1837 // <dd> Integer.
1838 // The second dimension of y, as it was declared in the calling
1839 // program. ldy2 >= MAX(n2, 1).
1840 //
1841 // <dt><b>isys</b>
1842 // <dd> Algorithm used; value dependent on hardware system. Currently, no
1843 // special options are supported; therefore, you must always specify
1844 // an isys argument as constant 0.
1845 //
1846 // isys = 0 or 1 depending on the amount of workspace the user
1847 // can provide to the routine.
1848 // </dl>
1849 // Output parameters:
1850 // <dl compact>
1851 // <dt><b>y</b>
1852 // <dd> Array of dimension (0:ldy-1, 0:ldy2-1, 0:n3-1).
1853 // <src>ccfft3d</src>: complex array.
1854 // <src>zzfft3d</src>: double complex array.
1855 //
1856 // Output array of transformed values. The output array may be
1857 // the same as the input array, in which case, the transform is
1858 // done in place; that is, the input array is overwritten with
1859 // the transformed values. In this case, it is necessary that
1860 // ldx = ldy, and ldx2 = ldy2.
1861 //
1862 // <dt><b>table</b>
1863 // <dd> Real array; dimension 30 + 2 * n1) + (30 + 2 * n2) + (30 + 2 * n3).
1864 //
1865 // Table of factors and trigonometric functions. If isign = 0,
1866 // the routine initializes table (table is output only). If
1867 // isign = +1 or -1, the values in table are assumed to be
1868 // initialized already by a prior call with isign = 0 (table is
1869 // input only).
1870 //
1871 // <dt><b>work</b>
1872 // <dd> Real array; dimension (n1 * n2 * n3).
1873 //
1874 // Work array. This is a scratch array used for intermediate
1875 // calculations. Its address space must be different from that
1876 // of the input and output arrays.
1877 //
1878 // </dl>
1879 // <group>
1880 static void ccfft3d(Int isign, Int n1, Int n2, Int n3, Float scale,
1881  Complex* x, Int ldx, Int ldx2, Complex* y, Int ldy,
1882  Int ldy2, Float* table, Float* work, Int isys);
1883 static void zzfft3d(Int isign, Int n1, Int n2, Int n3, Double scale,
1884  DComplex* x, Int ldx, Int ldx2, DComplex* y, Int
1885  ldy, Int ldy2, Double* table, Double* work, Int
1886  isys);
1887 // </group>
1888 
1889 // These are C++ wrapper functions for the 3D real-to-complex and
1890 // complex-to-real transform routines in the SGI/Cray Scientific Library
1891 // (SCSL). The purpose of these definitions is to overload the functions so
1892 // that C++ users can access the functions in SCSL with identical function
1893 // names.
1894 //
1895 // <note role=warning>
1896 // Currently, the SCSL is available only on SGI machines.
1897 // </note>
1898 //
1899 // <src>scfft3d/dzfft3d</src> computes the three-dimensional Fast Fourier
1900 // Transform (FFT) of the real matrix X, and it stores the results in the
1901 // complex matrix Y. <src>csfft3d/zdfft3d</src> computes the corresponding
1902 // inverse transform.
1903 //
1904 // In FFT applications, it is customary to use zero-based subscripts; the
1905 // formulas are simpler that way. First, the function of <src>SCFFT3D</src> is
1906 // described. Suppose the arrays are dimensioned as follows:
1907 //
1908 // <srcblock>
1909 // REAL X(0:ldx-1, 0:ldx2-1, 0:n3-1)
1910 // COMPLEX Y(0:ldy-1, 0:ldy2-1, 0:n3-1)
1911 // </srcblock>
1912 //
1913 // <src>scfft3d</src> computes the formula:
1914 //
1915 // <srcblock>
1916 // Y(k1,k2,k3) =
1917 // n1-1 n2-1 n3-1
1918 // scale * Sum Sum Sum [X(j1,j2,j3)*w1**(j1*k1)*w2**(j2*k2)*w3**(j3*k3)]
1919 // j1=0 j2=0 j3=0
1920 //
1921 // for k1 = 0, ..., n1/2,
1922 // k2 = 0, ..., n2 - 1,
1923 // k3 = 0, ..., n3 - 1,
1924 //
1925 // where:
1926 // w1 = exp(isign*2*pi*i/n1),
1927 // w2 = exp(isign*2*pi*i/n2),
1928 // w3 = exp(isign*2*pi*i/n3),
1929 // i = + sqrt(-1)
1930 // pi = 3.14159...
1931 // isign = +1 or -1
1932 // </srcblock>
1933 //
1934 // Different authors use different conventions for which of the
1935 // transforms, isign = +1 or isign = -1, is the forward or inverse
1936 // transform, and what the scale factor should be in either case. You
1937 // can make these routines compute any of the various possible
1938 // definitions, however, by choosing the appropriate values for isign and
1939 // scale.
1940 //
1941 // The relevant fact from FFT theory is this: If you take the FFT with
1942 // any particular values of isign and scale, the mathematical inverse
1943 // function is computed by taking the FFT with -isign and 1/(n1 * n2 * n3
1944 // * scale). In particular, if you use isign = +1 and scale = 1.0 for
1945 // the forward FFT, you can compute the inverse FFT by isign = -1 and
1946 //
1947 // <srcblock>
1948 // scale = 1.0/(n1*n2*n3).
1949 // </srcblock>
1950 //
1951 // <src>scfft3d</src> is very similar in function to <src>ccfft3d</src>, but
1952 // it takes the real-to-complex transform in the first dimension, followed by
1953 // the complex-to-complex transform in the second and third dimensions.
1954 //
1955 // <src>csfft3d</src> does the reverse. It takes the complex-to-complex FFT
1956 // in the third and second dimensions, followed by the complex-to-real FFT in
1957 // the first dimension.
1958 //
1959 // See the <src>scfftm</src> man page for more information about
1960 // real-to-complex and complex-to-real FFTs. The three dimensional analog of
1961 // the conjugate formula is as follows:
1962 //
1963 // <srcblock>
1964 // Y = conjg Y
1965 // k ,k ,k n1 - k , n2 - k , n3 - k
1966 // 1 2 3 1 2 3
1967 //
1968 // for n1/2 < k <= n1 - 1
1969 // 1
1970 //
1971 // 0 <= k <= n2 - 1
1972 // 2
1973 //
1974 // 0 <= k <= n3 - 1
1975 // 3
1976 // where the notation conjg(z) represents the complex conjugate of z.
1977 // </srcblock>
1978 //
1979 // Thus, you have to compute only (slightly more than) half out the
1980 // output values, namely:
1981 //
1982 // <srcblock>
1983 // Y
1984 // k ,k ,k
1985 // 1 2 3
1986 //
1987 // for 0 <= k <= n1/2
1988 // 1
1989 //
1990 // 0 <= k <= n2 - 1
1991 // 2
1992 //
1993 // 0 <= k <= n3 - 1
1994 // </srcblock>
1995 //
1996 // <h3>Algorithm</h3>
1997 // <src>scfft3d</src> uses a routine similar to <src>scfftm</src> to do
1998 // multiple FFTs first on all columns of the input matrix, then uses a routine
1999 // similar to <src>ccfftm</src> on all rows of the result, and then on all
2000 // planes of that result. See <src>scfftm</src> and <src>ccfftm</src> for
2001 // more information about the algorithms used.
2002 //
2003 // </example>
2004 // The following examples are for Origin series only.
2005 //
2006 // Example 1: Initialize the TABLE array in preparation for doing a
2007 // three-dimensional FFT of size 128 by 128 by 128. In this case only
2008 // the isign, n1, n2, n3, and table arguments are used; you can use dummy
2009 // arguments or zeros for other arguments.
2010 //
2011 // <srcblock>
2012 // REAL TABLE ((15 + 128) + 2(15+128) + 2( 15 + 128))
2013 // CALL SCFFT3D (0, 128, 128, 128, 0.0, DUMMY, 1, 1, DUMMY, 1, 1,
2014 // & TABLE, DUMMY, 0)
2015 // </srcblock>
2016 //
2017 // Example 2: X is a real array of size (0:128, 0:128, 0:128). The
2018 // first 128 elements of each dimension contain data; for performance
2019 // reasons, the extra element forces the leading dimensions to be odd
2020 // numbers. Y is a complex array of dimension (0:64, 0:128, 0:128).
2021 // Take the three-dimensional FFT of X and store it in Y. Initialize the
2022 // TABLE array, as in example 1.
2023 //
2024 // <srcblock>
2025 // REAL X(0:128, 0:128, 0:128)
2026 // COMPLEX Y(0:64, 0:128, 0:128)
2027 // REAL TABLE ((15+128) + 2(15 + 128) + 2(15 + 128))
2028 // REAL WORK(128*128*128)
2029 // ...
2030 // CALL SCFFT3D(0, 128, 128, 128, 1.0, X, 129, 129,
2031 // & Y, 65, 129, TABLE, WORK, 0)
2032 // CALL SCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
2033 // & Y, 65, 129, TABLE, WORK, 0)
2034 // </srcblock>
2035 //
2036 // Example 3: With X and Y as in example 2, take the inverse FFT of Y
2037 // and store it back in X. The scale factor 1/(128**3) is used. Assume
2038 // that the TABLE array is initialized already.
2039 //
2040 // <srcblock>
2041 // CALL CSFFT3D(-1, 128, 128, 128, 1.0/128.0**3, Y, 65, 129,
2042 // & X, 130, 129, TABLE, WORK, 0)
2043 // </srcblock>
2044 //
2045 // Example 4: Perform the same computation as in example 2, but assume
2046 // that the lower bound of each array is 1, rather than 0. No change is
2047 // made in the subroutine calls.
2048 //
2049 // <srcblock>
2050 // REAL X(129, 129, 129)
2051 // COMPLEX Y(65, 129, 129)
2052 // REAL TABLE ((15+128) + 2(15 + 128) + 2(15 + 128))
2053 // REAL WORK(128*128*128)
2054 // ...
2055 // CALL SCFFT3D(0, 128, 128, 128, 1.0, X, 129, 129,
2056 // & Y, 65, 129, TABLE, WORK, 0)
2057 // CALL SCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
2058 // & X, 129, 129, TABLE, WORK, 0)
2059 // </srcblock>
2060 //
2061 // Example 5: Perform the same computation as in example 4, but
2062 // equivalence the input and output arrays to save storage space. Assume
2063 // that the TABLE array is initialized already.
2064 //
2065 // <srcblock>
2066 // REAL X(130, 129, 129)
2067 // COMPLEX Y(65, 129, 129)
2068 // EQUIVALENCE (X(1, 1, 1), Y(1, 1, 1))
2069 // ...
2070 // CALL SCFFT3D(1, 128, 128, 128, 1.0, X, 130, 129,
2071 // & Y, 65, 129, TABLE, WORK, 0)
2072 // </srcblock>
2073 // </example>
2074 //
2075 // Input parameters:
2076 // <dl compact>
2077 // <dt><b>isign</b>
2078 // <dd> Integer.
2079 // Specifies whether to initialize the table array or to do the
2080 // forward or inverse Fourier transform, as follows:
2081 //
2082 // If isign = 0, the routine initializes the table array and
2083 // returns. In this case, the only arguments used or checked
2084 // are isign, n1, n2, n3, and table.
2085 //
2086 // If isign = +1 or -1, the value of isign is the sign of the
2087 // exponent used in the FFT formula.
2088 //
2089 // <dt><b>n1</b>
2090 // <dd> Integer.
2091 // Transform size in the first dimension. If n1 is not
2092 // positive, <src>scfft3d</src> returns without computing a transform.
2093 //
2094 // <dt><b>n2</b>
2095 // <dd> Integer.
2096 // Transform size in the second dimension. If n2 is not
2097 // positive, <src>scfft3d</src> returns without computing a transform.
2098 //
2099 // <dt><b>n3</b>
2100 // <dd> Integer.
2101 // Transform size in the third dimension. If n3 is not
2102 // positive, <src>scfft3d</src> returns without computing a transform.
2103 //
2104 // <dt><b>scale</b>
2105 // <dd> Scale factor.
2106 // <src>scfft3d</src>: real.
2107 // <src>dzfft3d</src>: double precision.
2108 // <src>csfft3d</src>: real.
2109 // <src>zdfft3d</src>: double precision.
2110 // Each element of the output array is multiplied by scale
2111 // after taking the Fourier transform, as defined previously.
2112 //
2113 // <dt><b>x</b>
2114 // <dd> Array of dimension (0:ldx-1, 0:ldx2-1, 0:n3-1).
2115 // <src>scfft3d</src>: real array.
2116 // <src>dzfft3d</src>: double precision array.
2117 // <src>csfft3d</src>: complex array.
2118 // <src>zdfft3d</src>: double complex array.
2119 //
2120 // Array of values to be transformed.
2121 //
2122 // <dt><b>ldx</b>
2123 // <dd> Integer.
2124 // The first dimension of x, as it was declared in the calling
2125 // program (the leading dimension of x).
2126 //
2127 // <src>scfft3d, dzfft3d</src>: ldx >= MAX(n1, 1).
2128 // <src>csfft3d, zdfft3d</src>: ldx >= MAX(n1/2 + 1, 1).
2129 //
2130 // <dt><b>ldx2</b>
2131 // <dd> Integer.
2132 // The second dimension of x, as it was declared in the calling
2133 // program. ldx2 >= MAX(n2, 1).
2134 //
2135 // <dt><b>ldy</b>
2136 // <dd> Integer.
2137 // The first dimension of y, as it was declared in the calling
2138 // program; that is, the leading dimension of y.
2139 //
2140 // <src>scfft3d, dzfft3d</src>: ldy >= MAX(n1/2 + 1, 1).
2141 // <src>csfft3d, zdfft3d</src>: ldy >= MAX(n1 + 2, 1).
2142 //
2143 // In the complex-to-real routine, two extra elements are in
2144 // the first dimension (that is, ldy >= n1 + 2, rather than
2145 // just ldy >= n1). These elements are needed for intermediate
2146 // storage during the computation. On exit, their value is
2147 // undefined.
2148 //
2149 // <dt><b>ldy2</b>
2150 // <dd> Integer.
2151 // The second dimension of y, as it was declared in the calling
2152 // program. ldy2 >= MAX(n2, 1).
2153 //
2154 // <dt><b>isys</b>
2155 // <dd> Algorithm used; value dependent on hardware system. Currently, no
2156 // special options are supported; therefore, you must always specify
2157 // an isys argument as constant 0.
2158 //
2159 // isys = 0 or 1 depending on the amount of workspace the user
2160 // can provide to the routine.
2161 // </dl>
2162 // Output parameters:
2163 // <dl compact>
2164 // <dt><b>y</b>
2165 // <dd> Array of dimension (0:ldy-1, 0:ldy2-1, 0:n3-1).
2166 // <src>scfft3d</src>: complex array.
2167 // <src>dzfft3d</src>: double complex array.
2168 // <src>csfft3d</src>: real array.
2169 // <src>zdfft3d</src>: double precision array.
2170 //
2171 // Output array of transformed values. The output array can be
2172 // the same as the input array, in which case, the transform is
2173 // done in place; that is, the input array is overwritten with
2174 // the transformed values. In this case, it is necessary that
2175 // the following equalities hold:
2176 //
2177 // <src>scfft3d, dzfft3d</src>: ldx = 2 * ldy, and ldx2 = ldy2.
2178 // <src>csfft3d, zdfft3d</src>: ldy = 2 * ldx, and ldx2 = ldy2.
2179 //
2180 // <dt><b>table</b>
2181 // <dd> Real array; dimension (15 + n1) + 2(15 + n2) + 2(15 + n3).
2182 //
2183 // Table of factors and trigonometric functions.
2184 //
2185 // This array must be initialized by a call to <src>scfft3d</src> or
2186 // <src>csfft3d</src> with isign = 0.
2187 //
2188 // If isign = 0, table is initialized to contain trigonometric
2189 // tables needed to compute a three-dimensional FFT of size n1
2190 // by n2 by n3. If isign = +1 or -1, the values in table are
2191 // assumed to be initialized already by a prior call with isign
2192 // = 0.
2193 //
2194 // <dt><b>work</b>
2195 // <dd> Real array; dimension n1 * n2 * n3.
2196 //
2197 // Work array. This is a scratch array used for intermediate
2198 // calculations. Its address space must be different from that
2199 // of the input and output arrays.
2200 //
2201 // </dl>
2202 // <group>
2203 static void scfft3d(Int isign, Int n1, Int n2, Int n3, Float scale,
2204  Float* x, Int ldx, Int ldx2, Complex* y, Int ldy,
2205  Int ldy2, Float* table, Float* work, Int isys);
2206 static void dzfft3d(Int isign, Int n1, Int n2, Int n3, Double scale,
2207  Double* x, Int ldx, Int ldx2, DComplex* y, Int
2208  ldy, Int ldy2, Double* table, Double* work, Int
2209  isys);
2210 static void csfft3d(Int isign, Int n1, Int n2, Int n3, Float scale,
2211  Complex* x, Int ldx, Int ldx2, Float* y, Int ldy,
2212  Int ldy2, Float* table, Float* work, Int isys);
2213 static void zdfft3d(Int isign, Int n1, Int n2, Int n3, Double scale,
2214  DComplex* x, Int ldx, Int ldx2, Double* y, Int
2215  ldy, Int ldy2, Double* table, Double* work, Int
2216  isys);
2217 // </group>
2218 };
2219 
2220 } //# NAMESPACE CASACORE - END
2221 
2222 #endif
static void dzfft(Int isign, Int n, Double scale, Double *x, DComplex *y, Double *table, Double *work, Int isys)
static void zzfftm(Int isign, Int n, Int lot, Double scale, DComplex *x, Int ldx, DComplex *y, Int ldy, Double *table, Double *work, Int isys)
static void dzfft3d(Int isign, Int n1, Int n2, Int n3, Double scale, Double *x, Int ldx, Int ldx2, DComplex *y, Int ldy, Int ldy2, Double *table, Double *work, Int isys)
static void scfft(Int isign, Int n, Double scale, Double *x, DComplex *y, Double *table, Double *work, Int isys)
static void zdfft(Int isign, Int n, Double scale, DComplex *x, Double *y, Double *table, Double *work, Int isys)
static void ccfft2d(Int isign, Int n1, Int n2, Float scale, Complex *x, Int ldx, Complex *y, Int ldy, Float *table, Float *work, Int isys)
These routines compute the two-dimensional complex Fast Fourier Transform (FFT) of the complex matrix...
static void zzfft3d(Int isign, Int n1, Int n2, Int n3, Double scale, DComplex *x, Int ldx, Int ldx2, DComplex *y, Int ldy, Int ldy2, Double *table, Double *work, Int isys)
static void zdfft3d(Int isign, Int n1, Int n2, Int n3, Double scale, DComplex *x, Int ldx, Int ldx2, Double *y, Int ldy, Int ldy2, Double *table, Double *work, Int isys)
static void scfft3d(Int isign, Int n1, Int n2, Int n3, Float scale, Float *x, Int ldx, Int ldx2, Complex *y, Int ldy, Int ldy2, Float *table, Float *work, Int isys)
These are C++ wrapper functions for the 3D real-to-complex and complex-to-real transform routines in ...
static void dzfft2d(Int isign, Int n1, Int n2, Double scale, Double *x, Int ldx, DComplex *y, Int ldy, Double *table, Double *work, Int isys)
static void csfftm(Int isign, Int n, Int lot, Float scale, Complex *x, Int ldx, Float *y, Int ldy, Float *table, Float *work, Int isys)
static void scfftm(Int isign, Int n, Int lot, Float scale, Float *x, Int ldx, Complex *y, Int ldy, Float *table, Float *work, Int isys)
scfftm/dzfftm computes the FFT of each column of the real matrix X, and it stores the results in the ...
static void ccfft(Int isign, Int n, Float scale, Complex *x, Complex *y, Float *table, Float *work, Int isys)
These routines compute the Fast Fourier Transform (FFT) of the complex vector x, and store the result...
static void ccfft3d(Int isign, Int n1, Int n2, Int n3, Float scale, Complex *x, Int ldx, Int ldx2, Complex *y, Int ldy, Int ldy2, Float *table, Float *work, Int isys)
These routines compute the three-dimensional complex FFT of the complex matrix X, and store the resul...
static void zdfftm(Int isign, Int n, Int lot, Double scale, DComplex *x, Int ldx, Double *y, Int ldy, Double *table, Double *work, Int isys)
static void zzfft(Int isign, Int n, Double scale, DComplex *x, DComplex *y, Double *table, Double *work, Int isys)
static void scfft2d(Int isign, Int n1, Int n2, Float scale, Float *x, Int ldx, Complex *y, Int ldy, Float *table, Float *work, Int isys)
scfft2d/dzfft2d computes the two-dimensional Fast Fourier Transform (FFT) of the real matrix x,...
static void scfft(Int isign, Int n, Float scale, Float *x, Complex *y, Float *table, Float *work, Int isys)
scfft/dzfft computes the FFT of the real array x, and it stores the results in the complex array y.
static void csfft2d(Int isign, Int n1, Int n2, Float scale, Complex *x, Int ldx, Float *y, Int ldy, Float *table, Float *work, Int isys)
static void zzfft2d(Int isign, Int n1, Int n2, Double scale, DComplex *x, Int ldx, DComplex *y, Int ldy, Double *table, Double *work, Int isys)
static void ccfftm(Int isign, Int n, Int lot, Float scale, Complex *x, Int ldx, Complex *y, Int ldy, Float *table, Float *work, Int isys)
ccfftm/zzfftm computes the FFT of each column of the complex matrix x, and stores the results in the ...
static void csfft3d(Int isign, Int n1, Int n2, Int n3, Float scale, Complex *x, Int ldx, Int ldx2, Float *y, Int ldy, Int ldy2, Float *table, Float *work, Int isys)
static void dzfftm(Int isign, Int n, Int lot, Double scale, Double *x, Int ldx, DComplex *y, Int ldy, Double *table, Double *work, Int isys)
static void ccfft(Int isign, Int n, Double scale, DComplex *x, DComplex *y, Double *table, Double *work, Int isys)
static void zdfft2d(Int isign, Int n1, Int n2, Double scale, DComplex *x, Int ldx, Double *y, Int ldy, Double *table, Double *work, Int isys)
static void csfft(Int isign, Int n, Double scale, DComplex *x, Double *y, Double *table, Double *work, Int isys)
static void csfft(Int isign, Int n, Float scale, Complex *x, Float *y, Float *table, Float *work, Int isys)
this file contains all the compiler specific defines
Definition: mainpage.dox:28
float Float
Definition: aipstype.h:54
int Int
Definition: aipstype.h:50
double Double
Definition: aipstype.h:55