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