forked from FFTW/fftw3
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmodern-fortran.texi
725 lines (617 loc) · 29.5 KB
/
modern-fortran.texi
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
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
@node Calling FFTW from Modern Fortran, Calling FFTW from Legacy Fortran, Distributed-memory FFTW with MPI, Top
@chapter Calling FFTW from Modern Fortran
@cindex Fortran interface
Fortran 2003 standardized ways for Fortran code to call C libraries,
and this allows us to support a direct translation of the FFTW C API
into Fortran. Compared to the legacy Fortran 77 interface
(@pxref{Calling FFTW from Legacy Fortran}), this direct interface
offers many advantages, especially compile-time type-checking and
aligned memory allocation. As of this writing, support for these C
interoperability features seems widespread, having been implemented in
nearly all major Fortran compilers (e.g. GNU, Intel, IBM,
Oracle/Solaris, Portland Group, NAG).
@cindex portability
This chapter documents that interface. For the most part, since this
interface allows Fortran to call the C interface directly, the usage
is identical to C translated to Fortran syntax. However, there are a
few subtle points such as memory allocation, wisdom, and data types
that deserve closer attention.
@menu
* Overview of Fortran interface::
* Reversing array dimensions::
* FFTW Fortran type reference::
* Plan execution in Fortran::
* Allocating aligned memory in Fortran::
* Accessing the wisdom API from Fortran::
* Defining an FFTW module::
@end menu
@c -------------------------------------------------------
@node Overview of Fortran interface, Reversing array dimensions, Calling FFTW from Modern Fortran, Calling FFTW from Modern Fortran
@section Overview of Fortran interface
FFTW provides a file @code{fftw3.f03} that defines Fortran 2003
interfaces for all of its C routines, except for the MPI routines
described elsewhere, which can be found in the same directory as
@code{fftw3.h} (the C header file). In any Fortran subroutine where
you want to use FFTW functions, you should begin with:
@cindex iso_c_binding
@example
use, intrinsic :: iso_c_binding
include 'fftw3.f03'
@end example
This includes the interface definitions and the standard
@code{iso_c_binding} module (which defines the equivalents of C
types). You can also put the FFTW functions into a module if you
prefer (@pxref{Defining an FFTW module}).
At this point, you can now call anything in the FFTW C interface
directly, almost exactly as in C other than minor changes in syntax.
For example:
@findex fftw_plan_dft_2d
@findex fftw_execute_dft
@findex fftw_destroy_plan
@example
type(C_PTR) :: plan
complex(C_DOUBLE_COMPLEX), dimension(1024,1000) :: in, out
plan = fftw_plan_dft_2d(1000,1024, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
...
call fftw_execute_dft(plan, in, out)
...
call fftw_destroy_plan(plan)
@end example
A few important things to keep in mind are:
@itemize @bullet
@item
@tindex fftw_complex
@ctindex C_PTR
@ctindex C_INT
@ctindex C_DOUBLE
@ctindex C_DOUBLE_COMPLEX
FFTW plans are @code{type(C_PTR)}. Other C types are mapped in the
obvious way via the @code{iso_c_binding} standard: @code{int} turns
into @code{integer(C_INT)}, @code{fftw_complex} turns into
@code{complex(C_DOUBLE_COMPLEX)}, @code{double} turns into
@code{real(C_DOUBLE)}, and so on. @xref{FFTW Fortran type reference}.
@item
Functions in C become functions in Fortran if they have a return value,
and subroutines in Fortran otherwise.
@item
The ordering of the Fortran array dimensions must be @emph{reversed}
when they are passed to the FFTW plan creation, thanks to differences
in array indexing conventions (@pxref{Multi-dimensional Array
Format}). This is @emph{unlike} the legacy Fortran interface
(@pxref{Fortran-interface routines}), which reversed the dimensions
for you. @xref{Reversing array dimensions}.
@item
@cindex alignment
@cindex SIMD
Using ordinary Fortran array declarations like this works, but may
yield suboptimal performance because the data may not be not aligned
to exploit SIMD instructions on modern proessors (@pxref{SIMD
alignment and fftw_malloc}). Better performance will often be obtained
by allocating with @samp{fftw_alloc}. @xref{Allocating aligned memory
in Fortran}.
@item
@findex fftw_execute
Similar to the legacy Fortran interface (@pxref{FFTW Execution in
Fortran}), we currently recommend @emph{not} using @code{fftw_execute}
but rather using the more specialized functions like
@code{fftw_execute_dft} (@pxref{New-array Execute Functions}).
However, you should execute the plan on the @code{same arrays} as the
ones for which you created the plan, unless you are especially
careful. @xref{Plan execution in Fortran}. To prevent
you from using @code{fftw_execute} by mistake, the @code{fftw3.f03}
file does not provide an @code{fftw_execute} interface declaration.
@item
@cindex flags
Multiple planner flags are combined with @code{ior} (equivalent to @samp{|} in C). e.g. @code{FFTW_MEASURE | FFTW_DESTROY_INPUT} becomes @code{ior(FFTW_MEASURE, FFTW_DESTROY_INPUT)}. (You can also use @samp{+} as long as you don't try to include a given flag more than once.)
@end itemize
@menu
* Extended and quadruple precision in Fortran::
@end menu
@node Extended and quadruple precision in Fortran, , Overview of Fortran interface, Overview of Fortran interface
@subsection Extended and quadruple precision in Fortran
@cindex precision
If FFTW is compiled in @code{long double} (extended) precision
(@pxref{Installation and Customization}), you may be able to call the
resulting @code{fftwl_} routines (@pxref{Precision}) from Fortran if
your compiler supports the @code{C_LONG_DOUBLE_COMPLEX} type code.
Because some Fortran compilers do not support
@code{C_LONG_DOUBLE_COMPLEX}, the @code{fftwl_} declarations are
segregated into a separate interface file @code{fftw3l.f03}, which you
should include @emph{in addition} to @code{fftw3.f03} (which declares
precision-independent @samp{FFTW_} constants):
@cindex iso_c_binding
@example
use, intrinsic :: iso_c_binding
include 'fftw3.f03'
include 'fftw3l.f03'
@end example
We also support using the nonstandard @code{__float128}
quadruple-precision type provided by recent versions of @code{gcc} on
32- and 64-bit x86 hardware (@pxref{Installation and Customization}),
using the corresponding @code{real(16)} and @code{complex(16)} types
supported by @code{gfortran}. The quadruple-precision @samp{fftwq_}
functions (@pxref{Precision}) are declared in a @code{fftw3q.f03}
interface file, which should be included in addition to
@code{fftw3.f03}, as above. You should also link with
@code{-lfftw3q -lquadmath -lm} as in C.
@c -------------------------------------------------------
@node Reversing array dimensions, FFTW Fortran type reference, Overview of Fortran interface, Calling FFTW from Modern Fortran
@section Reversing array dimensions
@cindex row-major
@cindex column-major
A minor annoyance in calling FFTW from Fortran is that FFTW's array
dimensions are defined in the C convention (row-major order), while
Fortran's array dimensions are the opposite convention (column-major
order). @xref{Multi-dimensional Array Format}. 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 @emph{reverse the
ordering of the dimensions}.
For example, consider the three-dimensional (@threedims{L,M,N}) arrays:
@example
complex(C_DOUBLE_COMPLEX), dimension(L,M,N) :: in, out
@end example
To plan a DFT for these arrays using @code{fftw_plan_dft_3d}, you could do:
@findex fftw_plan_dft_3d
@example
plan = fftw_plan_dft_3d(N,M,L, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
@end example
That is, from FFTW's perspective this is a @threedims{N,M,L} array.
@emph{No data transposition need occur}, as this is @emph{only
notation}. Similarly, to use the more generic routine
@code{fftw_plan_dft} with the same arrays, you could do:
@example
integer(C_INT), dimension(3) :: n = [N,M,L]
plan = fftw_plan_dft_3d(3, n, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
@end example
Note, by the way, that this is different from the legacy Fortran
interface (@pxref{Fortran-interface routines}), which automatically
reverses the order of the array dimension for you. Here, you are
calling the C interface directly, so there is no ``translation'' layer.
@cindex r2c/c2r multi-dimensional array format
An important thing to keep in mind is the implication of this for
multidimensional real-to-complex transforms (@pxref{Multi-Dimensional
DFTs of Real Data}). In C, a multidimensional real-to-complex DFT
chops the last dimension roughly in half (@threedims{N,M,L} real input
goes to @threedims{N,M,L/2+1} complex output). In Fortran, because
the array dimension notation is reversed, the @emph{first} dimension of
the complex data is chopped roughly in half. For example consider the
@samp{r2c} transform of @threedims{L,M,N} real input in Fortran:
@findex fftw_plan_dft_r2c_3d
@findex fftw_execute_dft_r2c
@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)
@end example
@cindex in-place
@cindex padding
Alternatively, for an in-place r2c transform, as described in the C
documentation we must @emph{pad} the @emph{first} 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
@emph{allocated} as a @threedims{2[L/2+1],M,N} array, even though only
@threedims{L,M,N} of it is actually used. In this example, we will
allocate the array as a pointer type, using @samp{fftw_alloc} to
ensure aligned memory for maximum performance (@pxref{Allocating
aligned memory in Fortran}); this also makes it easy to reference the
same memory as both a real array and a complex array.
@findex fftw_alloc_complex
@findex c_f_pointer
@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)
@end example
@c -------------------------------------------------------
@node FFTW Fortran type reference, Plan execution in Fortran, Reversing array dimensions, Calling FFTW from Modern Fortran
@section FFTW Fortran type reference
The following are the most important type correspondences between the
C interface and Fortran:
@itemize @bullet
@item
@tindex fftw_plan
Plans (@code{fftw_plan} and variants) are @code{type(C_PTR)} (i.e. an
opaque pointer).
@item
@tindex fftw_complex
@cindex precision
@ctindex C_DOUBLE
@ctindex C_FLOAT
@ctindex C_LONG_DOUBLE
@ctindex C_DOUBLE_COMPLEX
@ctindex C_FLOAT_COMPLEX
@ctindex C_LONG_DOUBLE_COMPLEX
The C floating-point types @code{double}, @code{float}, and @code{long
double} correspond to @code{real(C_DOUBLE)}, @code{real(C_FLOAT)}, and
@code{real(C_LONG_DOUBLE)}, respectively. The C complex types
@code{fftw_complex}, @code{fftwf_complex}, and @code{fftwl_complex}
correspond in Fortran to @code{complex(C_DOUBLE_COMPLEX)},
@code{complex(C_FLOAT_COMPLEX)}, and
@code{complex(C_LONG_DOUBLE_COMPLEX)}, respectively.
Just as in C
(@pxref{Precision}), the FFTW subroutines and types are prefixed with
@samp{fftw_}, @code{fftwf_}, and @code{fftwl_} for the different precisions, and link to different libraries (@code{-lfftw3}, @code{-lfftw3f}, and @code{-lfftw3l} on Unix), but use the @emph{same} include file @code{fftw3.f03} and the @emph{same} constants (all of which begin with @samp{FFTW_}). The exception is @code{long double} precision, for which you should @emph{also} include @code{fftw3l.f03} (@pxref{Extended and quadruple precision in Fortran}).
@item
@tindex ptrdiff_t
@ctindex C_INT
@ctindex C_INTPTR_T
@ctindex C_SIZE_T
@findex fftw_malloc
The C integer types @code{int} and @code{unsigned} (used for planner
flags) become @code{integer(C_INT)}. The C integer type @code{ptrdiff_t} (e.g. in the @ref{64-bit Guru Interface}) becomes @code{integer(C_INTPTR_T)}, and @code{size_t} (in @code{fftw_malloc} etc.) becomes @code{integer(C_SIZE_T)}.
@item
@tindex fftw_r2r_kind
@ctindex C_FFTW_R2R_KIND
The @code{fftw_r2r_kind} type (@pxref{Real-to-Real Transform Kinds})
becomes @code{integer(C_FFTW_R2R_KIND)}. The various constant values
of the C enumerated type (@code{FFTW_R2HC} etc.) become simply integer
constants of the same names in Fortran.
@item
@ctindex FFTW_DESTROY_INPUT
@cindex in-place
@findex fftw_flops
Numeric array pointer arguments (e.g. @code{double *})
become @code{dimension(*), intent(out)} arrays of the same type, or
@code{dimension(*), intent(in)} if they are pointers to constant data
(e.g. @code{const int *}). There are a few exceptions where numeric
pointers refer to scalar outputs (e.g. for @code{fftw_flops}), in which
case they are @code{intent(out)} scalar arguments in Fortran too.
For the new-array execute functions (@pxref{New-array Execute Functions}),
the input arrays are declared @code{dimension(*), intent(inout)}, since
they can be modified in the case of in-place or @code{FFTW_DESTROY_INPUT}
transforms.
@item
@findex fftw_alloc_real
@findex c_f_pointer
Pointer @emph{return} values (e.g @code{double *}) become
@code{type(C_PTR)}. (If they are pointers to arrays, as for
@code{fftw_alloc_real}, you can convert them back to Fortran array
pointers with the standard intrinsic function @code{c_f_pointer}.)
@item
@cindex guru interface
@tindex fftw_iodim
@tindex fftw_iodim64
@cindex 64-bit architecture
The @code{fftw_iodim} type in the guru interface (@pxref{Guru vector
and transform sizes}) becomes @code{type(fftw_iodim)} in Fortran, a
derived data type (the Fortran analogue of C's @code{struct}) with
three @code{integer(C_INT)} components: @code{n}, @code{is}, and
@code{os}, with the same meanings as in C. The @code{fftw_iodim64} type in the 64-bit guru interface (@pxref{64-bit Guru Interface}) is the same, except that its components are of type @code{integer(C_INTPTR_T)}.
@item
@ctindex C_FUNPTR
Using the wisdom import/export functions from Fortran is a bit tricky,
and is discussed in @ref{Accessing the wisdom API from Fortran}. In
brief, the @code{FILE *} arguments map to @code{type(C_PTR)}, @code{const char *} to @code{character(C_CHAR), dimension(*), intent(in)} (null-terminated!), and the generic read-char/write-char functions map to @code{type(C_FUNPTR)}.
@end itemize
@cindex portability
You may be wondering if you need to search-and-replace
@code{real(kind(0.0d0))} (or whatever your favorite Fortran spelling
of ``double precision'' is) with @code{real(C_DOUBLE)} everywhere in
your program, and similarly for @code{complex} and @code{integer}
types. The answer is no; you can still use your existing types. As
long as these types match their C counterparts, things should work
without a hitch. The worst that can happen, e.g. in the (unlikely)
event of a system where @code{real(kind(0.0d0))} is different from
@code{real(C_DOUBLE)}, is that the compiler will give you a
type-mismatch error. That is, if you don't use the
@code{iso_c_binding} kinds you need to accept at least the theoretical
possibility of having to change your code in response to compiler
errors on some future machine, but you don't need to worry about
silently compiling incorrect code that yields runtime errors.
@c -------------------------------------------------------
@node Plan execution in Fortran, Allocating aligned memory in Fortran, FFTW Fortran type reference, Calling FFTW from Modern Fortran
@section Plan execution in Fortran
In C, in order to use a plan, one normally calls @code{fftw_execute},
which executes the plan to perform the transform on the input/output
arrays passed when the plan was created (@pxref{Using Plans}). The
corresponding subroutine call in modern Fortran is:
@example
call fftw_execute(plan)
@end example
@findex fftw_execute
However, we have had reports that this causes problems with some
recent optimizing Fortran compilers. The problem is, because the
input/output arrays are not passed as explicit arguments to
@code{fftw_execute}, the semantics of Fortran (unlike C) allow the
compiler to assume that the input/output arrays are not changed by
@code{fftw_execute}. As a consequence, certain compilers end up
repositioning the call to @code{fftw_execute}, assuming incorrectly
that it does nothing to the arrays.
There are various workarounds to this, but the safest and simplest
thing is to not use @code{fftw_execute} in Fortran. Instead, use the
functions described in @ref{New-array Execute Functions}, which take
the input/output arrays as explicit arguments. For example, if the
plan is for a complex-data DFT and was created for the arrays
@code{in} and @code{out}, you would do:
@example
call fftw_execute_dft(plan, in, out)
@end example
@findex fftw_execute_dft
There are a few things to be careful of, however:
@itemize @bullet
@item
@findex fftw_execute_dft_r2c
@findex fftw_execute_dft_c2r
@findex fftw_execute_r2r
You must use the correct type of execute function, matching the way
the plan was created. Complex DFT plans should use
@code{fftw_execute_dft}, Real-input (r2c) DFT plans should use use
@code{fftw_execute_dft_r2c}, and real-output (c2r) DFT plans should
use @code{fftw_execute_dft_c2r}. The various r2r plans should use
@code{fftw_execute_r2r}. Fortunately, if you use the wrong one you
will get a compile-time type-mismatch error (unlike legacy Fortran).
@item
You should normally pass the same input/output arrays that were used when
creating the plan. This is always safe.
@item
@emph{If} you pass @emph{different} input/output arrays compared to
those used when creating the plan, you must abide by all the
restrictions of the new-array execute functions (@pxref{New-array
Execute Functions}). The most tricky of these is the
requirement that the new arrays have the same alignment as the
original arrays; the best (and possibly only) way to guarantee this
is to use the @samp{fftw_alloc} functions to allocate your arrays (@pxref{Allocating aligned memory in Fortran}). Alternatively, you can
use the @code{FFTW_UNALIGNED} flag when creating the
plan, in which case the plan does not depend on the alignment, but
this may sacrifice substantial performance on architectures (like x86)
with SIMD instructions (@pxref{SIMD alignment and fftw_malloc}).
@ctindex FFTW_UNALIGNED
@end itemize
@c -------------------------------------------------------
@node Allocating aligned memory in Fortran, Accessing the wisdom API from Fortran, Plan execution in Fortran, Calling FFTW from Modern Fortran
@section Allocating aligned memory in Fortran
@cindex alignment
@findex fftw_alloc_real
@findex fftw_alloc_complex
In order to obtain maximum performance in FFTW, you should store your
data in arrays that have been specially aligned in memory (@pxref{SIMD
alignment and fftw_malloc}). Enforcing alignment also permits you to
safely use the new-array execute functions (@pxref{New-array Execute
Functions}) to apply a given plan to more than one pair of in/out
arrays. Unfortunately, standard Fortran arrays do @emph{not} provide
any alignment guarantees. The @emph{only} way to allocate aligned
memory in standard Fortran is to allocate it with an external C
function, like the @code{fftw_alloc_real} and
@code{fftw_alloc_complex} functions. Fortunately, Fortran 2003 provides
a simple way to associate such allocated memory with a standard Fortran
array pointer that you can then use normally.
We therefore recommend allocating all your input/output arrays using
the following technique:
@enumerate
@item
Declare a @code{pointer}, @code{arr}, to your array of the desired type
and dimensions. For example, @code{real(C_DOUBLE), pointer :: a(:,:)}
for a 2d real array, or @code{complex(C_DOUBLE_COMPLEX), pointer ::
a(:,:,:)} for a 3d complex array.
@item
The number of elements to allocate must be an
@code{integer(C_SIZE_T)}. You can either declare a variable of this
type, e.g. @code{integer(C_SIZE_T) :: sz}, to store the number of
elements to allocate, or you can use the @code{int(..., C_SIZE_T)}
intrinsic function. e.g. set @code{sz = L * M * N} or use
@code{int(L * M * N, C_SIZE_T)} for an @threedims{L,M,N} array.
@item
Declare a @code{type(C_PTR) :: p} to hold the return value from
FFTW's allocation routine. Set @code{p = fftw_alloc_real(sz)} for a real array, or @code{p = fftw_alloc_complex(sz)} for a complex array.
@item
@findex c_f_pointer
Associate your pointer @code{arr} with the allocated memory @code{p}
using the standard @code{c_f_pointer} subroutine: @code{call
c_f_pointer(p, arr, [...dimensions...])}, where
@code{[...dimensions...])} are an array of the dimensions of the array
(in the usual Fortran order). e.g. @code{call c_f_pointer(p, arr,
[L,M,N])} for an @threedims{L,M,N} array. (Alternatively, you can
omit the dimensions argument if you specified the shape explicitly
when declaring @code{arr}.) You can now use @code{arr} as a usual
multidimensional array.
@item
When you are done using the array, deallocate the memory by @code{call
fftw_free(p)} on @code{p}.
@end enumerate
For example, here is how we would allocate an @twodims{L,M} 2d real array:
@example
real(C_DOUBLE), pointer :: arr(:,:)
type(C_PTR) :: p
p = fftw_alloc_real(int(L * M, C_SIZE_T))
call c_f_pointer(p, arr, [L,M])
@emph{...use arr and arr(i,j) as usual...}
call fftw_free(p)
@end example
and here is an @threedims{L,M,N} 3d complex array:
@example
complex(C_DOUBLE_COMPLEX), pointer :: arr(:,:,:)
type(C_PTR) :: p
p = fftw_alloc_complex(int(L * M * N, C_SIZE_T))
call c_f_pointer(p, arr, [L,M,N])
@emph{...use arr and arr(i,j,k) as usual...}
call fftw_free(p)
@end example
See @ref{Reversing array dimensions} for an example allocating a
single array and associating both real and complex array pointers with
it, for in-place real-to-complex transforms.
@c -------------------------------------------------------
@node Accessing the wisdom API from Fortran, Defining an FFTW module, Allocating aligned memory in Fortran, Calling FFTW from Modern Fortran
@section Accessing the wisdom API from Fortran
@cindex wisdom
@cindex saving plans to disk
As explained in @ref{Words of Wisdom-Saving Plans}, FFTW provides a
``wisdom'' API for saving plans to disk so that they can be recreated
quickly. The C API for exporting (@pxref{Wisdom Export}) and
importing (@pxref{Wisdom Import}) wisdom is somewhat tricky to use
from Fortran, however, because of differences in file I/O and string
types between C and Fortran.
@menu
* Wisdom File Export/Import from Fortran::
* Wisdom String Export/Import from Fortran::
* Wisdom Generic Export/Import from Fortran::
@end menu
@c =========>
@node Wisdom File Export/Import from Fortran, Wisdom String Export/Import from Fortran, Accessing the wisdom API from Fortran, Accessing the wisdom API from Fortran
@subsection Wisdom File Export/Import from Fortran
@findex fftw_import wisdom_from_filename
@findex fftw_export_wisdom_to_filename
The easiest way to export and import wisdom is to do so using
@code{fftw_export_wisdom_to_filename} and
@code{fftw_wisdom_from_filename}. The only trick is that these
require you to pass a C string, which is an array of type
@code{CHARACTER(C_CHAR)} that is terminated by @code{C_NULL_CHAR}.
You can call them like this:
@example
integer(C_INT) :: ret
ret = fftw_export_wisdom_to_filename(C_CHAR_'my_wisdom.dat' // C_NULL_CHAR)
if (ret .eq. 0) stop 'error exporting wisdom to file'
ret = fftw_import_wisdom_from_filename(C_CHAR_'my_wisdom.dat' // C_NULL_CHAR)
if (ret .eq. 0) stop 'error importing wisdom from file'
@end example
Note that prepending @samp{C_CHAR_} is needed to specify that the
literal string is of kind @code{C_CHAR}, and we null-terminate the
string by appending @samp{// C_NULL_CHAR}. These functions return an
@code{integer(C_INT)} (@code{ret}) which is @code{0} if an error
occurred during export/import and nonzero otherwise.
It is also possible to use the lower-level routines
@code{fftw_export_wisdom_to_file} and
@code{fftw_import_wisdom_from_file}, which accept parameters of the C
type @code{FILE*}, expressed in Fortran as @code{type(C_PTR)}.
However, you are then responsible for creating the @code{FILE*}
yourself. You can do this by using @code{iso_c_binding} to define
Fortran intefaces for the C library functions @code{fopen} and
@code{fclose}, which is a bit strange in Fortran but workable.
@c =========>
@node Wisdom String Export/Import from Fortran, Wisdom Generic Export/Import from Fortran, Wisdom File Export/Import from Fortran, Accessing the wisdom API from Fortran
@subsection Wisdom String Export/Import from Fortran
@findex fftw_export_wisdom_to_string
Dealing with FFTW's C string export/import is a bit more painful. In
particular, the @code{fftw_export_wisdom_to_string} function requires
you to deal with a dynamically allocated C string. To get its length,
you must define an interface to the C @code{strlen} function, and to
deallocate it you must define an interface to C @code{free}:
@example
use, intrinsic :: iso_c_binding
interface
integer(C_INT) function strlen(s) bind(C, name='strlen')
import
type(C_PTR), value :: s
end function strlen
subroutine free(p) bind(C, name='free')
import
type(C_PTR), value :: p
end subroutine free
end interface
@end example
Given these definitions, you can then export wisdom to a Fortran
character array:
@example
character(C_CHAR), pointer :: s(:)
integer(C_SIZE_T) :: slen
type(C_PTR) :: p
p = fftw_export_wisdom_to_string()
if (.not. c_associated(p)) stop 'error exporting wisdom'
slen = strlen(p)
call c_f_pointer(p, s, [slen+1])
...
call free(p)
@end example
@findex c_associated
@findex c_f_pointer
Note that @code{slen} is the length of the C string, but the length of
the array is @code{slen+1} because it includes the terminating null
character. (You can omit the @samp{+1} if you don't want Fortran to
know about the null character.) The standard @code{c_associated} function
checks whether @code{p} is a null pointer, which is returned by
@code{fftw_export_wisdom_to_string} if there was an error.
@findex fftw_import_wisdom_from_string
To import wisdom from a string, use
@code{fftw_import_wisdom_from_string} as usual; note that the argument
of this function must be a @code{character(C_CHAR)} that is terminated
by the @code{C_NULL_CHAR} character, like the @code{s} array above.
@c =========>
@node Wisdom Generic Export/Import from Fortran, , Wisdom String Export/Import from Fortran, Accessing the wisdom API from Fortran
@subsection Wisdom Generic Export/Import from Fortran
The most generic wisdom export/import functions allow you to provide
an arbitrary callback function to read/write one character at a time
in any way you want. However, your callback function must be written
in a special way, using the @code{bind(C)} attribute to be passed to a
C interface.
@findex fftw_export_wisdom
In particular, to call the generic wisdom export function
@code{fftw_export_wisdom}, you would write a callback subroutine of the form:
@example
subroutine my_write_char(c, p) bind(C)
use, intrinsic :: iso_c_binding
character(C_CHAR), value :: c
type(C_PTR), value :: p
@emph{...write c...}
end subroutine my_write_char
@end example
Given such a subroutine (along with the corresponding interface definition), you could then export wisdom using:
@findex c_funloc
@example
call fftw_export_wisdom(c_funloc(my_write_char), p)
@end example
@findex c_loc
@findex c_f_pointer
The standard @code{c_funloc} intrinsic converts a Fortran
@code{bind(C)} subroutine into a C function pointer. The parameter
@code{p} is a @code{type(C_PTR)} to any arbitrary data that you want
to pass to @code{my_write_char} (or @code{C_NULL_PTR} if none). (Note
that you can get a C pointer to Fortran data using the intrinsic
@code{c_loc}, and convert it back to a Fortran pointer in
@code{my_write_char} using @code{c_f_pointer}.)
Similarly, to use the generic @code{fftw_import_wisdom}, you would
define a callback function of the form:
@findex fftw_import_wisdom
@example
integer(C_INT) function my_read_char(p) bind(C)
use, intrinsic :: iso_c_binding
type(C_PTR), value :: p
character :: c
@emph{...read a character c...}
my_read_char = ichar(c, C_INT)
end function my_read_char
....
integer(C_INT) :: ret
ret = fftw_import_wisdom(c_funloc(my_read_char), p)
if (ret .eq. 0) stop 'error importing wisdom'
@end example
Your function can return @code{-1} if the end of the input is reached.
Again, @code{p} is an arbitrary @code{type(C_PTR} that is passed
through to your function. @code{fftw_import_wisdom} returns @code{0}
if an error occurred and nonzero otherwise.
@c -------------------------------------------------------
@node Defining an FFTW module, , Accessing the wisdom API from Fortran, Calling FFTW from Modern Fortran
@section Defining an FFTW module
Rather than using the @code{include} statement to include the
@code{fftw3.f03} interface file in any subroutine where you want to
use FFTW, you might prefer to define an FFTW Fortran module. FFTW
does not install itself as a module, primarily because
@code{fftw3.f03} can be shared between different Fortran compilers while
modules (in general) cannot. However, it is trivial to define your
own FFTW module if you want. Just create a file containing:
@example
module FFTW3
use, intrinsic :: iso_c_binding
include 'fftw3.f03'
end module
@end example
Compile this file into a module as usual for your compiler (e.g. with
@code{gfortran -c} you will get a file @code{fftw3.mod}). Now,
instead of @code{include 'fftw3.f03'}, whenever you want to use FFTW
routines you can just do:
@example
use FFTW3
@end example
as usual for Fortran modules. (You still need to link to the FFTW
library, of course.)