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
casacore::SCSL::csfft2d
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)
casacore::SCSL::csfft3d
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)
casacore::SCSL::zzfft3d
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)
casacore::SCSL::scfft
static void scfft(Int isign, Int n, Double scale, Double *x, DComplex *y, Double *table, Double *work, Int isys)
Complexfwd_global_functions_Complexfwd::casacore::DComplex
std::complex< Double > DComplex
Definition: Complexfwd.h:50
casacore::SCSL::dzfftm
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)
casacore::SCSL::scfftm
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 ...
casacore::SCSL::zzfft2d
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)
casacore::SCSL::ccfft
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...
casacore::SCSL::ccfft3d
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...
casacore::Float
float Float
Definition: aipstype.h:54
casacore::Double
double Double
Definition: aipstype.h:55
casacore::SCSL::zdfftm
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)
casacore::SCSL::dzfft2d
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)
casacore::SCSL::scfft
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.
casacore::SCSL::csfft
static void csfft(Int isign, Int n, Double scale, DComplex *x, Double *y, Double *table, Double *work, Int isys)
casacore::SCSL::ccfft
static void ccfft(Int isign, Int n, Double scale, DComplex *x, DComplex *y, Double *table, Double *work, Int isys)
casacore::SCSL::ccfftm
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 ...
casacore::SCSL
Definition: SCSL.h:50
casacore::Int
int Int
Definition: aipstype.h:50
casacore
this file contains all the compiler specific defines
Definition: mainpage.dox:28
casacore::SCSL::dzfft
static void dzfft(Int isign, Int n, Double scale, Double *x, DComplex *y, Double *table, Double *work, Int isys)
casacore::SCSL::dzfft3d
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)
casacore::SCSL::ccfft2d
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...
casacore::SCSL::zdfft
static void zdfft(Int isign, Int n, Double scale, DComplex *x, Double *y, Double *table, Double *work, Int isys)
casacore::SCSL::scfft2d
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,...
casacore::SCSL::scfft3d
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 ...
Complexfwd_global_functions_Complexfwd::casacore::Complex
std::complex< Float > Complex
Definition: Complexfwd.h:49
casacore::SCSL::zzfft
static void zzfft(Int isign, Int n, Double scale, DComplex *x, DComplex *y, Double *table, Double *work, Int isys)
casacore::SCSL::zdfft2d
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)
casacore::SCSL::zzfftm
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)
casacore::SCSL::zdfft3d
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)
casacore::SCSL::csfftm
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)
casacore::SCSL::csfft
static void csfft(Int isign, Int n, Float scale, Complex *x, Float *y, Float *table, Float *work, Int isys)