(git:08825b9)
Loading...
Searching...
No Matches
fftw3_lib.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
8 USE iso_c_binding, ONLY: c_associated, &
9 c_char, &
10 c_double, &
11 c_double_complex, &
12 c_int, &
13 c_ptr
14#if defined(__FFTW3)
15 USE iso_c_binding, ONLY: &
16 c_float, &
17 c_float_complex, &
18 c_funptr, &
19 c_int32_t, &
20 c_intptr_t, &
21 c_loc, &
22 c_null_char, &
23 c_size_t, c_f_pointer
24 USE mathconstants, ONLY: z_zero
25#endif
26 USE fft_kinds, ONLY: dp
27 USE fft_plan, ONLY: fft_plan_type
28
29!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
30
31#include "../../base/base_uses.f90"
32
33 IMPLICIT NONE
34 PRIVATE
35
38 PUBLIC :: fftw_alloc, fftw_dealloc
39
40 INTEGER, SAVE :: fftw_plan_type = -1
41
42#if defined(__FFTW3)
43#include "fftw3.f03"
44#endif
45
46 INTERFACE fftw_alloc
47 MODULE PROCEDURE :: fftw_alloc_complex_1d
48 MODULE PROCEDURE :: fftw_alloc_complex_2d
49 MODULE PROCEDURE :: fftw_alloc_complex_3d
50 END INTERFACE fftw_alloc
51
52 INTERFACE fftw_dealloc
53 MODULE PROCEDURE :: fftw_dealloc_complex_1d
54 MODULE PROCEDURE :: fftw_dealloc_complex_2d
55 MODULE PROCEDURE :: fftw_dealloc_complex_3d
56 END INTERFACE fftw_dealloc
57
58CONTAINS
59
60! Concatenate the components of the dimensions passed to this function to use it if FFTW3 is not used
61 SUBROUTINE fftw_alloc_complex_1d(array, n)
62 COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:), CONTIGUOUS, POINTER, INTENT(OUT) :: array
63 INTEGER, DIMENSION(1), INTENT(IN) :: n
64
65#if defined(__FFTW3)
66 TYPE(c_ptr) :: data_ptr
67 data_ptr = fftw_alloc_complex(int(product(n), kind=c_size_t))
68 CALL c_f_pointer(data_ptr, array, n)
69#else
70! Just allocate the array
71 ALLOCATE (array(n(1)))
72#endif
73
74 END SUBROUTINE fftw_alloc_complex_1d
75
76 SUBROUTINE fftw_dealloc_complex_1d(array)
77 COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:), CONTIGUOUS, POINTER, INTENT(INOUT) :: array
78
79#if defined(__FFTW3)
80 CALL fftw_free(c_loc(array))
81 NULLIFY (array)
82#else
83! Just deallocate the array
84 DEALLOCATE (array)
85#endif
86
87 END SUBROUTINE fftw_dealloc_complex_1d
88! Concatenate the components of the dimensions passed to this function to use it if FFTW3 is not used
89 SUBROUTINE fftw_alloc_complex_2d(array, n)
90 COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:, :), CONTIGUOUS, POINTER, INTENT(OUT) :: array
91 INTEGER, DIMENSION(2), INTENT(IN) :: n
92
93#if defined(__FFTW3)
94 TYPE(c_ptr) :: data_ptr
95 data_ptr = fftw_alloc_complex(int(product(n), kind=c_size_t))
96 CALL c_f_pointer(data_ptr, array, n)
97#else
98! Just allocate the array
99 ALLOCATE (array(n(1), n(2)))
100#endif
101
102 END SUBROUTINE fftw_alloc_complex_2d
103
104 SUBROUTINE fftw_dealloc_complex_2d(array)
105 COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:, :), CONTIGUOUS, POINTER, INTENT(INOUT) :: array
106
107#if defined(__FFTW3)
108 CALL fftw_free(c_loc(array))
109 NULLIFY (array)
110#else
111! Just deallocate the array
112 DEALLOCATE (array)
113#endif
114
115 END SUBROUTINE fftw_dealloc_complex_2d
116! Concatenate the components of the dimensions passed to this function to use it if FFTW3 is not used
117 SUBROUTINE fftw_alloc_complex_3d(array, n)
118 COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:, :, :), CONTIGUOUS, POINTER, INTENT(OUT) :: array
119 INTEGER, DIMENSION(3), INTENT(IN) :: n
120
121#if defined(__FFTW3)
122 TYPE(c_ptr) :: data_ptr
123 data_ptr = fftw_alloc_complex(int(product(n), kind=c_size_t))
124 CALL c_f_pointer(data_ptr, array, n)
125#else
126! Just allocate the array
127 ALLOCATE (array(n(1), n(2), n(3)))
128#endif
129
130 END SUBROUTINE fftw_alloc_complex_3d
131
132 SUBROUTINE fftw_dealloc_complex_3d(array)
133 COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:, :, :), CONTIGUOUS, POINTER, INTENT(INOUT) :: array
134
135#if defined(__FFTW3)
136 CALL fftw_free(c_loc(array))
137 NULLIFY (array)
138#else
139! Just deallocate the array
140 DEALLOCATE (array)
141#endif
142
143 END SUBROUTINE fftw_dealloc_complex_3d
144
145#if defined(__FFTW3)
146! **************************************************************************************************
147!> \brief A workaround that allows us to compile with -Werror=unused-parameter
148! **************************************************************************************************
149 SUBROUTINE dummy_routine_to_call_mark_used()
150 mark_used(fftw_r2hc)
151 mark_used(fftw_hc2r)
152 mark_used(fftw_dht)
153 mark_used(fftw_redft00)
154 mark_used(fftw_redft01)
155 mark_used(fftw_redft10)
156 mark_used(fftw_redft11)
157 mark_used(fftw_rodft00)
158 mark_used(fftw_rodft01)
159 mark_used(fftw_rodft10)
160 mark_used(fftw_rodft11)
161 mark_used(fftw_forward)
162 mark_used(fftw_backward)
163 mark_used(fftw_measure)
164 mark_used(fftw_destroy_input)
165 mark_used(fftw_unaligned)
166 mark_used(fftw_conserve_memory)
167 mark_used(fftw_exhaustive)
168 mark_used(fftw_preserve_input)
169 mark_used(fftw_patient)
170 mark_used(fftw_estimate)
171 mark_used(fftw_wisdom_only)
172 mark_used(fftw_estimate_patient)
173 mark_used(fftw_believe_pcost)
174 mark_used(fftw_no_dft_r2hc)
175 mark_used(fftw_no_nonthreaded)
176 mark_used(fftw_no_buffering)
177 mark_used(fftw_no_indirect_op)
178 mark_used(fftw_allow_large_generic)
179 mark_used(fftw_no_rank_splits)
180 mark_used(fftw_no_vrank_splits)
181 mark_used(fftw_no_vrecurse)
182 mark_used(fftw_no_simd)
183 mark_used(fftw_no_slow)
184 mark_used(fftw_no_fixed_radix_large_n)
185 mark_used(fftw_allow_pruning)
186 END SUBROUTINE dummy_routine_to_call_mark_used
187#endif
188
189! **************************************************************************************************
190!> \brief ...
191!> \param wisdom_file ...
192!> \param ionode ...
193!> \param export_wisdom Export the wisdom to wisdom_file (default: no export). It should only be
194!> requested for a wisdom file specified explicitly in the input, since the
195!> default wisdom file (/etc/fftw/wisdom) is usually not writable.
196! **************************************************************************************************
197 SUBROUTINE fftw3_do_cleanup(wisdom_file, ionode, export_wisdom)
198
199 CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
200 LOGICAL :: ionode
201 LOGICAL, INTENT(IN), OPTIONAL :: export_wisdom
202
203#if defined(__FFTW3)
204 CHARACTER(LEN=1, KIND=C_CHAR), DIMENSION(:), ALLOCATABLE :: wisdom_file_name_c
205 INTEGER :: file_name_length, i
206 INTEGER(KIND=C_INT) :: isuccess
207 LOGICAL :: my_export_wisdom
208
209 my_export_wisdom = .false.
210 IF (PRESENT(export_wisdom)) my_export_wisdom = export_wisdom
211
212 ! Write out FFTW3 wisdom to file, if requested. Only the ionode updates the wisdom.
213 ! FFTW opens the file itself and reports a failure.
214 IF (ionode .AND. my_export_wisdom) THEN
215 file_name_length = len_trim(wisdom_file)
216 ALLOCATE (wisdom_file_name_c(file_name_length + 1))
217 DO i = 1, file_name_length
218 wisdom_file_name_c(i) = wisdom_file(i:i)
219 END DO
220 wisdom_file_name_c(file_name_length + 1) = c_null_char
221 isuccess = fftw_export_wisdom_to_filename(wisdom_file_name_c)
222 IF (isuccess == 0) THEN
223 CALL cp_warn(__location__, "Error exporting wisdom to file "//trim(wisdom_file)//". "// &
224 "Wisdom was not exported.")
225 END IF
226 END IF
227
228 CALL fftw_cleanup()
229#else
230 mark_used(wisdom_file)
231 mark_used(ionode)
232 mark_used(export_wisdom)
233#endif
234
235 END SUBROUTINE fftw3_do_cleanup
236
237! **************************************************************************************************
238!> \brief ...
239!> \param wisdom_file ...
240! **************************************************************************************************
241 SUBROUTINE fftw3_do_init(wisdom_file, plan_style)
242
243 CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
244 INTEGER, INTENT(IN) :: plan_style
245
246#if defined(__FFTW3)
247 CHARACTER(LEN=1, KIND=C_CHAR), DIMENSION(:), ALLOCATABLE :: wisdom_file_name_c
248 INTEGER :: file_name_length, i
249 INTEGER(KIND=C_INT) :: isuccess
250 LOGICAL :: file_exists
251
252 isuccess = fftw_init_threads()
253 IF (isuccess == 0) THEN
254 cpabort("Error initializing FFTW with threads")
255 END IF
256
257 ! Read FFTW wisdom (if available). All nodes read the file here.
258 ! FFTW opens the file itself and reports a failure.
259 INQUIRE (file=wisdom_file, exist=file_exists)
260 IF (file_exists) THEN
261 file_name_length = len_trim(wisdom_file)
262 ALLOCATE (wisdom_file_name_c(file_name_length + 1))
263 DO i = 1, file_name_length
264 wisdom_file_name_c(i) = wisdom_file(i:i)
265 END DO
266 wisdom_file_name_c(file_name_length + 1) = c_null_char
267 isuccess = fftw_import_wisdom_from_filename(wisdom_file_name_c)
268 IF (isuccess == 0) THEN
269 CALL cp_warn(__location__, "Error importing wisdom from file "//trim(wisdom_file)//". "// &
270 "Maybe the file was created with a different configuration than CP2K is "// &
271 "run with. CP2K continues without importing wisdom.")
272 END IF
273 END IF
274
275 SELECT CASE (plan_style)
276 CASE (1)
277 fftw_plan_type = fftw_estimate
278 CASE (2)
279 fftw_plan_type = fftw_measure
280 CASE (3)
281 fftw_plan_type = fftw_patient
282 CASE (4)
283 fftw_plan_type = fftw_exhaustive
284 CASE DEFAULT
285 cpabort("Invalid FFTW_PLAN_TYPE")
286 END SELECT
287#else
288 mark_used(wisdom_file)
289 mark_used(plan_style)
290#endif
291
292 END SUBROUTINE fftw3_do_init
293
294! **************************************************************************************************
295!> \brief ...
296!> \param DATA ...
297!> \param max_length ...
298!> \par History
299!> JGH 23-Jan-2006 : initial version
300!> Adapted for new interface
301!> IAB 09-Jan-2009 : Modified to cache plans in fft_plan_type
302!> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
303!> IAB 09-Oct-2009 : Added OpenMP directives to 1D FFT, and planning routines
304!> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
305!> IAB 11-Sep-2012 : OpenMP parallel 3D FFT (Ruyman Reyes, PRACE)
306!> \author JGH
307! **************************************************************************************************
308 SUBROUTINE fftw3_get_lengths(DATA, max_length)
309
310 INTEGER, DIMENSION(*) :: data
311 INTEGER, INTENT(INOUT) :: max_length
312
313 INTEGER :: h, i, j, k, m, maxn, maxn_elevens, &
314 maxn_fives, maxn_sevens, &
315 maxn_thirteens, maxn_threes, &
316 maxn_twos, ndata, nmax, number
317 INTEGER, ALLOCATABLE, DIMENSION(:) :: dlocal, idx
318
319!------------------------------------------------------------------------------
320! compute ndata
321!! FFTW can do arbitrary(?) lengths, maybe you want to limit them to some
322!! powers of small prime numbers though...
323
324 maxn_twos = 15
325 maxn_threes = 3
326 maxn_fives = 2
327 maxn_sevens = 1
328 maxn_elevens = 1
329 maxn_thirteens = 0
330 maxn = 37748736
331
332 ndata = 0
333 DO h = 0, maxn_twos
334 nmax = huge(0)/2**h
335 DO i = 0, maxn_threes
336 DO j = 0, maxn_fives
337 DO k = 0, maxn_sevens
338 DO m = 0, maxn_elevens
339 number = (3**i)*(5**j)*(7**k)*(11**m)
340
341 IF (number > nmax) cycle
342
343 number = number*2**h
344 IF (number >= maxn) cycle
345
346 ndata = ndata + 1
347 END DO
348 END DO
349 END DO
350 END DO
351 END DO
352
353 ALLOCATE (dlocal(ndata), idx(ndata))
354
355 ndata = 0
356 dlocal(:) = 0
357 DO h = 0, maxn_twos
358 nmax = huge(0)/2**h
359 DO i = 0, maxn_threes
360 DO j = 0, maxn_fives
361 DO k = 0, maxn_sevens
362 DO m = 0, maxn_elevens
363 number = (3**i)*(5**j)*(7**k)*(11**m)
364
365 IF (number > nmax) cycle
366
367 number = number*2**h
368 IF (number >= maxn) cycle
369
370 ndata = ndata + 1
371 dlocal(ndata) = number
372 END DO
373 END DO
374 END DO
375 END DO
376 END DO
377
378 CALL sortint(dlocal, ndata, idx)
379 ndata = min(ndata, max_length)
380 DATA(1:ndata) = dlocal(1:ndata)
381 max_length = ndata
382
383 DEALLOCATE (dlocal, idx)
384
385 END SUBROUTINE fftw3_get_lengths
386
387! **************************************************************************************************
388!> \brief ...
389!> \param iarr ...
390!> \param n ...
391!> \param index ...
392! **************************************************************************************************
393 SUBROUTINE sortint(iarr, n, index)
394
395 INTEGER, INTENT(IN) :: n
396 INTEGER, INTENT(INOUT) :: iarr(1:n)
397 INTEGER, INTENT(OUT) :: index(1:n)
398
399 INTEGER, PARAMETER :: m = 7, nstack = 50
400
401 INTEGER :: a, i, ib, ir, istack(1:nstack), itemp, &
402 j, jstack, k, l, temp
403
404!------------------------------------------------------------------------------
405
406 DO i = 1, n
407 index(i) = i
408 END DO
409 jstack = 0
410 l = 1
411 ir = n
412 DO WHILE (.true.)
413 IF (ir - l < m) THEN
414 DO j = l + 1, ir
415 a = iarr(j)
416 ib = index(j)
417 DO i = j - 1, 0, -1
418 IF (i == 0) EXIT
419 IF (iarr(i) <= a) EXIT
420 iarr(i + 1) = iarr(i)
421 index(i + 1) = index(i)
422 END DO
423 iarr(i + 1) = a
424 index(i + 1) = ib
425 END DO
426 IF (jstack == 0) RETURN
427 ir = istack(jstack)
428 l = istack(jstack - 1)
429 jstack = jstack - 2
430 ELSE
431 k = (l + ir)/2
432 temp = iarr(k)
433 iarr(k) = iarr(l + 1)
434 iarr(l + 1) = temp
435 itemp = index(k)
436 index(k) = index(l + 1)
437 index(l + 1) = itemp
438 IF (iarr(l + 1) > iarr(ir)) THEN
439 temp = iarr(l + 1)
440 iarr(l + 1) = iarr(ir)
441 iarr(ir) = temp
442 itemp = index(l + 1)
443 index(l + 1) = index(ir)
444 index(ir) = itemp
445 END IF
446 IF (iarr(l) > iarr(ir)) THEN
447 temp = iarr(l)
448 iarr(l) = iarr(ir)
449 iarr(ir) = temp
450 itemp = index(l)
451 index(l) = index(ir)
452 index(ir) = itemp
453 END IF
454 IF (iarr(l + 1) > iarr(l)) THEN
455 temp = iarr(l + 1)
456 iarr(l + 1) = iarr(l)
457 iarr(l) = temp
458 itemp = index(l + 1)
459 index(l + 1) = index(l)
460 index(l) = itemp
461 END IF
462 i = l + 1
463 j = ir
464 a = iarr(l)
465 ib = index(l)
466 DO WHILE (.true.)
467 i = i + 1
468 DO WHILE (iarr(i) < a)
469 i = i + 1
470 END DO
471 j = j - 1
472 DO WHILE (iarr(j) > a)
473 j = j - 1
474 END DO
475 IF (j < i) EXIT
476 temp = iarr(i)
477 iarr(i) = iarr(j)
478 iarr(j) = temp
479 itemp = index(i)
480 index(i) = index(j)
481 index(j) = itemp
482 END DO
483 iarr(l) = iarr(j)
484 iarr(j) = a
485 index(l) = index(j)
486 index(j) = ib
487 jstack = jstack + 2
488 IF (jstack > nstack) cpabort("Nstack too small in sortint")
489 IF (ir - i + 1 >= j - l) THEN
490 istack(jstack) = ir
491 istack(jstack - 1) = i
492 ir = j - 1
493 ELSE
494 istack(jstack) = j - 1
495 istack(jstack - 1) = l
496 l = i
497 END IF
498 END IF
499
500 END DO
501
502 END SUBROUTINE sortint
503
504! **************************************************************************************************
505
506! **************************************************************************************************
507!> \brief ...
508!> \param plan ...
509!> \param fft_rank ...
510!> \param dim_n ...
511!> \param dim_istride ...
512!> \param dim_ostride ...
513!> \param hm_rank ...
514!> \param hm_n ...
515!> \param hm_istride ...
516!> \param hm_ostride ...
517!> \param zin ...
518!> \param zout ...
519!> \param fft_direction ...
520!> \param fftw_plan_type ...
521!> \param valid ...
522! **************************************************************************************************
523 SUBROUTINE fftw3_create_guru_plan(plan, fft_rank, dim_n, &
524 dim_istride, dim_ostride, hm_rank, &
525 hm_n, hm_istride, hm_ostride, &
526 zin, zout, fft_direction, fftw_plan_type, &
527 valid)
528
529 TYPE(c_ptr), INTENT(INOUT) :: plan
530 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
531 INTEGER, INTENT(IN) :: dim_n(2), dim_istride(2), dim_ostride(2), &
532 hm_n(2), hm_istride(2), hm_ostride(2), fft_rank, &
533 fft_direction, fftw_plan_type, hm_rank
534 LOGICAL, INTENT(OUT) :: valid
535
536#if defined(__FFTW3)
537 TYPE(fftw_iodim) :: dim(2), hm(2)
538 INTEGER :: i
539
540 DO i = 1, 2
541 dim(i) = fftw_iodim(dim_n(i), dim_istride(i), dim_ostride(i))
542 hm(i) = fftw_iodim(hm_n(i), hm_istride(i), hm_ostride(i))
543 END DO
544
545 plan = fftw_plan_guru_dft(fft_rank, &
546 dim, hm_rank, hm, &
547 zin, zout, &
548 fft_direction, fftw_plan_type)
549
550 valid = c_associated(plan)
551
552#else
553 mark_used(plan)
554 mark_used(fft_rank)
555 mark_used(dim_n)
556 mark_used(dim_istride)
557 mark_used(dim_ostride)
558 mark_used(hm_rank)
559 mark_used(hm_n)
560 mark_used(hm_istride)
561 mark_used(hm_ostride)
562 mark_used(fft_direction)
563 mark_used(fftw_plan_type)
564 !MARK_USED does not work with assumed size arguments
565 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
566 valid = .false.
567
568#endif
569
570 END SUBROUTINE fftw3_create_guru_plan
571
572! **************************************************************************************************
573
574! **************************************************************************************************
575!> \brief Attempt to create a plan with the guru interface for a 2d sub-space.
576!> If this fails, fall back to the FFTW3 threaded 3D transform instead
577!> of the hand-optimised version.
578!> \return ...
579! **************************************************************************************************
580 FUNCTION fftw3_is_guru_supported() RESULT(guru_supported)
581 LOGICAL :: guru_supported
582#if defined(__FFTW3)
583 INTEGER :: dim_n(2), dim_istride(2), dim_ostride(2), &
584 howmany_n(2), howmany_istride(2), howmany_ostride(2)
585 TYPE(c_ptr) :: test_plan
586 COMPLEX(KIND=dp), DIMENSION(1, 1, 1) :: zin
587
588 dim_n(1) = 1
589 dim_n(2) = 1
590 dim_istride(1) = 1
591 dim_istride(2) = 1
592 dim_ostride(1) = 1
593 dim_ostride(2) = 1
594 howmany_n(1) = 1
595 howmany_n(2) = 1
596 howmany_istride(1) = 1
597 howmany_istride(2) = 1
598 howmany_ostride(1) = 1
599 howmany_ostride(2) = 1
600 zin = z_zero
601 CALL fftw3_create_guru_plan(test_plan, 1, &
602 dim_n, dim_istride, dim_ostride, &
603 2, howmany_n, howmany_istride, howmany_ostride, &
604 zin, zin, &
605 fftw_forward, fftw_estimate, guru_supported)
606 IF (guru_supported) THEN
607 CALL fftw_destroy_plan(test_plan)
608 END IF
609
610#else
611 guru_supported = .false.
612#endif
613
614 END FUNCTION fftw3_is_guru_supported
615
616! **************************************************************************************************
617
618! **************************************************************************************************
619!> \brief ...
620!> \param nrows ...
621!> \param nt ...
622!> \param rows_per_thread ...
623!> \param rows_per_thread_r ...
624!> \param th_planA ...
625!> \param th_planB ...
626! **************************************************************************************************
627 SUBROUTINE fftw3_compute_rows_per_th(nrows, nt, rows_per_thread, rows_per_thread_r, &
628 th_planA, th_planB)
629
630 INTEGER, INTENT(IN) :: nrows, nt
631 INTEGER, INTENT(OUT) :: rows_per_thread, rows_per_thread_r, &
632 th_plana, th_planb
633
634 IF (mod(nrows, nt) == 0) THEN
635 rows_per_thread = nrows/nt
636 rows_per_thread_r = 0
637 th_plana = nt
638 th_planb = 0
639 ELSE
640 rows_per_thread = nrows/nt + 1
641 rows_per_thread_r = nrows/nt
642 th_plana = mod(nrows, nt)
643 th_planb = nt - th_plana
644 END IF
645
646 END SUBROUTINE fftw3_compute_rows_per_th
647
648! **************************************************************************************************
649
650! **************************************************************************************************
651!> \brief ...
652!> \param plan ...
653!> \param plan_r ...
654!> \param dim_n ...
655!> \param dim_istride ...
656!> \param dim_ostride ...
657!> \param hm_n ...
658!> \param hm_istride ...
659!> \param hm_ostride ...
660!> \param input ...
661!> \param output ...
662!> \param fft_direction ...
663!> \param fftw_plan_type ...
664!> \param rows_per_th ...
665!> \param rows_per_th_r ...
666! **************************************************************************************************
667 SUBROUTINE fftw3_create_3d_plans(plan, plan_r, dim_n, dim_istride, dim_ostride, &
668 hm_n, hm_istride, hm_ostride, &
669 input, output, &
670 fft_direction, fftw_plan_type, rows_per_th, &
671 rows_per_th_r)
672
673 TYPE(c_ptr), INTENT(INOUT) :: plan, plan_r
674 INTEGER, INTENT(INOUT) :: dim_n(2), dim_istride(2), &
675 dim_ostride(2), hm_n(2), &
676 hm_istride(2), hm_ostride(2)
677 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: input, output
678 INTEGER, INTENT(INOUT) :: fft_direction, fftw_plan_type
679 INTEGER, INTENT(IN) :: rows_per_th, rows_per_th_r
680
681 LOGICAL :: valid
682
683! First plans will have an additional row
684
685 hm_n(2) = rows_per_th
686 CALL fftw3_create_guru_plan(plan, 1, &
687 dim_n, dim_istride, dim_ostride, &
688 2, hm_n, hm_istride, hm_ostride, &
689 input, output, &
690 fft_direction, fftw_plan_type, valid)
691
692 IF (.NOT. valid) THEN
693 cpabort("fftw3_create_plan")
694 END IF
695
696 !!!! Remainder
697 hm_n(2) = rows_per_th_r
698 CALL fftw3_create_guru_plan(plan_r, 1, &
699 dim_n, dim_istride, dim_ostride, &
700 2, hm_n, hm_istride, hm_ostride, &
701 input, output, &
702 fft_direction, fftw_plan_type, valid)
703 IF (.NOT. valid) THEN
704 cpabort("fftw3_create_plan (remaining)")
705 END IF
706
707 END SUBROUTINE fftw3_create_3d_plans
708
709! **************************************************************************************************
710
711! **************************************************************************************************
712!> \brief ...
713!> \param plan ...
714!> \param zin ...
715!> \param zout ...
716! **************************************************************************************************
717 SUBROUTINE fftw3_create_plan_3d(plan, zin, zout)
718
719 TYPE(fft_plan_type), INTENT(INOUT) :: plan
720 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin
721 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zout
722#if defined(__FFTW3)
723 INTEGER :: n1, n2, n3
724 INTEGER :: nt
725 INTEGER :: rows_per_th
726 INTEGER :: rows_per_th_r
727 INTEGER :: fft_direction
728 INTEGER :: th_plana, th_planb
729 COMPLEX(KIND=dp), ALLOCATABLE :: tmp(:)
730
731 ! GURU Interface
732 INTEGER :: dim_n(2), dim_istride(2), dim_ostride(2), &
733 howmany_n(2), howmany_istride(2), howmany_ostride(2)
734
735 IF (plan%fsign == +1) THEN
736 fft_direction = fftw_forward
737 ELSE
738 fft_direction = fftw_backward
739 END IF
740
741 n1 = plan%n_3d(1)
742 n2 = plan%n_3d(2)
743 n3 = plan%n_3d(3)
744
745 nt = 1
746!$OMP PARALLEL DEFAULT(NONE) SHARED(nt)
747!$OMP MASTER
748!$ nt = omp_get_num_threads()
749!$OMP END MASTER
750!$OMP END PARALLEL
751
752 IF ((.NOT. fftw3_is_guru_supported()) .OR. &
753 (.NOT. fftw_plan_type == fftw_estimate) .OR. &
754 (n1 < 256 .AND. n2 < 256 .AND. n3 < 256 .AND. nt == 1)) THEN
755 ! If the plan type is MEASURE, PATIENT and EXHAUSTIVE or
756 ! the grid size is small (and we are single-threaded) then
757 ! FFTW3 does a better job than handmade optimization
758 ! so plan a single 3D FFT which will execute using all the threads
759
760 plan%separated_plans = .false.
761!$ CALL fftw_plan_with_nthreads(nt)
762
763 IF (plan%fft_in_place) THEN
764 plan%fftw_plan = fftw_plan_dft_3d(n3, n2, n1, zin, zin, fft_direction, fftw_plan_type)
765 ELSE
766 plan%fftw_plan = fftw_plan_dft_3d(n3, n2, n1, zin, zout, fft_direction, fftw_plan_type)
767 END IF
768 ELSE
769 ALLOCATE (tmp(n1*n2*n3))
770 ! ************************* PLANS WITH TRANSPOSITIONS ****************************
771 ! In the cases described above, we manually thread each stage of the 3D FFT.
772 !
773 ! The following plans replace the 3D FFT call by running 1D FFTW across all
774 ! 3 directions of the array.
775 !
776 ! Output of FFTW is transposed to ensure that the next round of FFTW access
777 ! contiguous information.
778 !
779 ! Assuming the input matrix is M(n3,n2,n1), FFTW/Transp are :
780 ! M(n3,n2,n1) -> fftw(x) -> M(n3,n1,n2) -> fftw(y) -> M(n1,n2,n3) -> fftw(z) -> M(n1,n2,n3)
781 ! Notice that last matrix is transposed in the Z axis. A DO-loop in the execute routine
782 ! will perform the final transposition. Performance evaluation showed that using an external
783 ! DO loop to do the final transposition performed better than directly transposing the output.
784 ! However, this might vary depending on the compiler/platform, so a potential tuning spot
785 ! is to perform the final transposition within the fftw library rather than using the external loop
786 ! See comments below in Z-FFT for how to transpose the output to avoid the final DO loop.
787 !
788 ! Doc. for the Guru interface is in http://www.fftw.org/doc/Guru-Interface.html
789 !
790 ! OpenMP : Work is distributed on the Z plane.
791 ! All transpositions are out-of-place to facilitate multi-threading
792 !
793 !!!! Plan for X : M(n3,n2,n1) -> fftw(x) -> M(n3,n1,n2)
794 CALL fftw3_compute_rows_per_th(n3, nt, rows_per_th, rows_per_th_r, &
795 th_plana, th_planb)
796
797 dim_n(1) = n1
798 dim_istride(1) = 1
799 dim_ostride(1) = n2
800 howmany_n(1) = n2
801 howmany_n(2) = rows_per_th
802 howmany_istride(1) = n1
803 howmany_istride(2) = n1*n2
804 howmany_ostride(1) = 1
805 howmany_ostride(2) = n1*n2
806 CALL fftw3_create_3d_plans(plan%fftw_plan_nx, plan%fftw_plan_nx_r, &
807 dim_n, dim_istride, dim_ostride, howmany_n, &
808 howmany_istride, howmany_ostride, &
809 zin, tmp, &
810 fft_direction, fftw_plan_type, rows_per_th, &
811 rows_per_th_r)
812
813 !!!! Plan for Y : M(n3,n1,n2) -> fftw(y) -> M(n1,n2,n3)
814 CALL fftw3_compute_rows_per_th(n3, nt, rows_per_th, rows_per_th_r, &
815 th_plana, th_planb)
816 dim_n(1) = n2
817 dim_istride(1) = 1
818 dim_ostride(1) = n3
819 howmany_n(1) = n1
820 howmany_n(2) = rows_per_th
821 howmany_istride(1) = n2
822 howmany_istride(2) = n1*n2
823 !!! transposed Z axis on output
824 howmany_ostride(1) = n2*n3
825 howmany_ostride(2) = 1
826
827 CALL fftw3_create_3d_plans(plan%fftw_plan_ny, plan%fftw_plan_ny_r, &
828 dim_n, dim_istride, dim_ostride, &
829 howmany_n, howmany_istride, howmany_ostride, &
830 tmp, zin, &
831 fft_direction, fftw_plan_type, rows_per_th, &
832 rows_per_th_r)
833
834 !!!! Plan for Z : M(n1,n2,n3) -> fftw(z) -> M(n1,n2,n3)
835 CALL fftw3_compute_rows_per_th(n1, nt, rows_per_th, rows_per_th_r, &
836 th_plana, th_planb)
837 dim_n(1) = n3
838 dim_istride(1) = 1
839 dim_ostride(1) = 1 ! To transpose: n2*n1
840 howmany_n(1) = n2
841 howmany_n(2) = rows_per_th
842 howmany_istride(1) = n3
843 howmany_istride(2) = n2*n3
844 howmany_ostride(1) = n3 ! To transpose: n1
845 howmany_ostride(2) = n2*n3 ! To transpose: 1
846
847 CALL fftw3_create_3d_plans(plan%fftw_plan_nz, plan%fftw_plan_nz_r, &
848 dim_n, dim_istride, dim_ostride, &
849 howmany_n, howmany_istride, howmany_ostride, &
850 zin, tmp, &
851 fft_direction, fftw_plan_type, rows_per_th, &
852 rows_per_th_r)
853
854 plan%separated_plans = .true.
855
856 DEALLOCATE (tmp)
857 END IF
858
859#else
860 mark_used(plan)
861 !MARK_USED does not work with assumed size arguments
862 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
863#endif
864
865 END SUBROUTINE fftw3_create_plan_3d
866
867! **************************************************************************************************
868
869! **************************************************************************************************
870!> \brief ...
871!> \param plan ...
872!> \param plan_r ...
873!> \param split_dim ...
874!> \param nt ...
875!> \param tid ...
876!> \param input ...
877!> \param istride ...
878!> \param output ...
879!> \param ostride ...
880! **************************************************************************************************
881 SUBROUTINE fftw3_workshare_execute_dft(plan, plan_r, split_dim, nt, tid, &
882 input, istride, output, ostride)
883
884 INTEGER, INTENT(IN) :: split_dim, nt, tid
885 INTEGER, INTENT(IN) :: istride, ostride
886 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: input, output
887 TYPE(c_ptr) :: plan, plan_r
888#if defined(__FFTW3)
889 INTEGER :: i_off, o_off
890 INTEGER :: th_plana, th_planb
891 INTEGER :: rows_per_thread, rows_per_thread_r
892
893 CALL fftw3_compute_rows_per_th(split_dim, nt, rows_per_thread, &
894 rows_per_thread_r, &
895 th_plana, th_planb)
896
897 IF (th_planb > 0) THEN
898 IF (tid < th_plana) THEN
899 i_off = (tid)*(istride*(rows_per_thread)) + 1
900 o_off = (tid)*(ostride*(rows_per_thread)) + 1
901 IF (rows_per_thread > 0) THEN
902 CALL fftw_execute_dft(plan, input(i_off), &
903 output(o_off))
904 END IF
905 ELSE IF ((tid - th_plana) < th_planb) THEN
906
907 i_off = (th_plana)*istride*(rows_per_thread) + &
908 (tid - th_plana)*istride*(rows_per_thread_r) + 1
909 o_off = (th_plana)*ostride*(rows_per_thread) + &
910 (tid - th_plana)*ostride*(rows_per_thread_r) + 1
911
912 CALL fftw_execute_dft(plan_r, input(i_off), &
913 output(o_off))
914 END IF
915
916 ELSE
917 i_off = (tid)*(istride*(rows_per_thread)) + 1
918 o_off = (tid)*(ostride*(rows_per_thread)) + 1
919
920 CALL fftw_execute_dft(plan, input(i_off), &
921 output(o_off))
922
923 END IF
924#else
925 mark_used(plan)
926 mark_used(plan_r)
927 mark_used(split_dim)
928 mark_used(nt)
929 mark_used(tid)
930 mark_used(istride)
931 mark_used(ostride)
932 !MARK_USED does not work with assumed size arguments
933 IF (.false.) then; do; IF (abs(input(1)) > abs(output(1))) exit; END do; END IF
934#endif
935
936 END SUBROUTINE fftw3_workshare_execute_dft
937
938! **************************************************************************************************
939
940! **************************************************************************************************
941!> \brief ...
942!> \param plan ...
943!> \param scale ...
944!> \param zin ...
945!> \param zout ...
946!> \param stat ...
947! **************************************************************************************************
948 SUBROUTINE fftw33d(plan, scale, zin, zout, stat)
949
950 TYPE(fft_plan_type), INTENT(IN) :: plan
951 REAL(kind=dp), INTENT(IN) :: scale
952 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT), TARGET:: zin
953 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT), TARGET:: zout
954 INTEGER, INTENT(OUT) :: stat
955#if defined(__FFTW3)
956 COMPLEX(KIND=dp), POINTER :: xout(:)
957 COMPLEX(KIND=dp), ALLOCATABLE :: tmp1(:)
958 INTEGER :: n1, n2, n3
959 INTEGER :: tid, nt
960 INTEGER :: i, j, k
961
962 n1 = plan%n_3d(1)
963 n2 = plan%n_3d(2)
964 n3 = plan%n_3d(3)
965
966 stat = 1
967
968 ! We use a POINTER to the output array to avoid duplicating code
969 IF (plan%fft_in_place) THEN
970 xout => zin(:n1*n2*n3)
971 ELSE
972 xout => zout(:n1*n2*n3)
973 END IF
974
975 ! Either compute the full 3D FFT using a multithreaded plan
976 IF (.NOT. plan%separated_plans) THEN
977 CALL fftw_execute_dft(plan%fftw_plan, zin, xout)
978 ELSE
979 ! Or use the 3 stage FFT scheme described in fftw3_create_plan_3d
980 ALLOCATE (tmp1(n1*n2*n3)) ! Temporary vector used for transpositions
981 !$OMP PARALLEL DEFAULT(NONE) PRIVATE(tid,nt,i,j,k) SHARED(zin,tmp1,n1,n2,n3,plan,xout)
982 tid = 0
983 nt = 1
984
985!$ tid = omp_get_thread_num()
986!$ nt = omp_get_num_threads()
987 CALL fftw3_workshare_execute_dft(plan%fftw_plan_nx, plan%fftw_plan_nx_r, &
988 n3, nt, tid, &
989 zin, n1*n2, tmp1, n1*n2)
990
991 !$OMP BARRIER
992 CALL fftw3_workshare_execute_dft(plan%fftw_plan_ny, plan%fftw_plan_ny_r, &
993 n3, nt, tid, &
994 tmp1, n1*n2, xout, 1)
995 !$OMP BARRIER
996 CALL fftw3_workshare_execute_dft(plan%fftw_plan_nz, plan%fftw_plan_nz_r, &
997 n1, nt, tid, &
998 xout, n2*n3, tmp1, n2*n3)
999 !$OMP BARRIER
1000
1001 !$OMP DO COLLAPSE(3)
1002 DO i = 1, n1
1003 DO j = 1, n2
1004 DO k = 1, n3
1005 xout((i - 1) + (j - 1)*n1 + (k - 1)*n1*n2 + 1) = &
1006 tmp1((k - 1) + (j - 1)*n3 + (i - 1)*n3*n2 + 1)
1007 END DO
1008 END DO
1009 END DO
1010 !$OMP END DO
1011
1012 !$OMP END PARALLEL
1013 END IF
1014
1015 IF (scale /= 1.0_dp) THEN
1016 CALL zdscal(n1*n2*n3, scale, xout, 1)
1017 END IF
1018
1019#else
1020 mark_used(plan)
1021 mark_used(scale)
1022 !MARK_USED does not work with assumed size arguments
1023 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
1024 stat = 0
1025
1026#endif
1027
1028 END SUBROUTINE fftw33d
1029
1030! **************************************************************************************************
1031
1032! **************************************************************************************************
1033!> \brief ...
1034!> \param plan ...
1035!> \param zin ...
1036!> \param zout ...
1037! **************************************************************************************************
1038 SUBROUTINE fftw3_create_plan_1d(plan, zin, zout)
1039 TYPE(fft_plan_type), INTENT(INOUT) :: plan
1040 COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: zin
1041 COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: zout
1042#if defined(__FFTW3)
1043 INTEGER :: istride, idist, ostride, odist, num_threads, num_rows
1044
1045 num_threads = 1
1046 plan%separated_plans = .false.
1047!$OMP PARALLEL DEFAULT(NONE), &
1048!$OMP SHARED(NUM_THREADS)
1049!$OMP MASTER
1050!$ num_threads = omp_get_num_threads()
1051!$OMP END MASTER
1052!$OMP END PARALLEL
1053
1054 num_rows = plan%m/num_threads
1055!$ plan%num_threads_needed = num_threads
1056
1057! Check for number of rows less than num_threads
1058!$ IF (plan%m < num_threads) THEN
1059!$ num_rows = 1
1060!$ plan%num_threads_needed = plan%m
1061!$ END IF
1062
1063! Check for total number of rows not divisible by num_threads
1064!$ IF (num_rows*plan%num_threads_needed /= plan%m) THEN
1065!$ plan%need_alt_plan = .TRUE.
1066!$ END IF
1067
1068!$ plan%num_rows = num_rows
1069 istride = 1
1070 idist = plan%ldx_in
1071 ostride = 1
1072 odist = plan%ldx_out
1073 IF (plan%trans_in) THEN
1074 istride = plan%ldx_in
1075 idist = 1
1076 END IF
1077 IF (plan%trans_out) THEN
1078 ostride = plan%ldx_out
1079 odist = 1
1080 END IF
1081
1082 IF (plan%fsign == +1) THEN
1083 CALL dfftw_plan_many_dft(plan%fftw_plan, 1, plan%n, num_rows, zin, 0, istride, idist, &
1084 zout, 0, ostride, odist, fftw_forward, fftw_plan_type)
1085 ELSE
1086 CALL dfftw_plan_many_dft(plan%fftw_plan, 1, plan%n, num_rows, zin, 0, istride, idist, &
1087 zout, 0, ostride, odist, fftw_backward, fftw_plan_type)
1088 END IF
1089
1090!$ IF (plan%need_alt_plan) THEN
1091!$ plan%alt_num_rows = plan%m - (plan%num_threads_needed - 1)*num_rows
1092!$ IF (plan%fsign == +1) THEN
1093!$ CALL dfftw_plan_many_dft(plan%alt_fftw_plan, 1, plan%n, plan%alt_num_rows, zin, 0, istride, idist, &
1094!$ zout, 0, ostride, odist, FFTW_FORWARD, fftw_plan_type)
1095!$ ELSE
1096!$ CALL dfftw_plan_many_dft(plan%alt_fftw_plan, 1, plan%n, plan%alt_num_rows, zin, 0, istride, idist, &
1097!$ zout, 0, ostride, odist, FFTW_BACKWARD, fftw_plan_type)
1098!$ END IF
1099!$ END IF
1100
1101#else
1102 mark_used(plan)
1103 !MARK_USED does not work with assumed size arguments
1104 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
1105#endif
1106
1107 END SUBROUTINE fftw3_create_plan_1d
1108
1109! **************************************************************************************************
1110!> \brief ...
1111!> \param plan ...
1112! **************************************************************************************************
1113 SUBROUTINE fftw3_destroy_plan(plan)
1114
1115 TYPE(fft_plan_type), INTENT(INOUT) :: plan
1116
1117#if defined(__FFTW3)
1118!$ IF (plan%need_alt_plan) THEN
1119!$ CALL fftw_destroy_plan(plan%alt_fftw_plan)
1120!$ END IF
1121
1122 IF (.NOT. plan%separated_plans) THEN
1123 CALL fftw_destroy_plan(plan%fftw_plan)
1124 ELSE
1125 ! If it is a separated plan then we have to destroy
1126 ! each dim plan individually
1127 CALL fftw_destroy_plan(plan%fftw_plan_nx)
1128 CALL fftw_destroy_plan(plan%fftw_plan_ny)
1129 CALL fftw_destroy_plan(plan%fftw_plan_nz)
1130 CALL fftw_destroy_plan(plan%fftw_plan_nx_r)
1131 CALL fftw_destroy_plan(plan%fftw_plan_ny_r)
1132 CALL fftw_destroy_plan(plan%fftw_plan_nz_r)
1133 END IF
1134
1135#else
1136 mark_used(plan)
1137#endif
1138
1139 END SUBROUTINE fftw3_destroy_plan
1140
1141! **************************************************************************************************
1142!> \brief ...
1143!> \param plan ...
1144!> \param zin ...
1145!> \param zout ...
1146!> \param scale ...
1147!> \param stat ...
1148! **************************************************************************************************
1149 SUBROUTINE fftw31dm(plan, zin, zout, scale, stat)
1150 TYPE(fft_plan_type), INTENT(IN) :: plan
1151 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT), &
1152 TARGET :: zin
1153 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT), &
1154 TARGET :: zout
1155 REAL(kind=dp), INTENT(IN) :: scale
1156 INTEGER, INTENT(OUT) :: stat
1157
1158 INTEGER :: in_offset, my_id, num_rows, out_offset, &
1159 scal_offset
1160 TYPE(c_ptr) :: fftw_plan
1161!------------------------------------------------------------------------------
1162
1163 my_id = 0
1164 num_rows = plan%m
1165
1166#if defined(__FFTW3)
1167 IF (plan%m <= 1) THEN
1168 stat = 1
1169 CALL fftw_execute_dft(plan%fftw_plan, zin(1), zout(1))
1170 IF (scale /= 1.0_dp) CALL zdscal(plan%n*plan%m, scale, zout, 1)
1171 RETURN
1172 END IF
1173#endif
1174
1175!$OMP PARALLEL DEFAULT(NONE), &
1176!$OMP PRIVATE(my_id,num_rows,in_offset,out_offset,scal_offset,fftw_plan), &
1177!$OMP SHARED(zin,zout), &
1178!$OMP SHARED(plan,scale,stat)
1179!$ my_id = omp_get_thread_num()
1180
1181!$ if (my_id < plan%num_threads_needed) then
1182
1183 fftw_plan = plan%fftw_plan
1184
1185 in_offset = 1
1186 out_offset = 1
1187 scal_offset = 1
1188
1189!$ in_offset = 1 + plan%num_rows*my_id*plan%n
1190!$ out_offset = 1 + plan%num_rows*my_id*plan%n
1191!$ IF (plan%trans_in) THEN
1192!$ in_offset = 1 + plan%num_rows*my_id
1193!$ END IF
1194!$ IF (plan%trans_out) THEN
1195!$ out_offset = 1 + plan%num_rows*my_id
1196!$ END IF
1197!$ scal_offset = 1 + plan%n*plan%num_rows*my_id
1198!$ IF (plan%need_alt_plan .AND. my_id == plan%num_threads_needed - 1) THEN
1199!$ num_rows = plan%alt_num_rows
1200!$ fftw_plan = plan%alt_fftw_plan
1201!$ ELSE
1202!$ num_rows = plan%num_rows
1203!$ END IF
1204
1205#if defined(__FFTW3)
1206!$OMP MASTER
1207 stat = 1
1208!$OMP END MASTER
1209 CALL fftw_execute_dft(fftw_plan, zin(in_offset), zout(out_offset))
1210!$ end if
1211! all threads need to meet at this barrier
1212!$OMP BARRIER
1213!$ if (my_id < plan%num_threads_needed) then
1214 IF (scale /= 1.0_dp) CALL zdscal(plan%n*num_rows, scale, zout(scal_offset:scal_offset), 1)
1215!$ end if
1216
1217#else
1218 mark_used(plan)
1219 mark_used(scale)
1220 !MARK_USED does not work with assumed size arguments
1221 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
1222 stat = 0
1223
1224!$ else
1225!$ end if
1226
1227#endif
1228
1229!$OMP END PARALLEL
1230
1231 END SUBROUTINE fftw31dm
1232
1233 END MODULE fftw3_lib
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Defines the basic variable types.
Definition fft_kinds.F:13
integer, parameter, public dp
Definition fft_kinds.F:18
Type to store data about a (1D or 3D) FFT, including FFTW plan.
Definition fft_plan.F:18
subroutine, public fftw3_do_cleanup(wisdom_file, ionode, export_wisdom)
...
Definition fftw3_lib.F:198
subroutine, public fftw3_destroy_plan(plan)
...
Definition fftw3_lib.F:1114
subroutine, public fftw3_get_lengths(data, max_length)
...
Definition fftw3_lib.F:309
subroutine, public fftw3_create_plan_1d(plan, zin, zout)
...
Definition fftw3_lib.F:1039
subroutine, public fftw3_create_plan_3d(plan, zin, zout)
...
Definition fftw3_lib.F:718
subroutine, public fftw33d(plan, scale, zin, zout, stat)
...
Definition fftw3_lib.F:949
subroutine, public fftw3_do_init(wisdom_file, plan_style)
...
Definition fftw3_lib.F:242
subroutine, public fftw31dm(plan, zin, zout, scale, stat)
...
Definition fftw3_lib.F:1150
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public z_zero