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