File: Reversing-array-dimensions.html

package info (click to toggle)
fftw3 3.3.8-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster
  • size: 28,428 kB
  • sloc: ansic: 259,592; ml: 5,474; sh: 4,442; perl: 1,648; makefile: 1,156; fortran: 110
file content (184 lines) | stat: -rw-r--r-- 9,238 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
<!-- This manual is for FFTW
(version 3.3.8, 24 May 2018).

Copyright (C) 2003 Matteo Frigo.

Copyright (C) 2003 Massachusetts Institute of Technology.

Permission is granted to make and distribute verbatim copies of this
manual provided the copyright notice and this permission notice are
preserved on all copies.

Permission is granted to copy and distribute modified versions of this
manual under the conditions for verbatim copying, provided that the
entire resulting derived work is distributed under the terms of a
permission notice identical to this one.

Permission is granted to copy and distribute translations of this manual
into another language, under the above conditions for modified versions,
except that this permission notice may be stated in a translation
approved by the Free Software Foundation. -->
<!-- Created by GNU Texinfo 6.3, http://www.gnu.org/software/texinfo/ -->
<head>
<title>FFTW 3.3.8: Reversing array dimensions</title>

<meta name="description" content="FFTW 3.3.8: Reversing array dimensions">
<meta name="keywords" content="FFTW 3.3.8: Reversing array dimensions">
<meta name="resource-type" content="document">
<meta name="distribution" content="global">
<meta name="Generator" content="makeinfo">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<link href="index.html#Top" rel="start" title="Top">
<link href="Concept-Index.html#Concept-Index" rel="index" title="Concept Index">
<link href="index.html#SEC_Contents" rel="contents" title="Table of Contents">
<link href="Calling-FFTW-from-Modern-Fortran.html#Calling-FFTW-from-Modern-Fortran" rel="up" title="Calling FFTW from Modern Fortran">
<link href="FFTW-Fortran-type-reference.html#FFTW-Fortran-type-reference" rel="next" title="FFTW Fortran type reference">
<link href="Extended-and-quadruple-precision-in-Fortran.html#Extended-and-quadruple-precision-in-Fortran" rel="prev" title="Extended and quadruple precision in Fortran">
<style type="text/css">
<!--
a.summary-letter {text-decoration: none}
blockquote.indentedblock {margin-right: 0em}
blockquote.smallindentedblock {margin-right: 0em; font-size: smaller}
blockquote.smallquotation {font-size: smaller}
div.display {margin-left: 3.2em}
div.example {margin-left: 3.2em}
div.lisp {margin-left: 3.2em}
div.smalldisplay {margin-left: 3.2em}
div.smallexample {margin-left: 3.2em}
div.smalllisp {margin-left: 3.2em}
kbd {font-style: oblique}
pre.display {font-family: inherit}
pre.format {font-family: inherit}
pre.menu-comment {font-family: serif}
pre.menu-preformatted {font-family: serif}
pre.smalldisplay {font-family: inherit; font-size: smaller}
pre.smallexample {font-size: smaller}
pre.smallformat {font-family: inherit; font-size: smaller}
pre.smalllisp {font-size: smaller}
span.nolinebreak {white-space: nowrap}
span.roman {font-family: initial; font-weight: normal}
span.sansserif {font-family: sans-serif; font-weight: normal}
ul.no-bullet {list-style: none}
-->
</style>


</head>

<body lang="en">
<a name="Reversing-array-dimensions"></a>
<div class="header">
<p>
Next: <a href="FFTW-Fortran-type-reference.html#FFTW-Fortran-type-reference" accesskey="n" rel="next">FFTW Fortran type reference</a>, Previous: <a href="Overview-of-Fortran-interface.html#Overview-of-Fortran-interface" accesskey="p" rel="prev">Overview of Fortran interface</a>, Up: <a href="Calling-FFTW-from-Modern-Fortran.html#Calling-FFTW-from-Modern-Fortran" accesskey="u" rel="up">Calling FFTW from Modern Fortran</a> &nbsp; [<a href="index.html#SEC_Contents" title="Table of contents" rel="contents">Contents</a>][<a href="Concept-Index.html#Concept-Index" title="Index" rel="index">Index</a>]</p>
</div>
<hr>
<a name="Reversing-array-dimensions-1"></a>
<h3 class="section">7.2 Reversing array dimensions</h3>

<a name="index-row_002dmajor-6"></a>
<a name="index-column_002dmajor-1"></a>
<p>A minor annoyance in calling FFTW from Fortran is that FFTW&rsquo;s array
dimensions are defined in the C convention (row-major order), while
Fortran&rsquo;s array dimensions are the opposite convention (column-major
order). See <a href="Multi_002ddimensional-Array-Format.html#Multi_002ddimensional-Array-Format">Multi-dimensional Array Format</a>.  This is just a
bookkeeping difference, with no effect on performance.  The only
consequence of this is that, whenever you create an FFTW plan for a
multi-dimensional transform, you must always <em>reverse the
ordering of the dimensions</em>.
</p>
<p>For example, consider the three-dimensional (L&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;N
) arrays:
</p>
<div class="example">
<pre class="example">  complex(C_DOUBLE_COMPLEX), dimension(L,M,N) :: in, out
</pre></div>

