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