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