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