(git:ec11232)
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 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 defined(__FFTW3_UNALIGNED)
742 fftw_plan_type = fftw_plan_type + fftw_unaligned
743#endif
744
745 IF (plan%fsign == +1) THEN
746 fft_direction = fftw_forward
747 ELSE
748 fft_direction = fftw_backward
749 END IF
750
751 n1 = plan%n_3d(1)
752 n2 = plan%n_3d(2)
753 n3 = plan%n_3d(3)
754
755 nt = 1
756!$OMP PARALLEL DEFAULT(NONE) SHARED(nt)
757!$OMP MASTER
758!$ nt = omp_get_num_threads()
759!$OMP END MASTER
760!$OMP END PARALLEL
761
762 IF ((.NOT. fftw3_is_guru_supported()) .OR. &
763 (.NOT. plan_style == 1) .OR. &
764 (n1 < 256 .AND. n2 < 256 .AND. n3 < 256 .AND. nt == 1)) THEN
765 ! If the plan type is MEASURE, PATIENT and EXHAUSTIVE or
766 ! the grid size is small (and we are single-threaded) then
767 ! FFTW3 does a better job than handmade optimization
768 ! so plan a single 3D FFT which will execute using all the threads
769
770 plan%separated_plans = .false.
771!$ CALL fftw_plan_with_nthreads(nt)
772
773 IF (plan%fft_in_place) THEN
774 plan%fftw_plan = fftw_plan_dft_3d(n3, n2, n1, zin, zin, fft_direction, fftw_plan_type)
775 ELSE
776 plan%fftw_plan = fftw_plan_dft_3d(n3, n2, n1, zin, zout, fft_direction, fftw_plan_type)
777 END IF
778 ELSE
779 ALLOCATE (tmp(n1*n2*n3))
780 ! ************************* PLANS WITH TRANSPOSITIONS ****************************
781 ! In the cases described above, we manually thread each stage of the 3D FFT.
782 !
783 ! The following plans replace the 3D FFT call by running 1D FFTW across all
784 ! 3 directions of the array.
785 !
786 ! Output of FFTW is transposed to ensure that the next round of FFTW access
787 ! contiguous information.
788 !
789 ! Assuming the input matrix is M(n3,n2,n1), FFTW/Transp are :
790 ! M(n3,n2,n1) -> fftw(x) -> M(n3,n1,n2) -> fftw(y) -> M(n1,n2,n3) -> fftw(z) -> M(n1,n2,n3)
791 ! Notice that last matrix is transposed in the Z axis. A DO-loop in the execute routine
792 ! will perform the final transposition. Performance evaluation showed that using an external
793 ! DO loop to do the final transposition performed better than directly transposing the output.
794 ! However, this might vary depending on the compiler/platform, so a potential tuning spot
795 ! is to perform the final transposition within the fftw library rather than using the external loop
796 ! See comments below in Z-FFT for how to transpose the output to avoid the final DO loop.
797 !
798 ! Doc. for the Guru interface is in http://www.fftw.org/doc/Guru-Interface.html
799 !
800 ! OpenMP : Work is distributed on the Z plane.
801 ! All transpositions are out-of-place to facilitate multi-threading
802 !
803 !!!! Plan for X : M(n3,n2,n1) -> fftw(x) -> M(n3,n1,n2)
804 CALL fftw3_compute_rows_per_th(n3, nt, rows_per_th, rows_per_th_r, &
805 th_plana, th_planb)
806
807 dim_n(1) = n1
808 dim_istride(1) = 1
809 dim_ostride(1) = n2
810 howmany_n(1) = n2
811 howmany_n(2) = rows_per_th
812 howmany_istride(1) = n1
813 howmany_istride(2) = n1*n2
814 howmany_ostride(1) = 1
815 howmany_ostride(2) = n1*n2
816 CALL fftw3_create_3d_plans(plan%fftw_plan_nx, plan%fftw_plan_nx_r, &
817 dim_n, dim_istride, dim_ostride, howmany_n, &
818 howmany_istride, howmany_ostride, &
819 zin, tmp, &
820 fft_direction, fftw_plan_type, rows_per_th, &
821 rows_per_th_r)
822
823 !!!! Plan for Y : M(n3,n1,n2) -> fftw(y) -> M(n1,n2,n3)
824 CALL fftw3_compute_rows_per_th(n3, nt, rows_per_th, rows_per_th_r, &
825 th_plana, th_planb)
826 dim_n(1) = n2
827 dim_istride(1) = 1
828 dim_ostride(1) = n3
829 howmany_n(1) = n1
830 howmany_n(2) = rows_per_th
831 howmany_istride(1) = n2
832 howmany_istride(2) = n1*n2
833 !!! transposed Z axis on output
834 howmany_ostride(1) = n2*n3
835 howmany_ostride(2) = 1
836
837 CALL fftw3_create_3d_plans(plan%fftw_plan_ny, plan%fftw_plan_ny_r, &
838 dim_n, dim_istride, dim_ostride, &
839 howmany_n, howmany_istride, howmany_ostride, &
840 tmp, zin, &
841 fft_direction, fftw_plan_type, rows_per_th, &
842 rows_per_th_r)
843
844 !!!! Plan for Z : M(n1,n2,n3) -> fftw(z) -> M(n1,n2,n3)
845 CALL fftw3_compute_rows_per_th(n1, nt, rows_per_th, rows_per_th_r, &
846 th_plana, th_planb)
847 dim_n(1) = n3
848 dim_istride(1) = 1
849 dim_ostride(1) = 1 ! To transpose: n2*n1
850 howmany_n(1) = n2
851 howmany_n(2) = rows_per_th
852 howmany_istride(1) = n3
853 howmany_istride(2) = n2*n3
854 howmany_ostride(1) = n3 ! To transpose: n1
855 howmany_ostride(2) = n2*n3 ! To transpose: 1
856
857 CALL fftw3_create_3d_plans(plan%fftw_plan_nz, plan%fftw_plan_nz_r, &
858 dim_n, dim_istride, dim_ostride, &
859 howmany_n, howmany_istride, howmany_ostride, &
860 zin, tmp, &
861 fft_direction, fftw_plan_type, rows_per_th, &
862 rows_per_th_r)
863
864 plan%separated_plans = .true.
865
866 DEALLOCATE (tmp)
867 END IF
868
869#else
870 mark_used(plan)
871 mark_used(plan_style)
872 !MARK_USED does not work with assumed size arguments
873 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
874#endif
875
876 END SUBROUTINE fftw3_create_plan_3d
877
878! **************************************************************************************************
879
880! **************************************************************************************************
881!> \brief ...
882!> \param plan ...
883!> \param plan_r ...
884!> \param split_dim ...
885!> \param nt ...
886!> \param tid ...
887!> \param input ...
888!> \param istride ...
889!> \param output ...
890!> \param ostride ...
891! **************************************************************************************************
892 SUBROUTINE fftw3_workshare_execute_dft(plan, plan_r, split_dim, nt, tid, &
893 input, istride, output, ostride)
894
895 INTEGER, INTENT(IN) :: split_dim, nt, tid
896 INTEGER, INTENT(IN) :: istride, ostride
897 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: input, output
898 TYPE(c_ptr) :: plan, plan_r
899#if defined(__FFTW3)
900 INTEGER :: i_off, o_off
901 INTEGER :: th_plana, th_planb
902 INTEGER :: rows_per_thread, rows_per_thread_r
903
904 CALL fftw3_compute_rows_per_th(split_dim, nt, rows_per_thread, &
905 rows_per_thread_r, &
906 th_plana, th_planb)
907
908 IF (th_planb > 0) THEN
909 IF (tid < th_plana) THEN
910 i_off = (tid)*(istride*(rows_per_thread)) + 1
911 o_off = (tid)*(ostride*(rows_per_thread)) + 1
912 IF (rows_per_thread > 0) THEN
913 CALL fftw_execute_dft(plan, input(i_off), &
914 output(o_off))
915 END IF
916 ELSE IF ((tid - th_plana) < th_planb) THEN
917
918 i_off = (th_plana)*istride*(rows_per_thread) + &
919 (tid - th_plana)*istride*(rows_per_thread_r) + 1
920 o_off = (th_plana)*ostride*(rows_per_thread) + &
921 (tid - th_plana)*ostride*(rows_per_thread_r) + 1
922
923 CALL fftw_execute_dft(plan_r, input(i_off), &
924 output(o_off))
925 END IF
926
927 ELSE
928 i_off = (tid)*(istride*(rows_per_thread)) + 1
929 o_off = (tid)*(ostride*(rows_per_thread)) + 1
930
931 CALL fftw_execute_dft(plan, input(i_off), &
932 output(o_off))
933
934 END IF
935#else
936 mark_used(plan)
937 mark_used(plan_r)
938 mark_used(split_dim)
939 mark_used(nt)
940 mark_used(tid)
941 mark_used(istride)
942 mark_used(ostride)
943 !MARK_USED does not work with assumed size arguments
944 IF (.false.) then; do; IF (abs(input(1)) > abs(output(1))) exit; END do; END IF
945#endif
946
947 END SUBROUTINE fftw3_workshare_execute_dft
948
949! **************************************************************************************************
950
951! **************************************************************************************************
952!> \brief ...
953!> \param plan ...
954!> \param scale ...
955!> \param zin ...
956!> \param zout ...
957!> \param stat ...
958! **************************************************************************************************
959 SUBROUTINE fftw33d(plan, scale, zin, zout, stat)
960
961 TYPE(fft_plan_type), INTENT(IN) :: plan
962 REAL(kind=dp), INTENT(IN) :: scale
963 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT), TARGET:: zin
964 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT), TARGET:: zout
965 INTEGER, INTENT(OUT) :: stat
966#if defined(__FFTW3)
967 COMPLEX(KIND=dp), POINTER :: xout(:)
968 COMPLEX(KIND=dp), ALLOCATABLE :: tmp1(:)
969 INTEGER :: n1, n2, n3
970 INTEGER :: tid, nt
971 INTEGER :: i, j, k
972
973 n1 = plan%n_3d(1)
974 n2 = plan%n_3d(2)
975 n3 = plan%n_3d(3)
976
977 stat = 1
978
979 ! We use a POINTER to the output array to avoid duplicating code
980 IF (plan%fft_in_place) THEN
981 xout => zin(:n1*n2*n3)
982 ELSE
983 xout => zout(:n1*n2*n3)
984 END IF
985
986 ! Either compute the full 3D FFT using a multithreaded plan
987 IF (.NOT. plan%separated_plans) THEN
988 CALL fftw_execute_dft(plan%fftw_plan, zin, xout)
989 ELSE
990 ! Or use the 3 stage FFT scheme described in fftw3_create_plan_3d
991 ALLOCATE (tmp1(n1*n2*n3)) ! Temporary vector used for transpositions
992 !$OMP PARALLEL DEFAULT(NONE) PRIVATE(tid,nt,i,j,k) SHARED(zin,tmp1,n1,n2,n3,plan,xout)
993 tid = 0
994 nt = 1
995
996!$ tid = omp_get_thread_num()
997!$ nt = omp_get_num_threads()
998 CALL fftw3_workshare_execute_dft(plan%fftw_plan_nx, plan%fftw_plan_nx_r, &
999 n3, nt, tid, &
1000 zin, n1*n2, tmp1, n1*n2)
1001
1002 !$OMP BARRIER
1003 CALL fftw3_workshare_execute_dft(plan%fftw_plan_ny, plan%fftw_plan_ny_r, &
1004 n3, nt, tid, &
1005 tmp1, n1*n2, xout, 1)
1006 !$OMP BARRIER
1007 CALL fftw3_workshare_execute_dft(plan%fftw_plan_nz, plan%fftw_plan_nz_r, &
1008 n1, nt, tid, &
1009 xout, n2*n3, tmp1, n2*n3)
1010 !$OMP BARRIER
1011
1012 !$OMP DO COLLAPSE(3)
1013 DO i = 1, n1
1014 DO j = 1, n2
1015 DO k = 1, n3
1016 xout((i - 1) + (j - 1)*n1 + (k - 1)*n1*n2 + 1) = &
1017 tmp1((k - 1) + (j - 1)*n3 + (i - 1)*n3*n2 + 1)
1018 END DO
1019 END DO
1020 END DO
1021 !$OMP END DO
1022
1023 !$OMP END PARALLEL
1024 END IF
1025
1026 IF (scale /= 1.0_dp) THEN
1027 CALL zdscal(n1*n2*n3, scale, xout, 1)
1028 END IF
1029
1030#else
1031 mark_used(plan)
1032 mark_used(scale)
1033 !MARK_USED does not work with assumed size arguments
1034 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
1035 stat = 0
1036
1037#endif
1038
1039 END SUBROUTINE fftw33d
1040
1041! **************************************************************************************************
1042
1043! **************************************************************************************************
1044!> \brief ...
1045!> \param plan ...
1046!> \param zin ...
1047!> \param zout ...
1048!> \param plan_style ...
1049! **************************************************************************************************
1050 SUBROUTINE fftw3_create_plan_1dm(plan, zin, zout, plan_style)
1051
1052 IMPLICIT NONE
1053
1054 TYPE(fft_plan_type), INTENT(INOUT) :: plan
1055 COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: zin
1056 COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: zout
1057 INTEGER, INTENT(IN) :: plan_style
1058#if defined(__FFTW3)
1059 INTEGER :: istride, idist, ostride, odist, num_threads, num_rows
1060
1061 INTEGER :: fftw_plan_type
1062 SELECT CASE (plan_style)
1063 CASE (1)
1064 fftw_plan_type = fftw_estimate
1065 CASE (2)
1066 fftw_plan_type = fftw_measure
1067 CASE (3)
1068 fftw_plan_type = fftw_patient
1069 CASE (4)
1070 fftw_plan_type = fftw_exhaustive
1071 CASE DEFAULT
1072 cpabort("fftw3_create_plan_1dm")
1073 END SELECT
1074
1075#if defined(__FFTW3_UNALIGNED)
1076 fftw_plan_type = fftw_plan_type + fftw_unaligned
1077#endif
1078 num_threads = 1
1079 plan%separated_plans = .false.
1080!$OMP PARALLEL DEFAULT(NONE), &
1081!$OMP SHARED(NUM_THREADS)
1082!$OMP MASTER
1083!$ num_threads = omp_get_num_threads()
1084!$OMP END MASTER
1085!$OMP END PARALLEL
1086
1087 num_rows = plan%m/num_threads
1088!$ plan%num_threads_needed = num_threads
1089
1090! Check for number of rows less than num_threads
1091!$ IF (plan%m < num_threads) THEN
1092!$ num_rows = 1
1093!$ plan%num_threads_needed = plan%m
1094!$ END IF
1095
1096! Check for total number of rows not divisible by num_threads
1097!$ IF (num_rows*plan%num_threads_needed /= plan%m) THEN
1098!$ plan%need_alt_plan = .TRUE.
1099!$ END IF
1100
1101!$ plan%num_rows = num_rows
1102 istride = 1
1103 idist = plan%n
1104 ostride = 1
1105 odist = plan%n
1106 IF (plan%fsign == +1 .AND. plan%trans) THEN
1107 istride = plan%m
1108 idist = 1
1109 ELSEIF (plan%fsign == -1 .AND. plan%trans) THEN
1110 ostride = plan%m
1111 odist = 1
1112 END IF
1113
1114 IF (plan%fsign == +1) THEN
1115 CALL dfftw_plan_many_dft(plan%fftw_plan, 1, plan%n, num_rows, zin, 0, istride, idist, &
1116 zout, 0, ostride, odist, fftw_forward, fftw_plan_type)
1117 ELSE
1118 CALL dfftw_plan_many_dft(plan%fftw_plan, 1, plan%n, num_rows, zin, 0, istride, idist, &
1119 zout, 0, ostride, odist, fftw_backward, fftw_plan_type)
1120 END IF
1121
1122!$ IF (plan%need_alt_plan) THEN
1123!$ plan%alt_num_rows = plan%m - (plan%num_threads_needed - 1)*num_rows
1124!$ IF (plan%fsign == +1) THEN
1125!$ CALL dfftw_plan_many_dft(plan%alt_fftw_plan, 1, plan%n, plan%alt_num_rows, zin, 0, istride, idist, &
1126!$ zout, 0, ostride, odist, FFTW_FORWARD, fftw_plan_type)
1127!$ ELSE
1128!$ CALL dfftw_plan_many_dft(plan%alt_fftw_plan, 1, plan%n, plan%alt_num_rows, zin, 0, istride, idist, &
1129!$ zout, 0, ostride, odist, FFTW_BACKWARD, fftw_plan_type)
1130!$ END IF
1131!$ END IF
1132
1133#else
1134 mark_used(plan)
1135 mark_used(plan_style)
1136 !MARK_USED does not work with assumed size arguments
1137 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
1138#endif
1139
1140 END SUBROUTINE fftw3_create_plan_1dm
1141
1142! **************************************************************************************************
1143!> \brief ...
1144!> \param plan ...
1145! **************************************************************************************************
1146 SUBROUTINE fftw3_destroy_plan(plan)
1147
1148 TYPE(fft_plan_type), INTENT(INOUT) :: plan
1149
1150#if defined(__FFTW3)
1151!$ IF (plan%need_alt_plan) THEN
1152!$ CALL fftw_destroy_plan(plan%alt_fftw_plan)
1153!$ END IF
1154
1155 IF (.NOT. plan%separated_plans) THEN
1156 CALL fftw_destroy_plan(plan%fftw_plan)
1157 ELSE
1158 ! If it is a separated plan then we have to destroy
1159 ! each dim plan individually
1160 CALL fftw_destroy_plan(plan%fftw_plan_nx)
1161 CALL fftw_destroy_plan(plan%fftw_plan_ny)
1162 CALL fftw_destroy_plan(plan%fftw_plan_nz)
1163 CALL fftw_destroy_plan(plan%fftw_plan_nx_r)
1164 CALL fftw_destroy_plan(plan%fftw_plan_ny_r)
1165 CALL fftw_destroy_plan(plan%fftw_plan_nz_r)
1166 END IF
1167
1168#else
1169 mark_used(plan)
1170#endif
1171
1172 END SUBROUTINE fftw3_destroy_plan
1173
1174! **************************************************************************************************
1175!> \brief ...
1176!> \param plan ...
1177!> \param zin ...
1178!> \param zout ...
1179!> \param scale ...
1180!> \param stat ...
1181! **************************************************************************************************
1182 SUBROUTINE fftw31dm(plan, zin, zout, scale, stat)
1183 TYPE(fft_plan_type), INTENT(IN) :: plan
1184 COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN), TARGET :: zin
1185 COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT), &
1186 TARGET :: zout
1187 REAL(kind=dp), INTENT(IN) :: scale
1188 INTEGER, INTENT(OUT) :: stat
1189
1190 INTEGER :: in_offset, my_id, num_rows, out_offset, &
1191 scal_offset
1192 TYPE(c_ptr) :: fftw_plan
1193
1194!------------------------------------------------------------------------------
1195
1196 my_id = 0
1197 num_rows = plan%m
1198
1199!$OMP PARALLEL DEFAULT(NONE), &
1200!$OMP PRIVATE(my_id,num_rows,in_offset,out_offset,scal_offset,fftw_plan), &
1201!$OMP SHARED(zin,zout), &
1202!$OMP SHARED(plan,scale,stat)
1203!$ my_id = omp_get_thread_num()
1204
1205!$ if (my_id < plan%num_threads_needed) then
1206
1207 fftw_plan = plan%fftw_plan
1208
1209 in_offset = 1
1210 out_offset = 1
1211 scal_offset = 1
1212
1213!$ in_offset = 1 + plan%num_rows*my_id*plan%n
1214!$ out_offset = 1 + plan%num_rows*my_id*plan%n
1215!$ IF (plan%fsign == +1 .AND. plan%trans) THEN
1216!$ in_offset = 1 + plan%num_rows*my_id
1217!$ ELSEIF (plan%fsign == -1 .AND. plan%trans) THEN
1218!$ out_offset = 1 + plan%num_rows*my_id
1219!$ END IF
1220!$ scal_offset = 1 + plan%n*plan%num_rows*my_id
1221!$ IF (plan%need_alt_plan .AND. my_id == plan%num_threads_needed - 1) THEN
1222!$ num_rows = plan%alt_num_rows
1223!$ fftw_plan = plan%alt_fftw_plan
1224!$ ELSE
1225!$ num_rows = plan%num_rows
1226!$ END IF
1227
1228#if defined(__FFTW3)
1229!$OMP MASTER
1230 stat = 1
1231!$OMP END MASTER
1232 CALL dfftw_execute_dft(fftw_plan, zin(in_offset:in_offset), zout(out_offset:out_offset))
1233!$ end if
1234! all theads need to meet at this barrier
1235!$OMP BARRIER
1236!$ if (my_id < plan%num_threads_needed) then
1237 IF (scale /= 1.0_dp) CALL zdscal(plan%n*num_rows, scale, zout(scal_offset:scal_offset), 1)
1238!$ end if
1239
1240#else
1241 mark_used(plan)
1242 mark_used(scale)
1243 !MARK_USED does not work with assumed size arguments
1244 IF (.false.) then; do; IF (abs(zin(1)) > abs(zout(1))) exit; END do; END IF
1245 stat = 0
1246
1247!$ else
1248!$ end if
1249
1250#endif
1251
1252!$OMP END PARALLEL
1253
1254 END SUBROUTINE fftw31dm
1255
1256 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:1147
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:960
fftw3_lib::fftw31dm
subroutine, public fftw31dm(plan, zin, zout, scale, stat)
...
Definition fftw3_lib.F:1183
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:1051
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