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