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 fft_lib
8 :
9 : USE fft_kinds, ONLY: dp
10 : USE fft_plan, ONLY: fft_plan_type
11 : USE fftsg_lib, ONLY: fftsg1dm,&
12 : fftsg3d,&
13 : fftsg_do_cleanup,&
14 : fftsg_do_init
15 : USE fftw3_lib, ONLY: &
16 : fft_alloc => fftw_alloc, fft_dealloc => fftw_dealloc, fftw31dm, fftw33d, &
17 : fftw3_create_plan_1d, fftw3_create_plan_3d, fftw3_destroy_plan, fftw3_do_cleanup, &
18 : fftw3_do_init, fftw3_get_lengths
19 : #include "../../base/base_uses.f90"
20 :
21 : IMPLICIT NONE
22 : PRIVATE
23 :
24 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fft_lib'
25 :
26 : INTEGER, SAVE :: fft_type = 0
27 :
28 : PUBLIC :: fft_do_cleanup, fft_do_init, fft_get_lengths, fft_create_plan_3d
29 : PUBLIC :: fft_create_plan_1d, fft_1d, fft_library, fft_3d, fft_destroy_plan
30 : PUBLIC :: fft_alloc, fft_dealloc
31 :
32 : CONTAINS
33 : ! **************************************************************************************************
34 : !> \brief Interface to FFT libraries
35 : !> \param fftlib ...
36 : !> \return ...
37 : !> \par History
38 : !> IAB 09-Jan-2009 : Modified to use fft_plan_type
39 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
40 : !> \author JGH
41 : ! **************************************************************************************************
42 11469 : FUNCTION fft_library(fftlib) RESULT(flib)
43 :
44 : CHARACTER(len=*), INTENT(IN) :: fftlib
45 : INTEGER :: flib
46 :
47 : SELECT CASE (fftlib)
48 : CASE DEFAULT
49 14 : flib = -1
50 : CASE ("FFTSG")
51 14 : flib = 1
52 : CASE ("FFTW3")
53 11469 : flib = 3
54 : END SELECT
55 :
56 11469 : END FUNCTION fft_library
57 :
58 : ! **************************************************************************************************
59 : !> \brief ...
60 : !> \param fftlib ...
61 : !> \param plan_style ...
62 : !> \param wisdom_file ...
63 : ! **************************************************************************************************
64 11469 : SUBROUTINE fft_do_init(fftlib, plan_style, wisdom_file)
65 : CHARACTER(LEN=*), INTENT(IN) :: fftlib
66 : INTEGER, INTENT(IN) :: plan_style
67 : CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
68 :
69 11469 : fft_type = fft_library(fftlib)
70 0 : SELECT CASE (fft_type)
71 : CASE DEFAULT
72 0 : CPABORT("fft_do_init")
73 : CASE (1)
74 14 : CALL fftsg_do_init()
75 : CASE (3)
76 11469 : CALL fftw3_do_init(wisdom_file, plan_style)
77 : END SELECT
78 :
79 11469 : END SUBROUTINE fft_do_init
80 :
81 : ! **************************************************************************************************
82 : !> \brief ...
83 : !> \param wisdom_file ...
84 : !> \param ionode ...
85 : ! **************************************************************************************************
86 11259 : SUBROUTINE fft_do_cleanup(wisdom_file, ionode)
87 : CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
88 : LOGICAL, INTENT(IN) :: ionode
89 :
90 0 : SELECT CASE (fft_type)
91 : CASE DEFAULT
92 0 : CPABORT("fft_do_cleanup")
93 : CASE (1)
94 14 : CALL fftsg_do_cleanup()
95 : CASE (3)
96 11259 : CALL fftw3_do_cleanup(wisdom_file, ionode)
97 : END SELECT
98 :
99 11259 : END SUBROUTINE fft_do_cleanup
100 :
101 : ! **************************************************************************************************
102 : !> \brief ...
103 : !> \param DATA ...
104 : !> \param max_length ...
105 : ! **************************************************************************************************
106 0 : SUBROUTINE fft_get_lengths(DATA, max_length)
107 : INTEGER, DIMENSION(*) :: DATA
108 : INTEGER, INTENT(INOUT) :: max_length
109 :
110 0 : CALL fftw3_get_lengths(DATA, max_length)
111 :
112 0 : END SUBROUTINE fft_get_lengths
113 :
114 : ! **************************************************************************************************
115 :
116 : ! **************************************************************************************************
117 : !> \brief ...
118 : !> \param plan ...
119 : !> \param fft_in_place ...
120 : !> \param fsign ...
121 : !> \param n ...
122 : !> \param zin ...
123 : !> \param zout ...
124 : ! **************************************************************************************************
125 62892 : SUBROUTINE fft_create_plan_3d(plan, fft_in_place, fsign, n, zin, zout)
126 :
127 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
128 : LOGICAL, INTENT(IN) :: fft_in_place
129 : INTEGER, INTENT(IN) :: fsign
130 : INTEGER, DIMENSION(3), INTENT(IN) :: n
131 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
132 :
133 62892 : plan%fsign = fsign
134 62892 : plan%fft_in_place = fft_in_place
135 251568 : plan%n_3d = n
136 62892 : !$ plan%need_alt_plan = .FALSE.
137 :
138 : ! Planning only needed for FFTW3
139 62892 : IF (fft_type == 3) THEN
140 62740 : CALL fftw3_create_plan_3d(plan, zin, zout)
141 62740 : plan%valid = .TRUE.
142 : END IF
143 :
144 62892 : END SUBROUTINE fft_create_plan_3d
145 :
146 : !
147 : ! really ugly, plan is intent out, because plan%fsign is also a status flag
148 : ! if something goes wrong, plan%fsign is set to zero, and the plan becomes invalid
149 : !
150 : ! **************************************************************************************************
151 : !> \brief ...
152 : !> \param plan ...
153 : !> \param scale ...
154 : !> \param zin ...
155 : !> \param zout ...
156 : !> \param stat ...
157 : ! **************************************************************************************************
158 663347 : SUBROUTINE fft_3d(plan, scale, zin, zout, stat)
159 : TYPE(fft_plan_type), INTENT(IN) :: plan
160 : REAL(KIND=dp), INTENT(IN) :: scale
161 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
162 : INTEGER, INTENT(OUT) :: stat
163 :
164 663347 : stat = plan%fsign
165 663347 : IF (plan%n_3d(1)*plan%n_3d(2)*plan%n_3d(3) > 0) THEN
166 0 : SELECT CASE (fft_type)
167 : CASE DEFAULT
168 0 : CPABORT("fft_3d")
169 : CASE (1)
170 2964 : CALL fftsg3d(plan%fft_in_place, stat, scale, plan%n_3d, zin, zout)
171 : CASE (3)
172 663347 : CALL fftw33d(plan, scale, zin, zout, stat)
173 : END SELECT
174 : END IF
175 : ! stat is set to zero on error, -1,+1 are OK
176 663347 : IF (stat == 0) THEN
177 0 : stat = 1
178 : ELSE
179 663347 : stat = 0
180 : END IF
181 :
182 663347 : END SUBROUTINE fft_3d
183 :
184 : ! **************************************************************************************************
185 :
186 : ! **************************************************************************************************
187 : !> \brief ...
188 : !> \param plan ...
189 : !> \param fsign ...
190 : !> \param trans_in ...
191 : !> \param trans_out ...
192 : !> \param ldx_in ...
193 : !> \param ldx_out ...
194 : !> \param n ...
195 : !> \param m ...
196 : !> \param zin ...
197 : !> \param zout ...
198 : ! **************************************************************************************************
199 463634 : SUBROUTINE fft_create_plan_1d(plan, fsign, trans_in, trans_out, ldx_in, ldx_out, n, m, zin, zout)
200 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
201 : INTEGER, INTENT(IN) :: fsign
202 : LOGICAL, INTENT(IN) :: trans_in, trans_out
203 : INTEGER, INTENT(IN) :: ldx_in, ldx_out, n, m
204 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: zin, zout
205 :
206 463634 : plan%fsign = fsign
207 463634 : plan%trans_in = trans_in
208 463634 : plan%trans_out = trans_out
209 463634 : IF (plan%trans_in) THEN
210 360605 : plan%ldx_in = ldx_in
211 360605 : plan%ldy_in = n
212 : ELSE
213 103029 : plan%ldx_in = ldx_in
214 103029 : plan%ldy_in = m
215 : END IF
216 463634 : IF (plan%trans_out) THEN
217 360605 : plan%ldx_out = ldx_out
218 360605 : plan%ldy_out = n
219 : ELSE
220 103029 : plan%ldx_out = ldx_out
221 103029 : plan%ldy_out = m
222 : END IF
223 :
224 463634 : plan%n = n
225 463634 : plan%m = m
226 463634 : !$ plan%need_alt_plan = .FALSE.
227 :
228 : ! Planning only needed for FFTW3
229 463634 : IF ((fft_type == 3) .AND. (n*m /= 0)) THEN
230 462278 : CALL fftw3_create_plan_1d(plan, zin, zout)
231 462278 : plan%valid = .TRUE.
232 : ELSE
233 1356 : plan%valid = .FALSE.
234 : END IF
235 :
236 463634 : END SUBROUTINE fft_create_plan_1d
237 :
238 : ! **************************************************************************************************
239 : !> \brief ...
240 : !> \param plan ...
241 : ! **************************************************************************************************
242 706910 : SUBROUTINE fft_destroy_plan(plan)
243 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
244 :
245 : ! Planning only needed for FFTW3
246 :
247 706910 : IF (plan%valid .AND. fft_type == 3) THEN
248 525018 : CALL fftw3_destroy_plan(plan)
249 525018 : plan%valid = .FALSE.
250 : END IF
251 :
252 706910 : END SUBROUTINE fft_destroy_plan
253 :
254 : ! **************************************************************************************************
255 : !> \brief ...
256 : !> \param plan ...
257 : !> \param zin ...
258 : !> \param zout ...
259 : !> \param scale ...
260 : !> \param stat ...
261 : ! **************************************************************************************************
262 19999300 : SUBROUTINE fft_1d(plan, zin, zout, scale, stat)
263 : TYPE(fft_plan_type), INTENT(IN) :: plan
264 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
265 : REAL(KIND=dp), INTENT(IN) :: scale
266 : INTEGER, INTENT(OUT) :: stat
267 :
268 19999300 : stat = plan%fsign
269 19999300 : IF (plan%n*plan%m > 0) THEN
270 0 : SELECT CASE (fft_type)
271 : CASE DEFAULT
272 0 : CPABORT("fft_1d")
273 : CASE (1)
274 : CALL fftsg1dm(stat, plan%trans_in, plan%trans_out, plan%n, plan%m, &
275 30238 : plan%ldx_in, plan%ldy_in, plan%ldx_out, plan%ldy_out, zin, zout, scale)
276 : CASE (3)
277 19999300 : CALL fftw31dm(plan, zin, zout, scale, stat)
278 : END SELECT
279 : END IF
280 : ! stat is set to zero on error, -1,+1 are OK
281 19999300 : IF (stat == 0) THEN
282 0 : stat = 1
283 : ELSE
284 19999300 : stat = 0
285 : END IF
286 :
287 19999300 : END SUBROUTINE fft_1d
288 :
289 : END MODULE fft_lib
|
