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