<p>To plan a DFT for these arrays using <code>fftw_plan_dft_3d</code>, you could do:
</p>
<a name="index-fftw_005fplan_005fdft_005f3d-2"></a>
<div class="example">
<pre class="example">  plan = fftw_plan_dft_3d(N,M,L, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
</pre></div>

<p>That is, from FFTW&rsquo;s perspective this is a N&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;L
 array.
<em>No data transposition need occur</em>, as this is <em>only
notation</em>.  Similarly, to use the more generic routine
<code>fftw_plan_dft</code> with the same arrays, you could do:
</p>
<div class="example">
<pre class="example">  integer(C_INT), dimension(3) :: n = [N,M,L]
  plan = fftw_plan_dft_3d(3, n, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
</pre></div>

<p>Note, by the way, that this is different from the legacy Fortran
interface (see <a href="Fortran_002dinterface-routines.html#Fortran_002dinterface-routines">Fortran-interface routines</a>), which automatically
reverses the order of the array dimension for you.  Here, you are
calling the C interface directly, so there is no &ldquo;translation&rdquo; layer.
</p>
<a name="index-r2c_002fc2r-multi_002ddimensional-array-format-2"></a>
<p>An important thing to keep in mind is the implication of this for
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
chops the last dimension roughly in half (N&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;L
 real input
goes to N&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;L/2+1
 complex output).  In Fortran, because
the array dimension notation is reversed, the <em>first</em> dimension of
the complex data is chopped roughly in half.  For example consider the
&lsquo;<samp>r2c</samp>&rsquo; transform of L&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;N
 real input in Fortran:
</p>
<a name="index-fftw_005fplan_005fdft_005fr2c_005f3d-2"></a>
<a name="index-fftw_005fexecute_005fdft_005fr2c-1"></a>
<div class="example">
<pre class="example">  type(C_PTR) :: plan
  real(C_DOUBLE), dimension(L,M,N) :: in
  complex(C_DOUBLE_COMPLEX), dimension(L/2+1,M,N) :: out
  plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
  ...
  call fftw_execute_dft_r2c(plan, in, out)
</pre></div>

<a name="index-in_002dplace-9"></a>
<a name="index-padding-5"></a>
<p>Alternatively, for an in-place r2c transform, as described in the C
documentation we must <em>pad</em> the <em>first</em> dimension of the
real input with an extra two entries (which are ignored by FFTW) so as
to leave enough space for the complex output. The input is
<em>allocated</em> as a 2[L/2+1]&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;N
 array, even though only
L&nbsp;&times;&nbsp;M&nbsp;&times;&nbsp;N
 of it is actually used.  In this example, we will
allocate the array as a pointer type, using &lsquo;<samp>fftw_alloc</samp>&rsquo; to
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
same memory as both a real array and a complex array.
</p>
<a name="index-fftw_005falloc_005fcomplex-4"></a>
<a name="index-c_005ff_005fpointer"></a>
<div class="example">
<pre class="example">  real(C_DOUBLE), pointer :: in(:,:,:)
  complex(C_DOUBLE_COMPLEX), pointer :: out(:,:,:)
  type(C_PTR) :: plan, data
  data = fftw_alloc_complex(int((L/2+1) * M * N, C_SIZE_T))
  call c_f_pointer(data, in, [2*(L/2+1),M,N])
  call c_f_pointer(data, out, [L/2+1,M,N])
  plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
  ...
  call fftw_execute_dft_r2c(plan, in, out)
  ...
  call fftw_destroy_plan(plan)
  call fftw_free(data)
</pre></div>

<hr>
<div class="header">
<p>
Next: <a href="FFTW-Fortran-type-reference.html#FFTW-Fortran-type-reference" accesskey="n" rel="next">FFTW Fortran type reference</a>, Previous: <a href="Overview-of-Fortran-interface.html#Overview-of-Fortran-interface" accesskey="p" rel="prev">Overview of Fortran interface</a>, Up: <a href="Calling-FFTW-from-Modern-Fortran.html#Calling-FFTW-from-Modern-Fortran" accesskey="u" rel="up">Calling FFTW from Modern Fortran</a> &nbsp; [<a href="index.html#SEC_Contents" title="Table of contents" rel="contents">Contents</a>][<a href="Concept-Index.html#Concept-Index" title="Index" rel="index">Index</a>]</p>
</div>



</body>
</html>