comparison src/fftw-3.3.3/doc/html/Reversing-array-dimensions.html @ 95:89f5e221ed7b

Add FFTW3
author Chris Cannam <cannam@all-day-breakfast.com>
date Wed, 20 Mar 2013 15:35:50 +0000
parents
children
comparison
equal deleted inserted replaced
94:d278df1123f9 95:89f5e221ed7b
1 <html lang="en">
2 <head>
3 <title>Reversing array dimensions - FFTW 3.3.3</title>
4 <meta http-equiv="Content-Type" content="text/html">
5 <meta name="description" content="FFTW 3.3.3">
6 <meta name="generator" content="makeinfo 4.13">
7 <link title="Top" rel="start" href="index.html#Top">
8 <link rel="up" href="Calling-FFTW-from-Modern-Fortran.html#Calling-FFTW-from-Modern-Fortran" title="Calling FFTW from Modern Fortran">
9 <link rel="prev" href="Overview-of-Fortran-interface.html#Overview-of-Fortran-interface" title="Overview of Fortran interface">
10 <link rel="next" href="FFTW-Fortran-type-reference.html#FFTW-Fortran-type-reference" title="FFTW Fortran type reference">
11 <link href="http://www.gnu.org/software/texinfo/" rel="generator-home" title="Texinfo Homepage">
12 <!--
13 This manual is for FFTW
14 (version 3.3.3, 25 November 2012).
15
16 Copyright (C) 2003 Matteo Frigo.
17
18 Copyright (C) 2003 Massachusetts Institute of Technology.
19
20 Permission is granted to make and distribute verbatim copies of
21 this manual provided the copyright notice and this permission
22 notice are preserved on all copies.
23
24 Permission is granted to copy and distribute modified versions of
25 this manual under the conditions for verbatim copying, provided
26 that the entire resulting derived work is distributed under the
27 terms of a permission notice identical to this one.
28
29 Permission is granted to copy and distribute translations of this
30 manual into another language, under the above conditions for
31 modified versions, except that this permission notice may be
32 stated in a translation approved by the Free Software Foundation.
33 -->
34 <meta http-equiv="Content-Style-Type" content="text/css">
35 <style type="text/css"><!--
36 pre.display { font-family:inherit }
37 pre.format { font-family:inherit }
38 pre.smalldisplay { font-family:inherit; font-size:smaller }
39 pre.smallformat { font-family:inherit; font-size:smaller }
40 pre.smallexample { font-size:smaller }
41 pre.smalllisp { font-size:smaller }
42 span.sc { font-variant:small-caps }
43 span.roman { font-family:serif; font-weight:normal; }
44 span.sansserif { font-family:sans-serif; font-weight:normal; }
45 --></style>
46 </head>
47 <body>
48 <div class="node">
49 <a name="Reversing-array-dimensions"></a>
50 <p>
51 Next:&nbsp;<a rel="next" accesskey="n" href="FFTW-Fortran-type-reference.html#FFTW-Fortran-type-reference">FFTW Fortran type reference</a>,
52 Previous:&nbsp;<a rel="previous" accesskey="p" href="Overview-of-Fortran-interface.html#Overview-of-Fortran-interface">Overview of Fortran interface</a>,
53 Up:&nbsp;<a rel="up" accesskey="u" href="Calling-FFTW-from-Modern-Fortran.html#Calling-FFTW-from-Modern-Fortran">Calling FFTW from Modern Fortran</a>
54 <hr>
55 </div>
56
57 <h3 class="section">7.2 Reversing array dimensions</h3>
58
59 <p><a name="index-row_002dmajor-517"></a><a name="index-column_002dmajor-518"></a>A minor annoyance in calling FFTW from Fortran is that FFTW's array
60 dimensions are defined in the C convention (row-major order), while
61 Fortran's array dimensions are the opposite convention (column-major
62 order). See <a href="Multi_002ddimensional-Array-Format.html#Multi_002ddimensional-Array-Format">Multi-dimensional Array Format</a>. This is just a
63 bookkeeping difference, with no effect on performance. The only
64 consequence of this is that, whenever you create an FFTW plan for a
65 multi-dimensional transform, you must always <em>reverse the
66 ordering of the dimensions</em>.
67
68 <p>For example, consider the three-dimensional (L&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;N) arrays:
69
70 <pre class="example"> complex(C_DOUBLE_COMPLEX), dimension(L,M,N) :: in, out
71 </pre>
72 <p>To plan a DFT for these arrays using <code>fftw_plan_dft_3d</code>, you could do:
73
74 <p><a name="index-fftw_005fplan_005fdft_005f3d-519"></a>
75 <pre class="example"> plan = fftw_plan_dft_3d(N,M,L, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
76 </pre>
77 <p>That is, from FFTW's perspective this is a N&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;L array.
78 <em>No data transposition need occur</em>, as this is <em>only
79 notation</em>. Similarly, to use the more generic routine
80 <code>fftw_plan_dft</code> with the same arrays, you could do:
81
82 <pre class="example"> integer(C_INT), dimension(3) :: n = [N,M,L]
83 plan = fftw_plan_dft_3d(3, n, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
84 </pre>
85 <p>Note, by the way, that this is different from the legacy Fortran
86 interface (see <a href="Fortran_002dinterface-routines.html#Fortran_002dinterface-routines">Fortran-interface routines</a>), which automatically
87 reverses the order of the array dimension for you. Here, you are
88 calling the C interface directly, so there is no &ldquo;translation&rdquo; layer.
89
90 <p><a name="index-r2c_002fc2r-multi_002ddimensional-array-format-520"></a>An important thing to keep in mind is the implication of this for
91 multidimensional real-to-complex transforms (see <a href="Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data">Multi-Dimensional DFTs of Real Data</a>). In C, a multidimensional real-to-complex DFT
92 chops the last dimension roughly in half (N&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;L real input
93 goes to N&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;L/2+1 complex output). In Fortran, because
94 the array dimension notation is reversed, the <em>first</em> dimension of
95 the complex data is chopped roughly in half. For example consider the
96 &lsquo;<samp><span class="samp">r2c</span></samp>&rsquo; transform of L&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;N real input in Fortran:
97
98 <p><a name="index-fftw_005fplan_005fdft_005fr2c_005f3d-521"></a><a name="index-fftw_005fexecute_005fdft_005fr2c-522"></a>
99 <pre class="example"> type(C_PTR) :: plan
100 real(C_DOUBLE), dimension(L,M,N) :: in
101 complex(C_DOUBLE_COMPLEX), dimension(L/2+1,M,N) :: out
102 plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
103 ...
104 call fftw_execute_dft_r2c(plan, in, out)
105 </pre>
106 <p><a name="index-in_002dplace-523"></a><a name="index-padding-524"></a>Alternatively, for an in-place r2c transform, as described in the C
107 documentation we must <em>pad</em> the <em>first</em> dimension of the
108 real input with an extra two entries (which are ignored by FFTW) so as
109 to leave enough space for the complex output. The input is
110 <em>allocated</em> as a 2[L/2+1]&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;N array, even though only
111 L&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;N of it is actually used. In this example, we will
112 allocate the array as a pointer type, using &lsquo;<samp><span class="samp">fftw_alloc</span></samp>&rsquo; to
113 ensure aligned memory for maximum performance (see <a href="Allocating-aligned-memory-in-Fortran.html#Allocating-aligned-memory-in-Fortran">Allocating aligned memory in Fortran</a>); this also makes it easy to reference the
114 same memory as both a real array and a complex array.
115
116 <p><a name="index-fftw_005falloc_005fcomplex-525"></a><a name="index-c_005ff_005fpointer-526"></a>
117 <pre class="example"> real(C_DOUBLE), pointer :: in(:,:,:)
118 complex(C_DOUBLE_COMPLEX), pointer :: out(:,:,:)
119 type(C_PTR) :: plan, data
120 data = fftw_alloc_complex(int((L/2+1) * M * N, C_SIZE_T))
121 call c_f_pointer(data, in, [2*(L/2+1),M,N])
122 call c_f_pointer(data, out, [L/2+1,M,N])
123 plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
124 ...
125 call fftw_execute_dft_r2c(plan, in, out)
126 ...
127 call fftw_destroy_plan(plan)
128 call fftw_free(data)
129 </pre>
130 <!-- -->
131 </body></html>
132