(git:30099c3)
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 IMPLICIT NONE
515
516 TYPE(c_ptr), INTENT(INOUT) :: plan
517 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
518 INTEGER, INTENT(IN) :: dim_n(2), dim_istride(2), dim_ostride(2), &
519 hm_n(2), hm_istride(2), hm_ostride(2), fft_rank, &
520 fft_direction, fftw_plan_type, hm_rank
521 LOGICAL, INTENT(OUT) :: valid
522
523#if defined(__FFTW3)
524 TYPE(fftw_iodim) :: dim(2), hm(2)
525 INTEGER :: i
526
527 DO i = 1, 2
528 dim(i) = fftw_iodim(dim_n(i), dim_istride(i), dim_ostride(i))
529 hm(i) = fftw_iodim(hm_n(i), hm_istride(i), hm_ostride(i))
530 END DO
531
532 plan = fftw_plan_guru_dft(fft_rank, &
533 dim, hm_rank, hm, &
534 zin, zout, &
535 fft_direction, fftw_plan_type)
536
537 valid = c_associated(plan)
538
539#else
540 mark_used(plan)
541 mark_used(fft_rank)
542 mark_used(dim_n)
543 mark_used(dim_istride)
544 mark_used(dim_ostride)
545 mark_used(hm_rank)
546 mark_used(hm_n)
547 mark_used(hm_istride)
548 mark_used(hm_ostride)
549 mark_used(fft_direction)
550 mark_used(fftw_plan_type)
551 !MARK_USED does not work with assumed size arguments
552 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
553 valid = .false.
554
555#endif
556
557 END SUBROUTINE fftw3_create_guru_plan
558
559! **************************************************************************************************
560
561! **************************************************************************************************
562!> \brief Attempt to create a plan with the guru interface for a 2d sub-space.
563!> If this fails, fall back to the FFTW3 threaded 3D transform instead
564!> of the hand-optimised version.
565!> \return ...
566! **************************************************************************************************
567 FUNCTION fftw3_is_guru_supported() RESULT(guru_supported)
568
569 IMPLICIT NONE
570
571 LOGICAL :: guru_supported
572#if defined(__FFTW3)
573 INTEGER :: dim_n(2), dim_istride(2), dim_ostride(2), &
574 howmany_n(2), howmany_istride(2), howmany_ostride(2)
575 TYPE(c_ptr) :: test_plan
576 COMPLEX(KIND=dp), DIMENSION(1, 1, 1) :: zin
577
578 dim_n(1) = 1
579 dim_n(2) = 1
580 dim_istride(1) = 1
581 dim_istride(2) = 1
582 dim_ostride(1) = 1
583 dim_ostride(2) = 1
584 howmany_n(1) = 1
585 howmany_n(2) = 1
586 howmany_istride(1) = 1
587 howmany_istride(2) = 1
588 howmany_ostride(1) = 1
589 howmany_ostride(2) = 1
590 zin = z_zero
591 CALL fftw3_create_guru_plan(test_plan, 1, &
592 dim_n, dim_istride, dim_ostride, &
593 2, howmany_n, howmany_istride, howmany_ostride, &
594 zin, zin, &
595 fftw_forward, fftw_estimate, guru_supported)
596 IF (guru_supported) THEN
597 CALL fftw_destroy_plan(test_plan)
598 END IF
599
600#else
601 guru_supported = .false.
602#endif
603
604 END FUNCTION fftw3_is_guru_supported
605
606! **************************************************************************************************
607
608! **************************************************************************************************
609!> \brief ...
610!> \param nrows ...
611!> \param nt ...
612!> \param rows_per_thread ...
613!> \param rows_per_thread_r ...
614!> \param th_planA ...
615!> \param th_planB ...
616! **************************************************************************************************
617 SUBROUTINE fftw3_compute_rows_per_th(nrows, nt, rows_per_thread, rows_per_thread_r, &
618 th_planA, th_planB)
619
620 INTEGER, INTENT(IN) :: nrows, nt
621 INTEGER, INTENT(OUT) :: rows_per_thread, rows_per_thread_r, &
622 th_plana, th_planb
623
624 IF (mod(nrows, nt) == 0) THEN
625 rows_per_thread = nrows/nt
626 rows_per_thread_r = 0
627 th_plana = nt
628 th_planb = 0
629 ELSE
630 rows_per_thread = nrows/nt + 1
631 rows_per_thread_r = nrows/nt
632 th_plana = mod(nrows, nt)
633 th_planb = nt - th_plana
634 END IF
635
636 END SUBROUTINE fftw3_compute_rows_per_th
637
638! **************************************************************************************************
639
640! **************************************************************************************************
641!> \brief ...
642!> \param plan ...
643!> \param plan_r ...
644!> \param dim_n ...
645!> \param dim_istride ...
646!> \param dim_ostride ...
647!> \param hm_n ...
648!> \param hm_istride ...
649!> \param hm_ostride ...
650!> \param input ...
651!> \param output ...
652!> \param fft_direction ...
653!> \param fftw_plan_type ...
654!> \param rows_per_th ...
655!> \param rows_per_th_r ...
656! **************************************************************************************************
657 SUBROUTINE fftw3_create_3d_plans(plan, plan_r, dim_n, dim_istride, dim_ostride, &
658 hm_n, hm_istride, hm_ostride, &
659 input, output, &
660 fft_direction, fftw_plan_type, rows_per_th, &
661 rows_per_th_r)
662
663 TYPE(c_ptr), INTENT(INOUT) :: plan, plan_r
664 INTEGER, INTENT(INOUT) :: dim_n(2), dim_istride(2), &
665 dim_ostride(2), hm_n(2), &
666 hm_istride(2), hm_ostride(2)
667 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: input, output
668 INTEGER, INTENT(INOUT) :: fft_direction, fftw_plan_type
669 INTEGER, INTENT(IN) :: rows_per_th, rows_per_th_r
670
671 LOGICAL :: valid
672
673! First plans will have an additional row
674
675 hm_n(2) = rows_per_th
676 CALL fftw3_create_guru_plan(plan, 1, &
677 dim_n, dim_istride, dim_ostride, &
678 2, hm_n, hm_istride, hm_ostride, &
679 input, output, &
680 fft_direction, fftw_plan_type, valid)
681
682 IF (.NOT. valid) THEN
683 cpabort("fftw3_create_plan")
684 END IF
685
686 !!!! Remainder
687 hm_n(2) = rows_per_th_r
688 CALL fftw3_create_guru_plan(plan_r, 1, &
689 dim_n, dim_istride, dim_ostride, &
690 2, hm_n, hm_istride, hm_ostride, &
691 input, output, &
692 fft_direction, fftw_plan_type, valid)
693 IF (.NOT. valid) THEN
694 cpabort("fftw3_create_plan (remaining)")
695 END IF
696
697 END SUBROUTINE fftw3_create_3d_plans
698
699! **************************************************************************************************
700
701! **************************************************************************************************
702!> \brief ...
703!> \param plan ...
704!> \param zin ...
705!> \param zout ...
706!> \param plan_style ...
707! **************************************************************************************************
708 SUBROUTINE fftw3_create_plan_3d(plan, zin, zout, plan_style)
709
710 TYPE(fft_plan_type), INTENT(INOUT) :: plan
711 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin
712 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zout
713 INTEGER :: plan_style
714#if defined(__FFTW3)
715 INTEGER :: n1, n2, n3
716 INTEGER :: nt
717 INTEGER :: rows_per_th
718 INTEGER :: rows_per_th_r
719 INTEGER :: fft_direction
720 INTEGER :: th_plana, th_planb
721 COMPLEX(KIND=dp), ALLOCATABLE :: tmp(:)
722
723 ! GURU Interface
724 INTEGER :: dim_n(2), dim_istride(2), dim_ostride(2), &
725 howmany_n(2), howmany_istride(2), howmany_ostride(2)
726
727 INTEGER :: fftw_plan_type
728 SELECT CASE (plan_style)
729 CASE (1)
730 fftw_plan_type = fftw_estimate
731 CASE (2)
732 fftw_plan_type = fftw_measure
733 CASE (3)
734 fftw_plan_type = fftw_patient
735 CASE (4)
736 fftw_plan_type = fftw_exhaustive
737 CASE DEFAULT
738 cpabort("fftw3_create_plan_3d")
739 END SELECT
740
741 IF (plan%fsign == +1) THEN
742 fft_direction = fftw_forward
743 ELSE
744 fft_direction = fftw_backward
745 END IF
746
747 n1 = plan%n_3d(1)
748 n2 = plan%n_3d(2)
749 n3 = plan%n_3d(3)
750
751 nt = 1
752!$OMP PARALLEL DEFAULT(NONE) SHARED(nt)
753!$OMP MASTER
754!$ nt = omp_get_num_threads()
755!$OMP END MASTER
756!$OMP END PARALLEL
757
758 IF ((.NOT. fftw3_is_guru_supported()) .OR. &
759 (.NOT. plan_style == 1) .OR. &
760 (n1 < 256 .AND. n2 < 256 .AND. n3 < 256 .AND. nt == 1)) THEN
761 ! If the plan type is MEASURE, PATIENT and EXHAUSTIVE or
762 ! the grid size is small (and we are single-threaded) then
763 ! FFTW3 does a better job than handmade optimization
764 ! so plan a single 3D FFT which will execute using all the threads
765
766 plan%separated_plans = .false.
767!$ CALL fftw_plan_with_nthreads(nt)
768
769 IF (plan%fft_in_place) THEN
770 plan%fftw_plan = fftw_plan_dft_3d(n3, n2, n1, zin, zin, fft_direction, fftw_plan_type)
771 ELSE
772 plan%fftw_plan = fftw_plan_dft_3d(n3, n2, n1, zin, zout, fft_direction, fftw_plan_type)
773 END IF
774 ELSE
775 ALLOCATE (tmp(n1*n2*n3))
776 ! ************************* PLANS WITH TRANSPOSITIONS ****************************
777 ! In the cases described above, we manually thread each stage of the 3D FFT.
778 !
779 ! The following plans replace the 3D FFT call by running 1D FFTW across all
780 ! 3 directions of the array.
781 !
782 ! Output of FFTW is transposed to ensure that the next round of FFTW access
783 ! contiguous information.
784 !
785 ! Assuming the input matrix is M(n3,n2,n1), FFTW/Transp are :
786 ! M(n3,n2,n1) -> fftw(x) -> M(n3,n1,n2) -> fftw(y) -> M(n1,n2,n3) -> fftw(z) -> M(n1,n2,n3)
787 ! Notice that last matrix is transposed in the Z axis. A DO-loop in the execute routine
788 ! will perform the final transposition. Performance evaluation showed that using an external
789 ! DO loop to do the final transposition performed better than directly transposing the output.
790 ! However, this might vary depending on the compiler/platform, so a potential tuning spot
791 ! is to perform the final transposition within the fftw library rather than using the external loop
792 ! See comments below in Z-FFT for how to transpose the output to avoid the final DO loop.
793 !
794 ! Doc. for the Guru interface is in http://www.fftw.org/doc/Guru-Interface.html
795 !
796 ! OpenMP : Work is distributed on the Z plane.
797 ! All transpositions are out-of-place to facilitate multi-threading
798 !
799 !!!! Plan for X : M(n3,n2,n1) -> fftw(x) -> M(n3,n1,n2)
800 CALL fftw3_compute_rows_per_th(n3, nt, rows_per_th, rows_per_th_r, &
801 th_plana, th_planb)
802
803 dim_n(1) = n1
804 dim_istride(1) = 1
805 dim_ostride(1) = n2
806 howmany_n(1) = n2
807 howmany_n(2) = rows_per_th
808 howmany_istride(1) = n1
809 howmany_istride(2) = n1*n2
810 howmany_ostride(1) = 1
811 howmany_ostride(2) = n1*n2
812 CALL fftw3_create_3d_plans(plan%fftw_plan_nx, plan%fftw_plan_nx_r, &
813 dim_n, dim_istride, dim_ostride, howmany_n, &
814 howmany_istride, howmany_ostride, &
815 zin, tmp, &
816 fft_direction, fftw_plan_type, rows_per_th, &
817 rows_per_th_r)
818
819 !!!! Plan for Y : M(n3,n1,n2) -> fftw(y) -> M(n1,n2,n3)
820 CALL fftw3_compute_rows_per_th(n3, nt, rows_per_th, rows_per_th_r, &
821 th_plana, th_planb)
822 dim_n(1) = n2
823 dim_istride(1) = 1
824 dim_ostride(1) = n3
825 howmany_n(1) = n1
826 howmany_n(2) = rows_per_th
827 howmany_istride(1) = n2
828 howmany_istride(2) = n1*n2
829 !!! transposed Z axis on output
830 howmany_ostride(1) = n2*n3
831 howmany_ostride(2) = 1
832
833 CALL fftw3_create_3d_plans(plan%fftw_plan_ny, plan%fftw_plan_ny_r, &
834 dim_n, dim_istride, dim_ostride, &
835 howmany_n, howmany_istride, howmany_ostride, &
836 tmp, zin, &
837 fft_direction, fftw_plan_type, rows_per_th, &
838 rows_per_th_r)
839
840 !!!! Plan for Z : M(n1,n2,n3) -> fftw(z) -> M(n1,n2,n3)
841 CALL fftw3_compute_rows_per_th(n1, nt, rows_per_th, rows_per_th_r, &
842 th_plana, th_planb)
843 dim_n(1) = n3
844 dim_istride(1) = 1
845 dim_ostride(1) = 1 ! To transpose: n2*n1
846 howmany_n(1) = n2
847 howmany_n(2) = rows_per_th
848 howmany_istride(1) = n3
849 howmany_istride(2) = n2*n3
850 howmany_ostride(1) = n3 ! To transpose: n1
851 howmany_ostride(2) = n2*n3 ! To transpose: 1
852
853 CALL fftw3_create_3d_plans(plan%fftw_plan_nz, plan%fftw_plan_nz_r, &
854 dim_n, dim_istride, dim_ostride, &
855 howmany_n, howmany_istride, howmany_ostride, &
856 zin, tmp, &
857 fft_direction, fftw_plan_type, rows_per_th, &
858 rows_per_th_r)
859
860 plan%separated_plans = .true.
861
862 DEALLOCATE (tmp)
863 END IF
864
865#else
866 mark_used(plan)
867 mark_used(plan_style)
868 !MARK_USED does not work with assumed size arguments
869 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
870#endif
871
872 END SUBROUTINE fftw3_create_plan_3d
873
874! **************************************************************************************************
875
876! **************************************************************************************************
877!> \brief ...
878!> \param plan ...
879!> \param plan_r ...
880!> \param split_dim ...
881!> \param nt ...
882!> \param tid ...
883!> \param input ...
884!> \param istride ...
885!> \param output ...
886!> \param ostride ...
887! **************************************************************************************************
888 SUBROUTINE fftw3_workshare_execute_dft(plan, plan_r, split_dim, nt, tid, &
889 input, istride, output, ostride)
890
891 INTEGER, INTENT(IN) :: split_dim, nt, tid
892 INTEGER, INTENT(IN) :: istride, ostride
893 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: input, output
894 TYPE(c_ptr) :: plan, plan_r
895#if defined(__FFTW3)
896 INTEGER :: i_off, o_off
897 INTEGER :: th_plana, th_planb
898 INTEGER :: rows_per_thread, rows_per_thread_r
899
900 CALL fftw3_compute_rows_per_th(split_dim, nt, rows_per_thread, &
901 rows_per_thread_r, &
902 th_plana, th_planb)
903
904 IF (th_planb > 0) THEN
905 IF (tid < th_plana) THEN
906 i_off = (tid)*(istride*(rows_per_thread)) + 1
907 o_off = (tid)*(ostride*(rows_per_thread)) + 1
908 IF (rows_per_thread > 0) THEN
909 CALL fftw_execute_dft(plan, input(i_off), &
910 output(o_off))
911 END IF
912 ELSE IF ((tid - th_plana) < th_planb) THEN
913
914 i_off = (th_plana)*istride*(rows_per_thread) + &
915 (tid - th_plana)*istride*(rows_per_thread_r) + 1
916 o_off = (th_plana)*ostride*(rows_per_thread) + &
917 (tid - th_plana)*ostride*(rows_per_thread_r) + 1
918
919 CALL fftw_execute_dft(plan_r, input(i_off), &
920 output(o_off))
921 END IF
922
923 ELSE
924 i_off = (tid)*(istride*(rows_per_thread)) + 1
925 o_off = (tid)*(ostride*(rows_per_thread)) + 1
926
927 CALL fftw_execute_dft(plan, input(i_off), &
928 output(o_off))
929
930 END IF
931#else
932 mark_used(plan)
933 mark_used(plan_r)
934 mark_used(split_dim)
935 mark_used(nt)
936 mark_used(tid)
937 mark_used(istride)
938 mark_used(ostride)
939 !MARK_USED does not work with assumed size arguments
940 IF (.false.) then; do; IF (abs(input(1)) > abs(output(1))) exit; END do; END IF
941#endif
942
943 END SUBROUTINE fftw3_workshare_execute_dft
944
945! **************************************************************************************************
946
947! **************************************************************************************************
948!> \brief ...
949!> \param plan ...
950!> \param scale ...
951!> \param zin ...
952!> \param zout ...
953!> \param stat ...
954! **************************************************************************************************
955 SUBROUTINE fftw33d(plan, scale, zin, zout, stat)
956
957 TYPE(fft_plan_type), INTENT(IN) :: plan
958 REAL(kind=dp), INTENT(IN) :: scale
959 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT), TARGET:: zin
960 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT), TARGET:: zout
961 INTEGER, INTENT(OUT) :: stat
962#if defined(__FFTW3)
963 COMPLEX(KIND=dp), POINTER :: xout(:)
964 COMPLEX(KIND=dp), ALLOCATABLE :: tmp1(:)
965 INTEGER :: n1, n2, n3
966 INTEGER :: tid, nt
967 INTEGER :: i, j, k
968
969 n1 = plan%n_3d(1)
970 n2 = plan%n_3d(2)
971 n3 = plan%n_3d(3)
972
973 stat = 1
974
975 ! We use a POINTER to the output array to avoid duplicating code
976 IF (plan%fft_in_place) THEN
977 xout => zin(:n1*n2*n3)
978 ELSE
979 xout => zout(:n1*n2*n3)
980 END IF
981
982 ! Either compute the full 3D FFT using a multithreaded plan
983 IF (.NOT. plan%separated_plans) THEN
984 CALL fftw_execute_dft(plan%fftw_plan, zin, xout)
985 ELSE
986 ! Or use the 3 stage FFT scheme described in fftw3_create_plan_3d
987 ALLOCATE (tmp1(n1*n2*n3)) ! Temporary vector used for transpositions
988 !$OMP PARALLEL DEFAULT(NONE) PRIVATE(tid,nt,i,j,k) SHARED(zin,tmp1,n1,n2,n3,plan,xout)
989 tid = 0
990 nt = 1
991
992!$ tid = omp_get_thread_num()
993!$ nt = omp_get_num_threads()
994 CALL fftw3_workshare_execute_dft(plan%fftw_plan_nx, plan%fftw_plan_nx_r, &
995 n3, nt, tid, &
996 zin, n1*n2, tmp1, n1*n2)
997
998 !$OMP BARRIER
999 CALL fftw3_workshare_execute_dft(plan%fftw_plan_ny, plan%fftw_plan_ny_r, &
1000 n3, nt, tid, &
1001 tmp1, n1*n2, xout, 1)
1002 !$OMP BARRIER
1003 CALL fftw3_workshare_execute_dft(plan%fftw_plan_nz, plan%fftw_plan_nz_r, &
1004 n1, nt, tid, &
1005 xout, n2*n3, tmp1, n2*n3)
1006 !$OMP BARRIER
1007
1008 !$OMP DO COLLAPSE(3)
1009 DO i = 1, n1
1010 DO j = 1, n2
1011 DO k = 1, n3
1012 xout((i - 1) + (j - 1)*n1 + (k - 1)*n1*n2 + 1) = &
1013 tmp1((k - 1) + (j - 1)*n3 + (i - 1)*n3*n2 + 1)
1014 END DO
1015 END DO
1016 END DO
1017 !$OMP END DO
1018
1019 !$OMP END PARALLEL
1020 END IF
1021
1022 IF (scale /= 1.0_dp) THEN
1023 CALL zdscal(n1*n2*n3, scale, xout, 1)
1024 END IF
1025
1026#else
1027 mark_used(plan)
1028 mark_used(scale)
1029 !MARK_USED does not work with assumed size arguments
1030 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
1031 stat = 0
1032
1033#endif
1034
1035 END SUBROUTINE fftw33d
1036
1037! **************************************************************************************************
1038
1039! **************************************************************************************************
1040!> \brief ...
1041!> \param plan ...
1042!> \param zin ...
1043!> \param zout ...
1044!> \param plan_style ...
1045! **************************************************************************************************
1046 SUBROUTINE fftw3_create_plan_1dm(plan, zin, zout, plan_style)
1047
1048 IMPLICIT NONE
1049
1050 TYPE(fft_plan_type), INTENT(INOUT) :: plan
1051 COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: zin
1052 COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: zout
1053 INTEGER, INTENT(IN) :: plan_style
1054#if defined(__FFTW3)
1055 INTEGER :: istride, idist, ostride, odist, num_threads, num_rows
1056
1057 INTEGER :: fftw_plan_type
1058 SELECT CASE (plan_style)
1059 CASE (1)
1060 fftw_plan_type = fftw_estimate
1061 CASE (2)
1062 fftw_plan_type = fftw_measure
1063 CASE (3)
1064 fftw_plan_type = fftw_patient
1065 CASE (4)
1066 fftw_plan_type = fftw_exhaustive
1067 CASE DEFAULT
1068 cpabort("fftw3_create_plan_1dm")
1069 END SELECT
1070
1071 num_threads = 1
1072 plan%separated_plans = .false.
1073!$OMP PARALLEL DEFAULT(NONE), &
1074!$OMP SHARED(NUM_THREADS)
1075!$OMP MASTER
1076!$ num_threads = omp_get_num_threads()
1077!$OMP END MASTER
1078!$OMP END PARALLEL
1079
1080 num_rows = plan%m/num_threads
1081!$ plan%num_threads_needed = num_threads
1082
1083! Check for number of rows less than num_threads
1084!$ IF (plan%m < num_threads) THEN
1085!$ num_rows = 1
1086!$ plan%num_threads_needed = plan%m
1087!$ END IF
1088
1089! Check for total number of rows not divisible by num_threads
1090!$ IF (num_rows*plan%num_threads_needed /= plan%m) THEN
1091!$ plan%need_alt_plan = .TRUE.
1092!$ END IF
1093
1094!$ plan%num_rows = num_rows
1095 istride = 1
1096 idist = plan%n
1097 ostride = 1
1098 odist = plan%n
1099 IF (plan%fsign == +1 .AND. plan%trans) THEN
1100 istride = plan%m
1101 idist = 1
1102 ELSEIF (plan%fsign == -1 .AND. plan%trans) THEN
1103 ostride = plan%m
1104 odist = 1
1105 END IF
1106
1107 IF (plan%fsign == +1) THEN
1108 CALL dfftw_plan_many_dft(plan%fftw_plan, 1, plan%n, num_rows, zin, 0, istride, idist, &
1109 zout, 0, ostride, odist, fftw_forward, fftw_plan_type)
1110 ELSE
1111 CALL dfftw_plan_many_dft(plan%fftw_plan, 1, plan%n, num_rows, zin, 0, istride, idist, &
1112 zout, 0, ostride, odist, fftw_backward, fftw_plan_type)
1113 END IF
1114
1115!$ IF (plan%need_alt_plan) THEN
1116!$ plan%alt_num_rows = plan%m - (plan%num_threads_needed - 1)*num_rows
1117!$ IF (plan%fsign == +1) THEN
1118!$ CALL dfftw_plan_many_dft(plan%alt_fftw_plan, 1, plan%n, plan%alt_num_rows, zin, 0, istride, idist, &
1119!$ zout, 0, ostride, odist, FFTW_FORWARD, fftw_plan_type)
1120!$ ELSE
1121!$ CALL dfftw_plan_many_dft(plan%alt_fftw_plan, 1, plan%n, plan%alt_num_rows, zin, 0, istride, idist, &
1122!$ zout, 0, ostride, odist, FFTW_BACKWARD, fftw_plan_type)
1123!$ END IF
1124!$ END IF
1125
1126#else
1127 mark_used(plan)
1128 mark_used(plan_style)
1129 !MARK_USED does not work with assumed size arguments
1130 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
1131#endif
1132
1133 END SUBROUTINE fftw3_create_plan_1dm
1134
1135! **************************************************************************************************
1136!> \brief ...
1137!> \param plan ...
1138! **************************************************************************************************
1139 SUBROUTINE fftw3_destroy_plan(plan)
1140
1141 TYPE(fft_plan_type), INTENT(INOUT) :: plan
1142
1143#if defined(__FFTW3)
1144!$ IF (plan%need_alt_plan) THEN
1145!$ CALL fftw_destroy_plan(plan%alt_fftw_plan)
1146!$ END IF
1147
1148 IF (.NOT. plan%separated_plans) THEN
1149 CALL fftw_destroy_plan(plan%fftw_plan)
1150 ELSE
1151 ! If it is a separated plan then we have to destroy
1152 ! each dim plan individually
1153 CALL fftw_destroy_plan(plan%fftw_plan_nx)
1154 CALL fftw_destroy_plan(plan%fftw_plan_ny)
1155 CALL fftw_destroy_plan(plan%fftw_plan_nz)
1156 CALL fftw_destroy_plan(plan%fftw_plan_nx_r)
1157 CALL fftw_destroy_plan(plan%fftw_plan_ny_r)
1158 CALL fftw_destroy_plan(plan%fftw_plan_nz_r)
1159 END IF
1160
1161#else
1162 mark_used(plan)
1163#endif
1164
1165 END SUBROUTINE fftw3_destroy_plan
1166
1167! **************************************************************************************************
1168!> \brief ...
1169!> \param plan ...
1170!> \param zin ...
1171!> \param zout ...
1172!> \param scale ...
1173!> \param stat ...
1174! **************************************************************************************************
1175 SUBROUTINE fftw31dm(plan, zin, zout, scale, stat)
1176 TYPE(fft_plan_type), INTENT(IN) :: plan
1177 COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN), TARGET :: zin
1178 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT), &
1179 TARGET :: zout
1180 REAL(kind=dp), INTENT(IN) :: scale
1181 INTEGER, INTENT(OUT) :: stat
1182
1183 INTEGER :: in_offset, my_id, num_rows, out_offset, &
1184 scal_offset
1185 TYPE(c_ptr) :: fftw_plan
1186
1187!------------------------------------------------------------------------------
1188
1189 my_id = 0
1190 num_rows = plan%m
1191
1192!$OMP PARALLEL DEFAULT(NONE), &
1193!$OMP PRIVATE(my_id,num_rows,in_offset,out_offset,scal_offset,fftw_plan), &
1194!$OMP SHARED(zin,zout), &
1195!$OMP SHARED(plan,scale,stat)
1196!$ my_id = omp_get_thread_num()
1197
1198!$ if (my_id < plan%num_threads_needed) then
1199
1200 fftw_plan = plan%fftw_plan
1201
1202 in_offset = 1
1203 out_offset = 1
1204 scal_offset = 1
1205
1206!$ in_offset = 1 + plan%num_rows*my_id*plan%n
1207!$ out_offset = 1 + plan%num_rows*my_id*plan%n
1208!$ IF (plan%fsign == +1 .AND. plan%trans) THEN
1209!$ in_offset = 1 + plan%num_rows*my_id
1210!$ ELSEIF (plan%fsign == -1 .AND. plan%trans) THEN
1211!$ out_offset = 1 + plan%num_rows*my_id
1212!$ END IF
1213!$ scal_offset = 1 + plan%n*plan%num_rows*my_id
1214!$ IF (plan%need_alt_plan .AND. my_id == plan%num_threads_needed - 1) THEN
1215!$ num_rows = plan%alt_num_rows
1216!$ fftw_plan = plan%alt_fftw_plan
1217!$ ELSE
1218!$ num_rows = plan%num_rows
1219!$ END IF
1220
1221#if defined(__FFTW3)
1222!$OMP MASTER
1223 stat = 1
1224!$OMP END MASTER
1225 CALL dfftw_execute_dft(fftw_plan, zin(in_offset:in_offset), zout(out_offset:out_offset))
1226!$ end if
1227! all theads need to meet at this barrier
1228!$OMP BARRIER
1229!$ if (my_id < plan%num_threads_needed) then
1230 IF (scale /= 1.0_dp) CALL zdscal(plan%n*num_rows, scale, zout(scal_offset:scal_offset), 1)
1231!$ end if
1232
1233#else
1234 mark_used(plan)
1235 mark_used(scale)
1236 !MARK_USED does not work with assumed size arguments
1237 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
1238 stat = 0
1239
1240!$ else
1241!$ end if
1242
1243#endif
1244
1245!$OMP END PARALLEL
1246
1247 END SUBROUTINE fftw31dm
1248
1249 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:1140
fftw3_lib::fftw3_create_plan_3d
subroutine, public fftw3_create_plan_3d(plan, zin, zout, plan_style)
...
Definition fftw3_lib.F:709
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:956
fftw3_lib::fftw31dm
subroutine, public fftw31dm(plan, zin, zout, scale, stat)
...
Definition fftw3_lib.F:1176
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:1047
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