home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / DOC / libfft_2.3.ASCII
Encoding:
Text File  |  1994-08-02  |  45.8 KB  |  1,347 lines
  1. #############################################################################
  2. #                                        #
  3. #    LIBFFT        Version 2.3  - Jan 12, 1992 -                 #
  4. #                                        #
  5. #    libfft    by  Jean-Pierre Panziera                    #
  6. #            Silicon Graphics Inc.                    #
  7. #            jpp@sgi.com                            #
  8. #                                        #
  9. #    libfft includes one, two and three dimensional FFTs            #
  10. #            complex <---> complex                    #
  11. #            real    <---> complex                    #
  12. #            single and double precision,                #
  13. #            2D and 3D are parallel                    #
  14. #                                        #
  15. #    libfft uses a "row-column" algorithm, and it is based on 1D FFTs    #
  16. #            from FFTPACK                        #
  17. #                                        #
  18. #    FFTPACK is a 1D FFTs package written by Paul N Swarztrauber        #
  19. #        at   National Center for Atmospheric Research  Boulder,CO        #
  20. #                                        #
  21. #############################################################################
  22. #                                        #
  23. #      Version 2.3 fixes a bug with SFFT2DU and DFFT2DU.            #
  24. #      Also the default value for FFT_MAXMEM is now set to 10000000         #
  25. #                                        #
  26. #############################################################################
  27. #                                        #
  28. #    The version 2.2 is equivalent to version 2.0 in terms of        #
  29. #    performances. However it adds a new family of function/routines        #
  30. #    for ( Real <---> Complex ) 1, 2 and 3 D FFTs.                #
  31. #    These new modules are identified by the letter U (e.g. SFFT2DU)        #
  32. #    They provide the same functionnality as the modules of version 2.0  #
  33. #    However, the storage scheme of the Fourier Transforms is much        #
  34. #    more simpler, for the same performance ...                #
  35. #                                        #
  36. #############################################################################
  37.  
  38.     All routines come with both Fortran and  C  interfaces.
  39.     In C two types "complex" and "zomplex" have been defined as structures of
  40.     two floating point variables ( re, im) ... see "fft.h". They are equivalent
  41.     to the "complex" and "double complex" Fortran types.
  42.  
  43.  
  44. #############################################################################
  45.  
  46. WARNING:
  47.  
  48.     The library FFTs are NOT normalized. The following successive double call :
  49.     call _fft_d( +1, ...)
  50.     call _fft_d( -1, ...)
  51.     will scale the input by N for 1D FFT, ( N1 x N2 for 2D FFT, and 
  52.      N1 x N2 x N3 for 3D FFT).
  53.  
  54. ############################################################################
  55.  
  56. The functions/routines names are:
  57.  
  58. _________________________________________________________________________
  59. |        |           |              |             |
  60. |            |   1 Dimension       |   2 Dimensions   |     3 Dimensions     |
  61. |_______________|__________________|__________________|__________________|
  62. |        |           |              |             |
  63. |  REAL        | sfft1di SFFT1DUI | sfft2di SFFT2DUI | sfft3di SFFT3DUI |
  64. |     <--->    | sfft1d  SFFT1DU  | sfft2d  SFFT2DU  | sfft3d  SFFT3DU     |
  65. |       COMPLEX | csfft1d       |              |             |
  66. |        | sprod1d SPROD1DU | sprod2d SPROD2DU | sprod3d SPROD3DU |
  67. |  Single    | sscal1d       | sscal2d          | sscal3d         |
  68. |     Precision | csscal1d       |              |             |
  69. |_______________|__________________|__________________|__________________|
  70. |        |           |              |             |
  71. |  REAL        | dfft1di DFFT1DUI | dfft2di DFFT2DUI | dfft3di DFFT3DUI |
  72. |     <--->    | dfft1d  DFFT1DU  | dfft2d  DFFT2DU  | dfft3d  DFFT3DU     |
  73. |       COMPLEX | zdfft1d       |              |             |
  74. |        | dprod1d DPROD1DU | dprod2d DPROD2DU | dprod3d DPROD3DU |
  75. |  Double    | dscal1d       | dscal2d          | dscal3d         |
  76. |     Precision | zdscal1d       |              |             |
  77. |_______________|__________________|__________________|__________________|
  78. |        |           |              |             |
  79. |  COMPLEX    | cfft1di       | cfft2di          |    cfft3di          |
  80. |     <--->    | cfft1d       | cfft2d          | cfft3d         |
  81. |    COMPLEX | cprod1d       | cprod2d          | cprod3d         |
  82. |        | cscal1d       | cscal2d          | cscal3d         |
  83. |    Single    |           |              |             |
  84. |_______________|__________________|__________________|__________________|
  85. |        |           |              |             |
  86. |  COMPLEX    | zfft1di       | zfft2di          |    zfft3di          |
  87. |     <--->    | zfft1d       | zfft2d          | zfft3d         |
  88. |    COMPLEX | zprod1d       | zprod2d          | zprod3d         |
  89. |        | zscal1d       | zscal2d          | zscal3d         |
  90. |    Double    |           |              |             |
  91. |_______________|__________________|__________________|__________________|
  92.  
  93. Note:    The new modules are written "upper-case" in the list above, It is 
  94.     advisable to use them instead of the "packed" versions (e.g. "SFFT2DU"
  95.     should be used instead of "sfft2d").
  96.  
  97. |------------------|
  98. |   Overview       |
  99. |------------------|
  100.  
  101.     The first Letter of the Function/Subroutine depends on the type:
  102.  
  103.     c____   : Single Precision Complex 
  104.     z____   : Double Precision Complex
  105.     s____   : Single Precision REAL
  106.     d____   : Double Precision REAL
  107.  
  108.     "_fft_di", "fft_dui"  functions/routines :
  109.  
  110.         Initialize the workspace array  
  111.  
  112.  
  113.  
  114.     "_fft_d", "_fft_du"  functions/routines :
  115.  
  116.         Fast Fourier Transformation.
  117.  
  118.         The library routines do the transforms IN PLACE. 
  119.         They work optimally if each of the dimension number 
  120.         n1, n2 or n3 is a multiple of 2, 3, 4 and 5.
  121.  
  122.         But they will always give you the right answer independent of
  123.         the value of N.
  124.  
  125.  
  126.  
  127.     "_prod_d", "_prod_du"  functions/routines :
  128.  
  129.         Performs the product of Fourier Transform of SAME sizes.
  130.  
  131.         The Array is multiplied term by term by those of a Filter.
  132.         This operation in the Fourier Domain is equivalent to a
  133.         Convolution in the time/space domain.
  134.  
  135.     "_scal_d"  functions/routines :
  136.  
  137.         Scale a Sequence by a "real" value.
  138.  
  139.         This functions/routines need to be used when going back and forth
  140.         between domains, since a direct FFT followed by an inverse FFT
  141.         scale the original array by a factor of N (1D), N1*N2 (2D) or
  142.         N1*N2*N3 (3D).
  143.  
  144.         
  145. PARALLEL routines:
  146.         Only the 2D and the 3D routines have been parallelized.
  147.  
  148. MEMORY allocation:
  149.         LIBFFT gives you a control on the amount of scratch arrays 
  150.         allocated to a maximum value that can be adjusted at runtime 
  151.         through an environment variable (FFT_MAXMEM). The following command:
  152.         % setenv FFT_MAXMEM 16777216
  153.         will increase this Maximum size to 16 MegaBytes. This may be usefull
  154.         to increase performances when running large FFTs like 250x300x256.
  155.         The default value for FFT_MAXMEM is 10^6 Bytes.
  156.  
  157.  
  158. ---------------------------
  159. |   COMPLEX <---> COMPLEX |
  160. ---------------------------
  161.  
  162.  1 DIMENSION :
  163.  =============
  164.  
  165.    C :
  166.         complex *cfft1di( int n, complex *workspace);
  167.  
  168.         zomplex *zfft1di( int n, zomplex *workspace);
  169.         
  170.         int cfft1d ( int job, int n, complex *sequence,
  171.              int inc, complex *workspace);
  172.  
  173.         int zfft1d ( int job, int n, zomplex *sequence,
  174.              int inc, zomplex *workspace);
  175.  
  176.     Fortran:
  177.         subroutine CFFT1DI( n, workspace)
  178.         integer        n
  179.         complex        workspace(*)
  180.  
  181.         subroutine ZFFT1DI( n, workspace)
  182.         integer        n
  183.         double complex  workspace(*)
  184.  
  185.         subroutine CFFT1D( job, n, sequence, inc, workspace)
  186.         integer job, n, inc
  187.         complex    sequence(*), workspace(*)
  188.  
  189.         subroutine ZFFT1D( job, n, sequence, inc, workspace)
  190.         integer job, n, inc
  191.         double complex    sequence(*), workspace(*)
  192.  
  193.  
  194.  "cfft1di", "zfft1di"  IMPORTANT:
  195.   
  196.         The "workspace" array must be at least (N + 15) elements long.
  197.         This array requires  N elements for the Sines and Cosines
  198.                 + 15 elements for the factor decomposition of N.
  199.  
  200.  "cfft1d", "zfft1d"
  201.  
  202.     job        = -1 or 1  for forward or backward Fast Fourier Transform
  203.     n        = number of elements in the 1 D FFT
  204.     sequence    = array containing the elements for the FFT
  205.     inc        = increment between two consecutive elements of the sequence
  206.     workspace    = Array containing the Sines/Cosines and factorization of N
  207.             workspace needs to be initialized by a call to _fft1di
  208.  
  209.  
  210. ---------------------------
  211. |   COMPLEX <---> COMPLEX |
  212. ---------------------------
  213.  
  214.  2 DIMENSIONS :
  215.  ==============
  216.  
  217.    C :
  218.         complex    *cfft2di( int n1, int n2, complex *workspace)
  219.  
  220.         zomplex    *zfft2di( int n1, int n2, zomplex *workspace)
  221.  
  222.         int        cfft2d( int job, int n1, int n2, complex *sequence, 
  223.                 int lda, complex *workspace)
  224.  
  225.         int        zfft2d( int job, int n1, int n2, zomplex *sequence, 
  226.                 int lda, zomplex *workspace)
  227.  
  228.     Fortran:
  229.         subroutine    CFFT2DI( n1, n2, workspace)
  230.         integer        n1, n2
  231.         complex        workspace(*)
  232.  
  233.         subroutine    ZFFT2DI( n1, n2, workspace)
  234.         integer        n1, n2
  235.         complex        workspace(*)
  236.  
  237.         subroutine CFFT2D( job, n1, n2, sequence, lda, workspace)
  238.         integer        job, n1, n2, lda
  239.         complex        sequence(lda,*), workspace(*)
  240.  
  241.         subroutine ZFFT2D( job, n1, n2, sequence, lda, workspace)
  242.         integer        job, n1, n2, lda
  243.         double complex        sequence(lda,*), workspace(*)
  244.  
  245.  
  246.   "cfft2di", "zfft2di"  IMPORTANT:  
  247.  
  248.     the "workspace" array must be at least 
  249.     (N1 + 15 + N2 + 15) elements long for Zfft2di and Cfft2di
  250.     this array needs N1+N2 elements for the Sines and Cosines
  251.         + 2 x 15 elements for the factorization of N1 and N2
  252.  
  253.  
  254.   "cfft2d", "zfft2d" :
  255.  
  256.     job        = -1 or 1  for forward or backward Fast Fourier Transform
  257.     n1        = number of elements in the First Dimension of the Sequence
  258.     n2        = number of elements in the Second Dimension of the Sequence
  259.     sequence    = Array containing the (lda*n2) elements for the FFT
  260.     lda        = Leading dimension for the Array, (lda >= n1)
  261.     workspace    = Array containing the Sin/Cos and factorization of N1 and N2
  262.             workspace needs to be initialized by a call to _fft2di
  263.  
  264. ---------------------------
  265. |   COMPLEX <---> COMPLEX |
  266. ---------------------------
  267.  
  268.  3 DIMENSIONS :
  269.  ==============
  270.  
  271.  
  272.    C :
  273.     complex    *cfft3di( int n1, int n2, int n3, complex *workspace);
  274.  
  275.     zomplex    *zfft3di( int n1, int n2, int n3, zomplex *workspace);
  276.  
  277.     int    cfft3d( int job, int n1, int n2, int n3, complex *sequence, 
  278.                 int ld1, int ld2, complex *workspace)
  279.  
  280.     int    zfft3d( int job, int n1, int n2, int n3, zomplex *sequence, 
  281.                 int ld1, int ld2, zomplex *workspace)
  282.  
  283.    Fortran:
  284.  
  285.     subroutine  CFFT3DI( n1, n2, n3, workspace)
  286.         integer        n1, n2, n3
  287.         complex        workspace(*)
  288.  
  289.     subroutine  ZFFT3DI( n1, n2, n3, workspace)
  290.         integer        n1, n2, n3
  291.         double complex    workspace(*)
  292.  
  293.     subroutine  CFFT3D( job, n1, n2, n3, sequence, ld1, ld2, workspace)
  294.         integer        job, n1, n2, ld1, ld2
  295.         complex        sequence(ld1,ld2,*), workspace(*)
  296.  
  297.     subroutine  ZFFT3D( job, n1, n2, n3, sequence, ld1, ld2, workspace)
  298.         integer        job, n1, n2, ld1, ld2
  299.         double complex    sequence(ld1,ld2,*), workspace(*)
  300.  
  301.  
  302.   "cfft3di", "zfft3di"  IMPORTANT:  
  303.  
  304.     the "workspace" array must be at least 
  305.     (N1 + 15 + N2 + 15 + N3 + 15 ) elements long for Zfft3di and Cfft3di
  306.     this array needs N1+N2+N3 elements for the Sines and Cosines
  307.         + 3 x 15 elements for the factorization of N1,N2 and N3.
  308.  
  309.  
  310.  "cfft3d", "zfft3d" :
  311.  
  312.     job        = -1 or 1  for forward or backward Fast Fourier Transform
  313.     n1,n2,n3    = number of elements in the First, Second and Third 
  314.             Dimension of the Sequence
  315.     sequence    = Array containing the (ld1*ld2*n3) elements for the FFT
  316.     ld1, ld2    = Leading dimensions for the Array, (ld1 >= n1 and ld2 >= n2)
  317.     workspace    = Array containing the Sin/Cos and factorization of N1, N2, N3
  318.             workspace needs to be initialized by a call to _fft3di
  319.  
  320. ---------------------------
  321. |   COMPLEX <---> COMPLEX |
  322. ---------------------------
  323.  
  324.  1 DIMENSION  :
  325.  ==============
  326.  
  327. Fortran:
  328.  
  329.   C
  330.   C    the N elements of sequence(j) for j = 1, ..., n
  331.   C        are the SUM for k = 1, ..., n of
  332.   C        sequence(k) * exp( job*iz*(j-1)*(k-1)*2*PI/N)
  333.   C            where iz = sqrt(-1)
  334.   C
  335.  
  336. C language:
  337.   /*
  338.     the N elements sequence[j] for j = 0, ..., (n-1)
  339.         are the SUM for k = 0, ..., (n-1) of
  340.         sequence[k] * exp( job*iz*j*k*2*PI / n)
  341.             where iz = sqrt( -1)
  342.   */
  343.  
  344.  2 DIMENSIONS :
  345.  ==============
  346.  
  347.   Fortran:
  348.     C    the N1xN2 elements of sequence(j,k) for j=1,...,n1 and k=1,...,n2
  349.     C      are the SUM for l = 1, ..., n1 and m = 1, ..., n2 of 
  350.     C        sequence(l,m) * exp( job*2*iz*PI*((j-1)*(l-1)/N1 + (k-1)*(m-1)/N2))
  351.     C            where iz = sqrt(-1)
  352.  
  353.   C language:
  354.     /*
  355.       the N1xN2 elements sequence[j+k*lda] for j=0,..,(n1-1) and k=0,..,(n2-1)
  356.         are the SUM for l = 0,..,(n1-1) and m = 0,..,(n2-1) of
  357.         sequence[l+m*lda] * exp( job*2*iz*PI*(j*l/N1 + k*m/N2))
  358.             where iz = sqrt( -1)
  359.     */
  360.  3 DIMENSIONS :
  361.  ==============
  362.   Fortran:
  363.  
  364.     C the N1xN2xN3 elements of sequence (j1,j2,j3) 
  365.     C   for j1 = 1, ..., n1, for j2 = 1, ..., n2   and   j3 = 1, ..., n3
  366.     C    are the SUM for k1, k2 and k3 of: 
  367.     C        sequence(k1,k2,k3) * exp( job*2*iz*PI*
  368.     C        ((j1-1)*(k1-1)/N1 + (j2-1)*(k2-1)/N2 + (j3-1)*(k3-1)/N3))
  369.     C                where iz = sqrt(-1)
  370.  
  371.   C language:
  372.     /*
  373.     the N1xN2xN3 elements "sequence[j0+j1*ld1+j2*ld2]"
  374.         are the SUM of
  375.         sequence[k0+k1*ld1+k2*ld2] 
  376.             * exp( job*2*iz*PI*(j1*k1/N1 + j2*k2/N2 + j3*k3/N3))
  377.             where iz = sqrt( -1)
  378.     */
  379.  
  380.  
  381. ---------------------------
  382. |   REAL  <--->  COMPLEX  |
  383. ---------------------------
  384.  
  385.  1 DIMENSION :
  386.  =============
  387.  
  388.    C :
  389.         float   *sfft1dui( int n, float *workspace);
  390.  
  391.         double  *dfft1dui( int n, double *workspace);
  392.         
  393.         int sfft1d ( int job, int n, float *sequence,
  394.              int inc, float *workspace);
  395.  
  396.         int dfft1d ( int job, int n, double *sequence,
  397.              int inc, double *workspace);
  398.  
  399.     Fortran:
  400.         subroutine SFFT1DUI( n, workspace)
  401.         integer        n
  402.         real        workspace(*)
  403.  
  404.         subroutine DFFT1DUI( n, workspace)
  405.         integer        n
  406.         double precision  workspace(*)
  407.  
  408.         subroutine SFFT1DU( job, n, sequence, inc, workspace)
  409.         integer job, n, inc
  410.         real    sequence(*), workspace(*)
  411.  
  412.         subroutine DFFT1DU( job, n, sequence, inc, workspace)
  413.         integer job, n, inc
  414.         double precision    sequence(*), workspace(*)
  415.  
  416.  
  417.  "sfft1dui", "dfft1dui"  IMPORTANT:
  418.   
  419.         The "workspace" array must be at least (N + 15) elements long.
  420.         This array requires  N elements for the Sines and Cosines
  421.                 + 15 elements for the factor decomposition of N.
  422.  
  423.  "sfft1du", "dfft1du"
  424.  
  425.     job        = -1 or 1  for forward or backward Fast Fourier Transform
  426.     n        = number of elements in the 1 D FFT
  427.     sequence    = array containing the elements for the FFT
  428.     inc        = increment between two consecutive elements of the sequence
  429.     workspace    = Array containing the Sines/Cosines and factorization of N
  430.             workspace needs to be initialized by a call to _fft1di
  431.  
  432.   NOTA BENE :  After the direct FFT, the sequence is 2*((n+2)/2) elements long.
  433.         The Fourier Transform is 1 or 2 elements LONGER than the 
  434.         original sequence.
  435.         PLEASE, pad your input array, so that no useful data is over 
  436.         written.
  437.  
  438. ---------------------------
  439. |   REAL  <--->  COMPLEX  |
  440. ---------------------------
  441.  
  442.  2 DIMENSIONS :
  443.  ==============
  444.  
  445.    C :
  446.         float     *sfft2dui( int n1, int n2, float  *workspace)
  447.  
  448.         double    *dfft2dui( int n1, int n2, double *workspace)
  449.  
  450.         int        sfft2du( int job, int n1, int n2, float  *sequence, 
  451.                 int lda, float  *workspace)
  452.  
  453.         int        dfft2du( int job, int n1, int n2, double *sequence, 
  454.                 int lda, double *workspace)
  455.  
  456.     Fortran:
  457.         subroutine    SFFT2DUI( n1, n2, workspace)
  458.         integer        n1, n2
  459.         real        workspace(*)
  460.  
  461.         subroutine    DFFT2DUI( n1, n2, workspace)
  462.         integer        n1, n2
  463.         double precision        workspace(*)
  464.  
  465.         subroutine SFFT2DU( job, n1, n2, sequence, lda, workspace)
  466.         integer        job, n1, n2, lda
  467.         real        sequence(lda,*), workspace(*)
  468.  
  469.         subroutine DFFT2DU( job, n1, n2, sequence, lda, workspace)
  470.         integer        job, n1, n2, lda
  471.         double precision        sequence(lda,*), workspace(*)
  472.  
  473.  
  474.   "sfft2dui", "dfft2dui"  IMPORTANT:  
  475.  
  476.     the "workspace" array must be at least 
  477.     (N1 + 15 + 2*N2 + 15 ) elements long for Zfft2di and Cfft2di
  478.     N1 elements are needed for the Sines and Cosines of the first dim.
  479.     But we need the (2 * N2 ) for the second dimension ( COMPLEX ).
  480.         + 2 x 15 elements for the factorization of N1 and N2
  481.  
  482.  
  483.   "sfft2du", "dfft2du" :
  484.  
  485.     job        = -1 or 1  for forward or backward Fast Fourier Transform
  486.     n1        = number of elements in the First Dimension of the Sequence
  487.     n2        = number of elements in the Second Dimension of the Sequence
  488.     sequence    = Array containing the (lda*n2) elements for the FFT
  489.     lda        = Leading dimension for the Array, (lda >= 2*((n1+2)/2) )
  490.     workspace    = Array containing the Sin/Cos and factorization of N1 and N2
  491.             workspace needs to be initialized by a call to _fft2di
  492.  
  493.   NOTA BENE :  During the direct FFT, after the FFT along the First dimension, 
  494.         the sequences are 2*((n1+2)/2) elements long.
  495.         PLEASE, make sure that   lda >= 2*((n1+2)/2)
  496.  
  497. ---------------------------
  498. |   REAL  <--->  COMPLEX  |
  499. ---------------------------
  500.  
  501.  3 DIMENSIONS :
  502.  ==============
  503.  
  504.  
  505.    C :
  506.     float     *sfft3dui( int n1, int n2, int n3, float  *workspace);
  507.  
  508.     double    *dfft3dui( int n1, int n2, int n3, double *workspace);
  509.  
  510.     int    sfft3du( int job, int n1, int n2, int n3, float  *sequence, 
  511.                 int ld1, int ld2, float  *workspace)
  512.  
  513.     int    dfft3du( int job, int n1, int n2, int n3, double *sequence, 
  514.                 int ld1, int ld2, double *workspace)
  515.  
  516.    Fortran:
  517.  
  518.     subroutine  SFFT3DUI( n1, n2, n3, workspace)
  519.         integer        n1, n2, n3
  520.         real        workspace(*)
  521.  
  522.     subroutine  DFFT3DUI( n1, n2, n3, workspace)
  523.         integer        n1, n2, n3
  524.         double precision    workspace(*)
  525.  
  526.     subroutine  SFFT3DU( job, n1, n2, n3, sequence, ld1, ld2, workspace)
  527.         integer        job, n1, n2, ld1, ld2
  528.         real        sequence(ld1,ld2,*), workspace(*)
  529.  
  530.     subroutine  DFFT3DU( job, n1, n2, n3, sequence, ld1, ld2, workspace)
  531.         integer        job, n1, n2, ld1, ld2
  532.         double precision    sequence(ld1,ld2,*), workspace(*)
  533.  
  534.  
  535.   "sfft3dui", "dfft3dui"  IMPORTANT:  
  536.  
  537.     the "workspace" array must be at least 
  538.     (N1 + 15 + 2*N2 + 15 + 2*N3 + 15 ) elements long for Zfft3di and Cfft3di
  539.     this array needs N1+2*N2+2*N3 elements for the Sines and Cosines
  540.         + 3 x 15 elements for the factorization of N1, N2 and N3.
  541.  
  542.  
  543.  "sfft3du", "dfft3du" :
  544.  
  545.     job        = -1 or 1  for forward or backward Fast Fourier Transform
  546.     n1,n2,n3    = number of elements in the First, Second and Third 
  547.             Dimension of the Sequence
  548.     sequence    = Array containing the (ld1*ld2*n3) elements for the FFT
  549.     ld1, ld2    = Leading dimensions for the Array, 
  550.             (ld1 >= 2*((n1+2)/2)  and ld2 >= n2)
  551.     workspace    = Array containing the Sin/Cos and factorization of N1, N2, N3
  552.             workspace needs to be initialized by a call to _fft3di
  553.  
  554.   NOTA BENE :  During the direct FFT, after the FFT along the First dimension, 
  555.         the sequences are 2*((n1+2)/2) elements long.
  556.         PLEASE, make sure that   ld1 >= 2*((n1+2)/2)
  557.  
  558. ---------------------------
  559. |   REAL  <--->  COMPLEX  |
  560. ---------------------------
  561.  
  562.  1 DIMENSION  :
  563.  ==============
  564.  
  565.     Given a sequence a(j), j = 0, ..., n-1 of  N real samples, 
  566.     and B(I) its (complex) FFT.
  567.  
  568.     Only the first half of the FFT need to be computed, the second half can be 
  569.     deduced by symmetry  :  B(I) = Conjugate { B(N-I) }.
  570.  
  571.     example n = 5 :
  572.     ---------------
  573.     a0        a1        a2        a3        a4
  574.     B0r,B0i        B1r,B1i    B2r,B2i        B3r,B3i    B4r,B4i
  575.  
  576.         B0i = 0            <--- Frequence 0
  577.         B4r = B1r, B4i = - B1i
  578.         B3r = B2r, B3i = - B2i
  579.  
  580.       the output from SFFT1DU or DFFT1DU will be:
  581.  
  582.     B0r,0.0        B1r,B1i    B2r,B2i
  583.  
  584.     example n = 4 :
  585.     ---------------
  586.     a0        a1        a2        a3
  587.     B0r,B0i        B1r,B1i    B2r,B2i        B3r,B3i
  588.  
  589.         B0i = 0                <--- Frequence 0
  590.         B2i = 0                <--- Nyquist's Frequence (N/2)
  591.         B3r = B1r, B3i = - B1i    <--- Symmetry due to real input
  592.  
  593.       the output from SFFT1DU or DFFT1DU will be:
  594.  
  595.     B0r,0.0        B1r,B1i    B2r,0.0
  596.  
  597. *****
  598.     In "sfft1du" and "dfft1du" (contrary to "sfft1d" and "dfft1d" ) no 
  599.     attempt is made to pack the result (we explicitly keep the zeroes).
  600.     The Fourier transform is therefore much simpler to interpret as a 
  601.     sequence of COMPLEX numbers. 
  602.     Formally, the definition "sfft1du" is the SAME as for "cfft1d".
  603.  
  604.  2 and 3 DIMENSIONS :
  605.  ====================
  606.     For a 2D the direct Fourier transform we compute :
  607.     first the    N2     ( REAL--->COMPLEX ) 1D FFTs of size N1,
  608.     then         (N1+2)/2     ( COMPLEX--->COMPLEX ) FFTs of size N2.
  609.  
  610.     For a 3D the direct Fourier transform we compute :
  611.     first the    N2xN3          (  REAL  --->COMPLEX ) 1D FFTs of size N1,
  612.     then         ((N1+2)/2)xN3   ( COMPLEX--->COMPLEX ) 1D FFTs of size N2, 
  613.     and finally  ((N1+2)/2)xN2   ( COMPLEX--->COMPLEX ) 1D FFTs of size N3.
  614.  
  615.     The computations are performed in the reverse order for the inverse FFTs.
  616.  
  617. -----------------------------------------
  618. |   REAL  <--->  COMPLEX             |
  619. -----------------------------------------
  620.  
  621.  2 DIMENSIONS :
  622.  ==============
  623.  
  624.     Original Real Array(lda,5) --- N1 = 4 and N2 = 5:
  625.     |--------|--------|--------|--------|--------|
  626.     | f(1,1) | f(1,2) | f(1,3) | f(1,4) | f(1,5) |
  627.     |--------|--------|--------|--------|--------|
  628.     | f(2,1) | f(2,2) | f(2,3) | f(2,4) | f(2,5) |
  629.     |--------|--------|--------|--------|--------|
  630.     | f(3,1) | f(3,2) | f(3,3) | f(3,4) | f(3,5) |
  631.     |--------|--------|--------|--------|--------|
  632.     | f(4,1) | f(4,2) | f(4,3) | f(4,4) | f(4,5) |
  633.     |--------|--------|--------|--------|--------|
  634.  
  635.     After First Dimension Real FFT:
  636.     |--------|--------|--------|--------|--------|
  637.     | Re1(1) | Re1(2) | Re1(3) | Re1(4) | Re1(5) |
  638.     |  0.0   |  0.0   |  0.0   |  0.0   |  0.0   |
  639.     |--------|--------|--------|--------|--------|
  640.     | Re2(1) | Re2(2) | Re2(3) | Re2(4) | Re2(5) |
  641.     | Im2(1) | Im2(2) | Im2(3) | Im2(4) | Im2(5) |
  642.     |--------|--------|--------|--------|--------|
  643.     | Re3(1) | Re3(2) | Re3(3) | Re3(4) | Re3(5) |
  644.     |  0.0   |  0.0   |  0.0   |  0.0   |  0.0   |
  645.     |--------|--------|--------|--------|--------|
  646.  
  647.     After Second Dimension FFT:
  648.     |--------|--------|--------|--------|--------|
  649.     |Fr(1,1) |Fr(1,2) |Fr(1,3) |Fr(1,4) |Fr(1,5) |
  650.     |  0.0   |  0.0   |  0.0   |  0.0   |  0.0   |
  651.     |--------|--------|--------|--------|--------|
  652.     |Fr(2,1) |Fr(2,2) |Fr(2,3) |Fr(2,4) |Fr(2,5) |
  653.     |Fi(2,1) |Fi(2,2) |Fi(2,3) |Fi(2,4) |Fi(2,5) |
  654.     |--------|--------|--------|--------|--------|
  655.     |Fr(3,1) |Fr(3,2) |Fr(3,3) |Fr(3,4) |Fr(3,5) |
  656.     |  0.0   |  0.0   |  0.0   |  0.0   |  0.0   |
  657.     |--------|--------|--------|--------|--------|
  658.  
  659.     The Real ---> Complex FFTs store only the first half of the Fourier
  660.     transform. 
  661.     For simplicity, the values corresponding to Frequency 0,and Nyquist 
  662.     (N/2 when is N is an even number), are treated as Complex numbers
  663.     even though their imaginary part is ZERO.
  664.     When processing the Lines of the matrix for the FFTs along the second
  665.     dimension, all lines are then considered as Complex arrays.
  666.  
  667.     NOTE: when using a matrix MxN, depending on M (odd or even) the result 
  668.     is and (M+1)xN or (M+2)xN matrix. 
  669.  
  670. ---------------------------------
  671. |   Fourier Transforms PRODUCT    |
  672. ---------------------------------
  673.  
  674.     It is well known that the convolution of two sequences in the space domain
  675.     is equivalent to the product of their Fourier transforms.
  676.     For convenience, the following modules perform a "product" of 2 Fourier
  677.     Transforms, one called "y", the second "filter".
  678.  
  679.  
  680.   C :
  681.   
  682.     void    cprod1d( int n, complex *y, int incy, complex *filter, int incx);
  683.  
  684.     void    zprod1d( int n, complex *y, int incy, complex *filter, int incx);
  685.  
  686.     void    sprod1du( int n, complex *y, int incy, complex *filter, int incx);
  687.  
  688.     void    dprod1du( int n, complex *y, int incy, complex *filter, int incx);
  689.  
  690.     void    cprod2d( int n1, int n2, complex *y,int ldy, 
  691.             complex *filter,int ldx);
  692.  
  693.     void    zprod2d( int n1, int n2, zomplex *y,int ldy, 
  694.             zomplex *filter,int ldx);
  695.  
  696.     void    sprod2du( int n1, int n2, float *y, int ldy, 
  697.             float *filter, int ldx);
  698.  
  699.     void    dprod2du( int n1, int n2, double *y, int ldy, 
  700.             double *filter,int ldx);
  701.  
  702.     void    cprod3d( int n1, int n2, int n3, complex *y, int ldy1, int ldy2, 
  703.             complex *filter, int ldx1, int ldx2);
  704.  
  705.     void    zprod3d( int n1, int n2, int n3, complex *y, int ldy1, int ldy2, 
  706.             complex *filter, int ldx1, int ldx2);
  707.  
  708.     void    sprod3du( int n1, int n2, int n3, complex *y, int ldy1, int ldy2, 
  709.             complex *filter, int ldx1, int ldx2);
  710.  
  711.     void    dprod3du( int n1, int n2, int n3, complex *y, int ldy1, int ldy2, 
  712.             complex *filter, int ldx1, int ldx2);
  713.  
  714.  
  715. ---------------------------------
  716. |   Fourier Transforms PRODUCT    |
  717. ---------------------------------
  718.  
  719.   Fortran:
  720.  
  721.  
  722.     subroutine CPROD1D( n, y, inca, filter, incf)
  723.         integer n, inca, incf
  724.         complex    y(*), filter(*)
  725.  
  726.     subroutine ZPROD1D( n, y, inca, filter, incf)
  727.         integer n, inca, incf
  728.         double complex    y(*), filter(*)
  729.  
  730.     subroutine SPROD1DU( n, y, inca, filter, incf)
  731.         integer n, inca, incf
  732.         real    y(*), filter(*)
  733.  
  734.     subroutine DPROD1DU( n, y, inca, filter, incf)
  735.         integer n, inca, incf
  736.         double precision    y(*), filter(*)
  737.  
  738.  
  739.     subroutine CPROD2D( n1, n2, y, lda, filter, ldf)
  740.         integer n1, n2, lda, ldf
  741.         complex    y(lda,*), filter(ldf,*)
  742.  
  743.     subroutine ZPROD2D( n1, n2, y, lda, filter, ldf)
  744.         integer n1, n2, lda, ldf
  745.         double complex    y(lda,*), filter(ldf,*)
  746.  
  747.     subroutine SPROD2DU( n1, n2, y, lda, filter, ldf)
  748.         integer n1, n2, lda, ldf
  749.         real    y(lda,*), filter(ldf,*)
  750.  
  751.     subroutine DPROD2DU( n1, n2, y, lda, filter, ldf)
  752.         integer n1, n2, lda, ldf
  753.         double precision    y(lda,*), filter(ldf,*)
  754.  
  755.  
  756.     subroutine CPROD3D( n1, n2, n3, y, ld1, ld2, filter, ldf1,ldf2)
  757.         integer n1, n2, n3, ld1, ld2, ldf1, ldf2
  758.         complex    y(ld1,ld2,*), filter(ldf1,ldf2,*)
  759.  
  760.     subroutine ZPROD3D( n1, n2, n3, y, ld1, ld2, filter, ldf1,ldf2)
  761.         integer n1, n2, n3, ld1, ld2, ldf1, ldf2
  762.         double complex    y(ld1,ld2,*), filter(ldf1,ldf2,*)
  763.  
  764.     subroutine SPROD3DU( n1, n2, n3, y, ld1,ld2, filter, ldf1,ldf2)
  765.         integer n1, n2, n3, ld1, ld2, ldf1, ldf2
  766.         real    y(ld1,ld2,*), filter(ldf1,ldf2,*)
  767.  
  768.     subroutine DPROD3DU( n1, n2, n3, y, ld1,ld2, filter, ldf1,ldf2)
  769.         integer n1, n2, n3, ld1, ld2, ldf1, ldf2
  770.         double precision    y(ld1,ld2,*), filter(ldf1,ldf2,*)
  771.  
  772.  
  773. -----------------------------------------
  774. |   Fourier Transforms      SCALING    |
  775. -----------------------------------------
  776.  
  777. The Fourier transforms of this library are not normalized, therefore after a 
  778. double call ( direct_transform, inverse_transform), the original sequence 
  779. is scaled by the total number of samples.
  780. To get absolute values, you need to scale back the result by 
  781.     alpha = 1/n_total 
  782.             Where  
  783.     n_total = N            ... for 1D FFTs,
  784.     n_total = N1 x N2        ... for 2D FFTs,
  785.     n_total = N1 x N2 x N3        ... for 3D FFTs.
  786.  
  787.   C:
  788.  
  789.     void cscal1d( int n, float alpha, complex y, int inc);
  790.     void zscal1d( int n, double alpha, zomplex y, int inc);
  791.     void sscal1d( int n, float alpha, float y, int inc);
  792.     void dscal1d( int n, double alpha, double y, int inc);
  793.  
  794.     void cscal2d( int nx, int ny, float alpha, complex y, int ld);
  795.     void zscal2d( int nx, int ny, double alpha, zomplex y, int ld);
  796.     void sscal2d( int nx, int ny, float alpha, float y, int ld);
  797.     void dscal2d( int nx, int ny, double alpha, double y, int ld);
  798.  
  799.     void cscal3d( int nx, int ny, int nz, float alpha, 
  800.             complex y, int ld1,int ld2);
  801.     void zscal3d( int nx, int ny, int nz, double alpha, 
  802.             zomplex y,int ld1,int ld2);
  803.     void sscal3d( int nx, int ny, int nz, float alpha, 
  804.             float y, int ld1,int ld2);
  805.     void dscal3d( int nx, int ny, int nz, double alpha, 
  806.             double y, int ld1,int ld2);
  807.  
  808.  
  809.  
  810.   Fortran:  
  811.  
  812.     subroutine CSCAL1D( n, alpha, sequence, inc)
  813.         integer    n, inc
  814.         complex    sequence(*)
  815.         real    alpha
  816.  
  817.     subroutine ZSCAL1D( n, alpha, sequence, inc)
  818.         integer        n, inc
  819.         double complex    sequence(*)
  820.         double precision    alpha
  821.  
  822.     subroutine SSCAL1D( n, alpha, sequence, inc)
  823.         integer n, inc
  824.         real    sequence(*), alpha
  825.  
  826.     subroutine DSCAL1D( n, alpha, sequence, inc)
  827.         integer        n, inc
  828.         double precision    sequence(*), alpha
  829.  
  830.  
  831.     subroutine CSCAL2D( nx, ny, alpha, array, ld)
  832.         integer        nx ,ny, ld
  833.         complex        array(ld,*)
  834.         real        alpha
  835.  
  836.     subroutine ZSCAL2D( nx, ny, alpha, array, ld)
  837.         integer        nx ,ny, ld
  838.         double complex    array(ld,*)
  839.         double precision  alpha
  840.     
  841.     subroutine SSCAL2D( nx, ny, alpha, array, ld)
  842.         integer        nx ,ny, ld
  843.         real        array(ld,*), alpha
  844.  
  845.     subroutine DSCAL2D( nx, ny, alpha, array, ld)
  846.         integer        nx ,ny, ld
  847.         double precision    array(ld,*), alpha
  848.  
  849.  
  850.     subroutine CSCAL3D( nx, ny, nz, alpha, array, ld1, ld2)
  851.         integer        nx ,ny, nz, ld1, ld2
  852.         complex        array(ld1,ld2,*)
  853.         real        alpha
  854.  
  855.     subroutine ZSCAL3D( nx, ny, nz, alpha, array, ld1, ld2)
  856.         integer        nx ,ny, nz, ld1, ld2
  857.         double complex    array(ld1,ld2,*)
  858.         double precision    alpha
  859.     
  860.     subroutine SSCAL3D( nx, ny, nz, alpha, array, ld1, ld2)
  861.         integer        nx ,ny, nz, ld1, ld2
  862.         real        array(ld1,ld2,*), alpha
  863.  
  864.     subroutine DSCAL3D( nx, ny, nz, alpha, array, ld1, ld2)
  865.         integer        nx ,ny, nz, ld1, ld2
  866.         double precision    array(ld1,ld2,*), alpha
  867.  
  868. -----------------------------------------
  869. |   REAL  <--->  COMPLEX  *** PACKED    |
  870. -----------------------------------------
  871.  
  872.  
  873.  
  874.     The modules described in the following pages, the PACKED FFTs, are
  875.     kept for compatibility with the previous release.
  876.  
  877.     Various users experience proved these modules to be difficult to use,
  878.  
  879.     the "not_so_packed" modules described above are to be prefered !
  880.     Both versions "packed" and "not_so_packed" have equivalent performances.
  881.  
  882.  
  883.  
  884.     OLD "packed" versions        NEWER, SIMPLER, BETTER versions
  885.  
  886.     SFFT1DI, SFFT1D        --->    SFFT1DUI, SFFT1DU
  887.     DFFT1DI, DFFT1D        --->    DFFT1DUI, DFFT1DU
  888.  
  889.     SFFT2DI, SFFT2D        --->    SFFT2DUI, SFFT2DU
  890.     DFFT2DI, DFFT2D        --->    DFFT2DUI, DFFT2DU
  891.  
  892.     SFFT3DI, SFFT3D        --->    SFFT3DUI, SFFT3DU
  893.     DFFT3DI, DFFT3D        --->    DFFT3DUI, DFFT3DU
  894.  
  895.  
  896.  
  897. -----------------------------------------
  898. |   REAL  <--->  COMPLEX  *** PACKED    |
  899. -----------------------------------------
  900.  
  901.  1 DIMENSION :
  902.  =============
  903.  
  904.    C :
  905.         float   *sfft1di( int n, float *workspace);
  906.  
  907.         double  *dfft1di( int n, double *workspace);
  908.         
  909.         int sfft1d ( int job, int n, float *sequence,
  910.              int inc, float *workspace);
  911.  
  912.         int dfft1d ( int job, int n, double *sequence,
  913.              int inc, double *workspace);
  914.  
  915.     Fortran:
  916.         subroutine SFFT1DI( n, workspace)
  917.         integer        n
  918.         real        workspace(*)
  919.  
  920.         subroutine DFFT1DI( n, workspace)
  921.         integer        n
  922.         double precision  workspace(*)
  923.  
  924.         subroutine SFFT1D( job, n, sequence, inc, workspace)
  925.         integer job, n, inc
  926.         real    sequence(*), workspace(*)
  927.  
  928.         subroutine DFFT1D( job, n, sequence, inc, workspace)
  929.         integer job, n, inc
  930.         double precision    sequence(*), workspace(*)
  931.  
  932.  
  933.  "sfft1di", "dfft1di"  IMPORTANT:
  934.   
  935.         The "workspace" array must be at least (N + 15) elements long.
  936.         This array requires  N elements for the Sines and Cosines
  937.                 + 15 elements for the factor decomposition of N.
  938.  
  939.  "sfft1d", "dfft1d"
  940.  
  941.     job        = -1 or 1  for forward or backward Fast Fourier Transform
  942.     n        = number of elements in the 1 D FFT
  943.     sequence    = array containing the elements for the FFT
  944.     inc        = increment between two consecutive elements of the sequence
  945.     workspace    = Array containing the Sines/Cosines and factorization of N
  946.             workspace needs to be initialized by a call to _fft1di
  947.  
  948.   NOTA BENE :   The packed FFTs are done "in place"
  949.  
  950. -----------------------------------------
  951. |   REAL  <--->  COMPLEX  *** PACKED    |
  952. -----------------------------------------
  953.  
  954.  1 DIMENSION :
  955.  =============
  956.  
  957. "csfft1d" or "zdfft1d":
  958.  
  959.   C:
  960.           int csfft1d( int job,int n,float *sequence,
  961.                 int iof,int inc,float *workspace)
  962.  
  963.           int zdfft1d( int job,int n,double *sequence,
  964.                 int iof,int inc,double *workspace)
  965.  
  966.  
  967.   Fortran:
  968.         subroutine csfft1d( job, n, sequence, iof, inc, workspace)
  969.         integer job, n, iof, inc
  970.         complex    sequence(*), workspace(*)
  971.  
  972.         subroutine zdfft1d( job, n, sequence, iof, inc, workspace)
  973.         integer job, n, iof, inc
  974.         double complex    sequence(*), workspace(*)
  975.  
  976.  
  977.     job        = -1 or 1  for forward or backward Fast Fourier Transform
  978.     n        = number of complex elements of the 1 D FFT
  979.     sequence    = array containing the REAL values of a complex sequence
  980.     iof        = offset between the Imaginary and the Real value 
  981.             of a complex number
  982.     inc        = increment between two consecutive Real values of a sequence
  983.     workspace    = Array containing the Sines/Cosines and factorization of N
  984.             workspace needs to be initialized by a call to _fft1di
  985.  
  986. IMPORTANT: Theses routines must be used when the Imaginary part of a 
  987.         complex number does not follow immediately its Real part. They are
  988.         more general than CFFT1D and ZFFT1D.
  989.         call cfft1d( job, n, sequence, inc, workspace)
  990.         is equivalent to:
  991.         call csfft1d(( job, n, sequence, 1, 2*inc, workspace)
  992.  
  993.         CSFFT1D's workspace needs to be initialized by SFFT1DI.
  994.         ZDFFT1D's workspace needs to be initialized by ZFFT1DI.
  995.  
  996.  
  997. -----------------------------------------
  998. |   REAL  <--->  COMPLEX  *** PACKED    |
  999. -----------------------------------------
  1000.  
  1001.  2 DIMENSIONS :
  1002.  ==============
  1003.  
  1004.    C :
  1005.         float     *sfft2di( int n1, int n2, float  *workspace)
  1006.  
  1007.         double    *dfft2di( int n1, int n2, double *workspace)
  1008.  
  1009.         int        sfft2d( int job, int n1, int n2, float  *sequence, 
  1010.                 int lda, float  *workspace)
  1011.  
  1012.         int        dfft2d( int job, int n1, int n2, double *sequence, 
  1013.                 int lda, double *workspace)
  1014.  
  1015.     Fortran:
  1016.         subroutine    SFFT2DI( n1, n2, workspace)
  1017.         integer        n1, n2
  1018.         real        workspace(*)
  1019.  
  1020.         subroutine    DFFT2DI( n1, n2, workspace)
  1021.         integer        n1, n2
  1022.         double precision        workspace(*)
  1023.  
  1024.         subroutine SFFT2D( job, n1, n2, sequence, lda, workspace)
  1025.         integer        job, n1, n2, lda
  1026.         real        sequence(lda,*), workspace(*)
  1027.  
  1028.         subroutine DFFT2D( job, n1, n2, sequence, lda, workspace)
  1029.         integer        job, n1, n2, lda
  1030.         double precision        sequence(lda,*), workspace(*)
  1031.  
  1032.  
  1033.   "sfft2di", "dfft2di"  IMPORTANT:  
  1034.  
  1035.     the "workspace" array must be at least 
  1036.     (N1 + 15 + 3 * (N2 + 15) ) elements long for Zfft2di and Cfft2di
  1037.     N1 elements are needed for the Sines and Cosines of the first dim.
  1038.     But we need the (2 * N2) for the second dimension ( COMPLEX ),
  1039.     and (1 * N2) for the second dimension ( REAL ),
  1040.         + 4 x 15 elements for the factorization of N1 and N2
  1041.  
  1042.  
  1043.   "sfft2d", "dfft2d" :
  1044.  
  1045.     job        = -1 or 1  for forward or backward Fast Fourier Transform
  1046.     n1        = number of elements in the First Dimension of the Sequence
  1047.     n2        = number of elements in the Second Dimension of the Sequence
  1048.     sequence    = Array containing the (lda*n2) elements for the FFT
  1049.     lda        = Leading dimension for the Array, (lda >= n1 )
  1050.     workspace    = Array containing the Sin/Cos and factorization of N1 and N2
  1051.             workspace needs to be initialized by a call to _fft2di
  1052.  
  1053.   NOTA BENE :   The packed FFTs are done "in place"
  1054.  
  1055. -----------------------------------------
  1056. |   REAL  <--->  COMPLEX  *** PACKED    |
  1057. -----------------------------------------
  1058.  
  1059.  3 DIMENSIONS :
  1060.  ==============
  1061.  
  1062.  
  1063.    C :
  1064.     float     *sfft3di( int n1, int n2, int n3, float  *workspace);
  1065.  
  1066.     double    *dfft3di( int n1, int n2, int n3, double *workspace);
  1067.  
  1068.     int    sfft3d( int job, int n1, int n2, int n3, float  *sequence, 
  1069.                 int ld1, int ld2, float  *workspace)
  1070.  
  1071.     int    dfft3d( int job, int n1, int n2, int n3, double *sequence, 
  1072.                 int ld1, int ld2, double *workspace)
  1073.  
  1074.    Fortran:
  1075.  
  1076.     subroutine  SFFT3DI( n1, n2, n3, workspace)
  1077.         integer        n1, n2, n3
  1078.         real        workspace(*)
  1079.  
  1080.     subroutine  DFFT3DI( n1, n2, n3, workspace)
  1081.         integer        n1, n2, n3
  1082.         double precision    workspace(*)
  1083.  
  1084.     subroutine  SFFT3D( job, n1, n2, n3, sequence, ld1, ld2, workspace)
  1085.         integer        job, n1, n2, ld1, ld2
  1086.         real        sequence(ld1,ld2,*), workspace(*)
  1087.  
  1088.     subroutine  DFFT3D( job, n1, n2, n3, sequence, ld1, ld2, workspace)
  1089.         integer        job, n1, n2, ld1, ld2
  1090.         double precision    sequence(ld1,ld2,*), workspace(*)
  1091.  
  1092.  
  1093.   "sfft3di", "dfft3di"  IMPORTANT:  
  1094.  
  1095.     the "workspace" array must be at least 
  1096.     (N1 + 15 + 3 * (N2 + 15) + 3*(N3 + 15) ) elements long for Zfft3di 
  1097.     and Cfft3di this array needs 
  1098.     N1+N2+N3 elements for the Sines and Cosines (REAL)
  1099.     2*N2 + 2*N3  for the sines (COMPLEX)
  1100.         + 7 x 15 elements for the factorization of N1, N2 and N3.
  1101.  
  1102.  
  1103.  "sfft3d", "dfft3d" :
  1104.  
  1105.     job        = -1 or 1  for forward or backward Fast Fourier Transform
  1106.     n1,n2,n3    = number of elements in the First, Second and Third 
  1107.             Dimension of the Sequence
  1108.     sequence    = Array containing the (ld1*ld2*n3) elements for the FFT
  1109.     ld1, ld2    = Leading dimensions for the Array, 
  1110.             (ld1 >= 2*((n1+2)/2)  and ld2 >= n2)
  1111.     workspace    = Array containing the Sin/Cos and factorization of N1, N2, N3
  1112.             workspace needs to be initialized by a call to _fft3di
  1113.  
  1114.   NOTA BENE :   The packed FFTs are done "in place"
  1115.  
  1116.  
  1117. -----------------------------------------
  1118. |   REAL  <--->  COMPLEX  *** PACKED    |
  1119. -----------------------------------------
  1120.  
  1121.  1 DIMENSION  :
  1122.  ==============
  1123.  
  1124.     Given a sequence a(j), j = 0, ..., n-1 of  N real samples, 
  1125.     and B(I) its (complex) FFT.
  1126.  
  1127.     Only the first half of the FFT need to be computed, the second half can be 
  1128.     deduced by symmetry  :  B(I) = Conjugate { B(N-I) }.
  1129.  
  1130.     example n = 5 :
  1131.     ---------------
  1132.     a0        a1        a2        a3        a4
  1133.     B0r,B0i        B1r,B1i    B2r,B2i        B3r,B3i    B4r,B4i
  1134.  
  1135.         B0i = 0            <--- Frequence 0
  1136.         B4r = B1r, B4i = - B1i
  1137.         B3r = B2r, B3i = - B2i
  1138.  
  1139.       the output from SFFT1D or DFFT1D will be:
  1140.  
  1141.     B0r,    B1r,B1i,    B2r,B2i
  1142.  
  1143.     example n = 4 :
  1144.     ---------------
  1145.     a0        a1        a2        a3
  1146.     B0r,B0i        B1r,B1i    B2r,B2i        B3r,B3i
  1147.  
  1148.         B0i = 0            <--- Frequence 0
  1149.         B2i = 0            <--- Nyquist's Frequence (N/2)
  1150.         B3r = B1r, B3i = - B1i
  1151.  
  1152.       Once PACKED, the output from SFFT1D or DFFT1D will be:
  1153.  
  1154.     B0r,    B1r,B1i,    B2r
  1155.  
  1156. *****
  1157.     PLEASE note that the first value (and the last when N is even )is real,
  1158.     while the others are complex numbers.
  1159.  
  1160. -----------------------------------------
  1161. |   REAL  <--->  COMPLEX  *** PACKED    |
  1162. -----------------------------------------
  1163.  
  1164.  2 DIMENSIONS :
  1165.  ==============
  1166.  
  1167.     Original Real Array(lda,5) --- N1 = 4 and N2 = 5:
  1168.     |--------|--------|--------|--------|--------|
  1169.     | f(1,1) | f(1,2) | f(1,3) | f(1,4) | f(1,5) |
  1170.     |--------|--------|--------|--------|--------|
  1171.     | f(2,1) | f(2,2) | f(2,3) | f(2,4) | f(2,5) |
  1172.     |--------|--------|--------|--------|--------|
  1173.     | f(3,1) | f(3,2) | f(3,3) | f(3,4) | f(3,5) |
  1174.     |--------|--------|--------|--------|--------|
  1175.     | f(4,1) | f(4,2) | f(4,3) | f(4,4) | f(4,5) |
  1176.     |--------|--------|--------|--------|--------|
  1177.  
  1178.     After First Dimension Real FFT:
  1179.     |--------|--------|--------|--------|--------|
  1180.     | Re1(1) | Re1(2) | Re1(3) | Re1(4) | Re1(5) |
  1181.     |--------|--------|--------|--------|--------|
  1182.     | Re2(1) | Re2(2) | Re2(3) | Re2(4) | Re2(5) |
  1183.     | Im2(1) | Im2(2) | Im2(3) | Im2(4) | Im2(5) |
  1184.     |--------|--------|--------|--------|--------|
  1185.     | Re3(1) | Re3(2) | Re3(3) | Re3(4) | Re3(5) |
  1186.     |--------|--------|--------|--------|--------|
  1187.  
  1188.     After Second Dimension FFT:
  1189.     |--------|-----------------|-----------------|
  1190.     | F(1,1) | F(1,2)   F(1,3) | F(1,4)   F(1,5) |
  1191.     |--------|--------|--------|--------|--------|
  1192.     | F(2,1) | F(2,2) | F(2,3) | F(2,4) | F(2,5) |
  1193.     | F(3,1) | F(3,2) | F(3,3) | F(3,4) | F(3,5) |
  1194.     |--------|-----------------|-----------------|
  1195.     | F(4,1) | F(4,2)   F(4,3) | F(4,4)   F(4,5) |
  1196.     |--------|--------|--------|--------|--------|
  1197.  
  1198.     The FIRST and LAST line FFT is a REAl FFT (SFFT1D) --- N=5 and INC=lda.
  1199.     Note that F(1,2) and F(1,3) are respectively the real and imaginary part of
  1200.     the same complex number.
  1201.     call SFFT1D( -1, 5, array(1,1), lda, workspace)
  1202.     call SFFT1D( -1, 5, array(4,1), lda, workspace)
  1203.  
  1204.     Lines 2 and 3 use a Complex FFT (CSFFT1D) --- N=5 , IOF=1, INC=lda.
  1205.     In this case CSFFT1D has to be called since lda may not be an even number.
  1206.     call CSFFT1D( -1, 5, array(2,1), 1, lda, workspace)
  1207.  
  1208. -----------------------------------------
  1209. |   REAL  <--->  COMPLEX  *** PACKED    |
  1210. -----------------------------------------
  1211.  
  1212.  2 DIMENSIONS :
  1213.  ==============
  1214.  
  1215.     Real(8x8) <--> Complex FFT OUTPUT:
  1216.  
  1217. Let's start with a real Array F(0:7,0:7):
  1218.  
  1219. |-----------------------------------------------------------------------|
  1220. | F(0,0) | F(0,1) | F(0,2) | F(0,3) | F(0,4) | F(0,5) | F(0,6) | F(0,7) |
  1221. |________|________|________|________|________|________|________|________|
  1222. | F(1,0) | F(1,1) | F(1,2) | F(1,3) | F(1,4) | F(1,5) | F(1,6) | F(1,7) |
  1223. |________|________|________|________|________|________|________|________|
  1224. | F(2,0) | F(2,1) | F(2,2) | F(2,3) | F(2,4) | F(2,5) | F(2,6) | F(2,7) |
  1225. |________|________|________|________|________|________|________|________|
  1226. | F(3,0) | F(3,1) | F(3,2) | F(3,3) | F(3,4) | F(3,5) | F(3,6) | F(3,7) |
  1227. |________|________|________|________|________|________|________|________|
  1228. | F(4,0) | F(4,1) | F(4,2) | F(4,3) | F(4,4) | F(4,5) | F(4,6) | F(4,7) |
  1229. |________|________|________|________|________|________|________|________|
  1230. | F(5,0) | F(5,1) | F(5,2) | F(5,3) | F(5,4) | F(5,5) | F(5,6) | F(5,7) |
  1231. |________|________|________|________|________|________|________|________|
  1232. | F(6,0) | F(6,1) | F(6,2) | F(6,3) | F(6,4) | F(6,5) | F(6,6) | F(6,7) |
  1233. |________|________|________|________|________|________|________|________|
  1234. | F(7,0) | F(7,1) | F(7,2) | F(7,3) | F(7,4) | F(7,5) | F(7,6) | F(7,7) |
  1235. |________|________|________|________|________|________|________|________|
  1236.  
  1237. Now After the First FFT (along columns) ...
  1238.  
  1239. |-----------------------------------------------------------------------|
  1240. | R(0,0) | R(0,1) | R(0,2) | R(0,3) | R(0,4) | R(0,5) | R(0,6) | R(0,7) |
  1241. |________|________|________|________|________|________|________|________|
  1242. | R(1,0) | R(1,1) | R(1,2) | R(1,3) | R(1,4) | R(1,5) | R(1,6) | R(1,7) |
  1243. |________|________|________|________|________|________|________|________|
  1244. | I(1,0) | I(1,1) | I(1,2) | I(1,3) | I(1,4) | I(1,5) | I(1,6) | I(1,7) |
  1245. |________|________|________|________|________|________|________|________|
  1246. | R(2,0) | R(2,1) | R(2,2) | R(2,3) | R(2,4) | R(2,5) | R(2,6) | R(2,7) |
  1247. |________|________|________|________|________|________|________|________|
  1248. | I(2,0) | I(2,1) | I(2,2) | I(2,3) | I(2,4) | I(2,5) | I(2,6) | I(2,7) |
  1249. |________|________|________|________|________|________|________|________|
  1250. | R(3,0) | R(3,1) | R(3,2) | R(3,3) | R(3,4) | R(3,5) | R(3,6) | R(3,7) |
  1251. |________|________|________|________|________|________|________|________|
  1252. | I(3,0) | I(3,1) | I(3,2) | I(3,3) | I(3,4) | I(3,5) | I(3,6) | I(3,7) |
  1253. |________|________|________|________|________|________|________|________|
  1254. | R(4,0) | R(4,1) | R(4,2) | R(4,3) | R(4,4) | R(4,5) | R(4,6) | R(4,7) |
  1255. |________|________|________|________|________|________|________|________|
  1256.  
  1257. Where R(i,j) represents the Real part of the 1D FFTs, and I(i,j) the
  1258. imaginary part.
  1259.  
  1260. After the second Direction FFTs (along lines) the final result is:
  1261.  
  1262. |-----------------------------------------------------------------------|
  1263. | R(0,0) | R(0,1) | I(0,1) | R(0,2) | I(0,2) | R(0,3) | I(0,3) | R(0,4) |
  1264. |________|________|________|________|________|________|________|________|
  1265. | R(1,0) | R(1,1) | R(1,2) | R(1,3) | R(1,4) | R(1,5) | R(1,6) | R(1,7) |
  1266. |________|________|________|________|________|________|________|________|
  1267. | I(1,0) | I(1,1) | I(1,2) | I(1,3) | I(1,4) | I(1,5) | I(1,6) | I(1,7) |
  1268. |________|________|________|________|________|________|________|________|
  1269. | R(2,0) | R(2,1) | R(2,2) | R(2,3) | R(2,4) | R(2,5) | R(2,6) | R(2,7) |
  1270. |________|________|________|________|________|________|________|________|
  1271. | I(2,0) | I(2,1) | I(2,2) | I(2,3) | I(2,4) | I(2,5) | I(2,6) | I(2,7) |
  1272. |________|________|________|________|________|________|________|________|
  1273. | R(3,0) | R(3,1) | R(3,2) | R(3,3) | R(3,4) | R(3,5) | R(3,6) | R(3,7) |
  1274. |________|________|________|________|________|________|________|________|
  1275. | I(3,0) | I(3,1) | I(3,2) | I(3,3) | I(3,4) | I(3,5) | I(3,6) | I(3,7) |
  1276. |________|________|________|________|________|________|________|________|
  1277. | R(4,0) | R(4,1) | I(4,1) | R(4,2) | I(4,2) | R(4,3) | I(4,3) | R(4,4) |
  1278. |________|________|________|________|________|________|________|________|
  1279.  
  1280. Suppose You want to UNPACK this result in (5x8) Complex Matrix to get a 
  1281. result equivalent to SPROD2D , you  would fill it up like this:
  1282.  
  1283. |-----------------------------------------------------------------------|
  1284. | R(0,0) | R(0,1) | R(0,2) | R(0,3) | R(0,4) | r(0,3) | r(0,2) | r(0,1) |
  1285. | 0.     | I(0,1) | I(0,2) | I(0,3) | 0.     |-i(0,3) |-i(0,2) |-i(0,1) |
  1286. |________|________|________|________|________|________|________|________|
  1287. | R(1,0) | R(1,1) | R(1,2) | R(1,3) | R(1,4) | R(1,5) | R(1,6) | R(1,7) |
  1288. | I(1,0) | I(1,1) | I(1,2) | I(1,3) | I(1,4) | I(1,5) | I(1,6) | I(1,7) |
  1289. |________|________|________|________|________|________|________|________|
  1290. | R(2,0) | R(2,1) | R(2,2) | R(2,3) | R(2,4) | R(2,5) | R(2,6) | R(2,7) |
  1291. | I(2,0) | I(2,1) | I(2,2) | I(2,3) | I(2,4) | I(2,5) | I(2,6) | I(2,7) |
  1292. |________|________|________|________|________|________|________|________|
  1293. | R(3,0) | R(3,1) | R(3,2) | R(3,3) | R(3,4) | R(3,5) | R(3,6) | R(3,7) |
  1294. | I(3,0) | I(3,1) | I(3,2) | I(3,3) | I(3,4) | I(3,5) | I(3,6) | I(3,7) |
  1295. |________|________|________|________|________|________|________|________|
  1296. | R(4,0) | R(4,1) | R(4,2) | R(4,3) | R(4,4) | r(4,3) | r(4,2) | r(4,1) |
  1297. | 0.     | I(4,1) | I(4,2) | I(4,3) | 0.     |-i(4,3) |-i(4,2) |-i(4,1) |
  1298. |________|________|________|________|________|________|________|________|
  1299.  
  1300. The UPPER case (R & I) represent original values, while the lower case (r & i) 
  1301. represent redundant data that can be deduced from symmetries.
  1302.  
  1303. |------------------|
  1304. |   3 Dimensions   |
  1305. |------------------|
  1306.  
  1307. Complex <--> Complex FFT OUTPUT:
  1308.  
  1309. Real <--> Complex FFT OUTPUT:
  1310.  
  1311.     1)- 2 D FFTs are first performed along the First and Second direction 
  1312.         ( see above the Real <---> Complex FFT OUTPUT in the 2 D case) for
  1313.         each XY slice.
  1314.     2)- The third dimension FFT is then performed.
  1315.  
  1316.     There is no real problem for the Complex YZ slices. However the First
  1317.     (and Last if N1 is even) YZ slices need to be carefully explained.
  1318.  
  1319.  Suppose a N1x5x4 FFT, the first XY slice, indexes (1,j,k) will look like 
  1320. the following after the First 2 Dimensions FFTs:
  1321.     |----------|---------------------|---------------------|
  1322.     | F(1,1,1) | F(1,2,1)   F(1,3,1) | F(1,4,1)   F(1,5,1) |
  1323.     |----------|---------------------|---------------------|
  1324.     | F(1,1,2) | F(1,2,2)   F(1,3,2) | F(1,4,2)   F(1,5,2) |
  1325.     |----------|---------------------|---------------------|
  1326.     | F(1,1,3) | F(1,2,3)   F(1,3,3) | F(1,4,3)   F(1,5,3) |
  1327.     |----------|---------------------|---------------------|
  1328.     | F(1,1,4) | F(1,2,4)   F(1,3,4) | F(1,4,4)   F(1,5,4) |
  1329.     |----------|---------------------|---------------------|
  1330.  
  1331. after the third Dimension FFT, the array wil finally be:
  1332.     |----------|---------------------|---------------------|
  1333.     | G(1,1,1) | G(1,2,1)   G(1,3,1) | G(1,4,1)   G(1,5,1) |
  1334.     |----------|---------------------|---------------------|
  1335.     | G(1,1,2) | G(1,2,2)   G(1,3,2) | G(1,4,2)   G(1,5,2) |
  1336.     |          |---------------------|---------------------|
  1337.     | G(1,1,3) | G(1,2,3)   G(1,3,3) | G(1,4,3)   G(1,5,3) |
  1338.     |----------|---------------------|---------------------|
  1339.     | G(1,1,4) | G(1,2,4)   G(1,3,4) | G(1,4,4)   G(1,5,4) |
  1340.     |----------|---------------------|---------------------|
  1341.  
  1342. This is actually achieved through 1 call to SFFT1D and 2 Calls to CSFFT1D:
  1343.  
  1344.     call SFFT1D( -1, n3, F(1,1,1), ld1*ld2, workspace)
  1345.     call CSFFT1D( -1, n3, F(1,2,1), ld1, ld1*ld2, workspace)
  1346.     call CSFFT1D( -1, n3, F(1,4,1), ld1, ld1*ld2, workspace)
  1347.