Branch data Line data Source code
1 : : !!****m* ABINIT/m_dfti
2 : : !! NAME
3 : : !! m_dfti
4 : : !!
5 : : !! FUNCTION
6 : : !! This module provides wrappers for the MKL DFTI routines: in-place and out-of-place version.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2009-2026 ABINIT group (MG)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! NOTES
15 : : !! 1) MPI parallelism is not supported
16 : : !! 2) For better performance the FFT divisions should contain small factors (/2, 3, 5, 7, 11, 13/)
17 : : !! see http://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_userguide_lnx/index.htm
18 : : !!
19 : : !! SOURCE
20 : :
21 : : #if defined HAVE_CONFIG_H
22 : : #include "config.h"
23 : : #endif
24 : :
25 : : #include "abi_common.h"
26 : :
27 : : #ifdef HAVE_DFTI
28 : :
29 : : ! Include MKL_DFTI fortran module
30 : : #include "mkl_dfti.f90"
31 : :
32 : : ! Macros for template files.
33 : : #define FFTLIB "DFTI"
34 : : #define FFT_PREF_fftrisc dfti_fftrisc
35 : : #define FFT_PREF_fftrisc_mixprec dfti_fftrisc_mixprec
36 : : #define FFT_PREF_fftpad dfti_fftpad
37 : : #define FFT_PREF_fftug_dp dfti_fftug_dp
38 : : #define FFT_PREF_fftur_dp dfti_fftur_dp
39 : : #define FFT_PREF_fftug dfti_fftug
40 : : #define FFT_PREF_fftur dfti_fftur
41 : : #define SPAWN_THREADS_HERE(ndat, nthreads) dfti_spawn_threads_here(ndat, nthreads)
42 : :
43 : : #define FFT_DOUBLE 1
44 : : #define FFT_SINGLE 2
45 : : #define FFT_MIXPREC 3
46 : :
47 : : #endif
48 : :
49 : : MODULE m_dfti
50 : :
51 : : use, intrinsic :: iso_c_binding
52 : : use defs_basis
53 : : use m_abicore
54 : : use m_errors
55 : : use m_xomp
56 : : use m_cgtools
57 : : use m_cplxtools
58 : : use m_fftcore
59 : : use m_fft_mesh
60 : : #ifdef HAVE_DFTI
61 : : use MKL_DFTI
62 : : #endif
63 : :
64 : : use m_fstrings, only : basename, strcat, int2char10, itoa, sjoin
65 : : use m_hide_blas, only : xcopy
66 : :
67 : : implicit none
68 : :
69 : : private
70 : :
71 : : ! Entry points for client code
72 : : public :: dfti_seqfourdp ! 3D FFT of lengths nx, ny, nz. Mainly used for densities or potentials.
73 : : public :: dfti_seqfourwf ! FFT transform of wavefunctions (high-level interface).
74 : : public :: dfti_fftrisc
75 : : public :: dfti_fftrisc_mixprec ! Mixed precision version of fftrisc: input/output in dp, computation done in sp.
76 : : public :: dfti_fftug ! G-->R, 3D zero-padded FFT of lengths nx, ny, nz. Mainly used for wavefunctions
77 : : public :: dfti_fftur ! R-->G, 3D zero-padded FFT of lengths nx, ny, nz. Mainly used for wavefunctions
78 : :
79 : : ! Low-level routines.
80 : : public :: dfti_r2c_op ! Real to complex transform (out-of-place version).
81 : : public :: dfti_c2r_op ! Complex to real transform (out-of-place version).
82 : : public :: dfti_c2c_op ! complex to complex transform (out-of-place version).
83 : : public :: dfti_c2c_ip ! complex to complex transform (in-place version).
84 : : public :: dfti_many_dft_op ! Driver routine for many out-of-place 3D complex-to-complex FFTs.
85 : : public :: dfti_many_dft_ip ! Driver routine for many in-place 3D complex-to-complex FFTs.
86 : : public :: dfti_fftpad ! Driver routines for zero-padded FFT of wavefunctions.
87 : :
88 : : !FIXME I don't know why gcc does not recognize this one
89 : : ! Perhaps I have to provide an interfaces for fofr(:,:,:,:)
90 : : public :: dfti_fftpad_dp ! Driver routines for zero-padded FFT of wavefunctions.
91 : : public :: dfti_fftug_dp ! Driver routines for zero-padded FFT of wavefunctions.
92 : : public :: dfti_use_lib_threads
93 : : !!***
94 : :
95 : : interface dfti_fftrisc
96 : : module procedure dfti_fftrisc_sp
97 : : module procedure dfti_fftrisc_dp
98 : : end interface dfti_fftrisc
99 : :
100 : : interface dfti_fftug
101 : : module procedure dfti_fftug_dp
102 : : module procedure dfti_fftug_spc
103 : : module procedure dfti_fftug_dpc
104 : : end interface dfti_fftug
105 : :
106 : : interface dfti_fftur
107 : : module procedure dfti_fftur_dp
108 : : module procedure dfti_fftur_spc
109 : : module procedure dfti_fftur_dpc
110 : : end interface dfti_fftur
111 : :
112 : : interface dfti_r2c_op
113 : : module procedure dfti_r2c_op_dp
114 : : module procedure dfti_r2c_op_dpc
115 : : end interface dfti_r2c_op
116 : :
117 : : interface dfti_c2r_op
118 : : module procedure dfti_c2r_op_dp
119 : : module procedure dfti_c2r_op_dpc
120 : : end interface dfti_c2r_op
121 : :
122 : : interface dfti_c2c_op
123 : : module procedure dfti_c2c_op_spc
124 : : module procedure dfti_c2c_op_dpc
125 : : end interface dfti_c2c_op
126 : :
127 : : interface dfti_c2c_ip
128 : : module procedure dfti_c2c_ip_spc
129 : : module procedure dfti_c2c_ip_dpc
130 : : end interface dfti_c2c_ip
131 : :
132 : : !interface dfti_many_dft_op
133 : : ! module procedure dfti_many_dft_op
134 : : ! module procedure dfti_many_dft_op
135 : : !end interface dfti_many_dft_op
136 : :
137 : : !interface dfti_many_dft_ip
138 : : ! module procedure dfti_many_dft_ip
139 : : ! module procedure dfti_many_dft_ip
140 : : !end interface dfti_many_dft_ip
141 : :
142 : : interface dfti_fftpad
143 : : module procedure dfti_fftpad_dp
144 : : module procedure dfti_fftpad_spc
145 : : module procedure dfti_fftpad_dpc
146 : : end interface dfti_fftpad
147 : :
148 : : logical,private,save :: USE_LIB_THREADS = .FALSE.
149 : :
150 : : #ifdef HAVE_DFTI
151 : : ! dfti_alloc_* allocates arrays aligned on DFTI_DEFAULT_ALIGNMENT boundaries.
152 : : integer(C_INT),private,parameter :: DFTI_DEFAULT_ALIGNMENT_SP = 64
153 : : integer(C_INT),private,parameter :: DFTI_DEFAULT_ALIGNMENT_DP = 64
154 : :
155 : : interface dfti_alloc_real
156 : : !module procedure dfti_alloc_real_sp
157 : : module procedure dfti_alloc_real_dp
158 : : end interface dfti_alloc_real
159 : :
160 : : interface dfti_alloc_complex
161 : : module procedure dfti_alloc_complex_spc
162 : : module procedure dfti_alloc_complex_dpc
163 : : end interface dfti_alloc_complex
164 : :
165 : : ! Fortran binding for MKL_malloc
166 : : interface mkl_malloc
167 : : type(C_PTR) function mkl_malloc(alloc_size, alignment) bind(C, name='MKL_malloc')
168 : : import
169 : : integer(C_SIZE_T), value :: alloc_size
170 : : integer(C_INT), value :: alignment
171 : : end function mkl_malloc
172 : : end interface mkl_malloc
173 : :
174 : : ! Fortran binding for MKL_free
175 : : interface dfti_free
176 : : subroutine mkl_free(cptr) bind(C, name='MKL_free')
177 : : import
178 : : type(C_PTR), value :: cptr
179 : : end subroutine mkl_free
180 : : end interface dfti_free
181 : : #endif
182 : :
183 : : !----------------------------------------------------------------------
184 : :
185 : : CONTAINS !===========================================================
186 : : !!***
187 : :
188 : : !!****f* m_dfti/dfti_seqfourdp
189 : : !! NAME
190 : : !! dfti_seqfourdp
191 : : !!
192 : : !! FUNCTION
193 : : !! Driver routine for 3D FFT of lengths nx, ny, nz. Mainly used for densities or potentials.
194 : : !! FFT Transform is out-of-place
195 : : !!
196 : : !! INPUTS
197 : : !! cplex=1 if fofr is real, 2 if fofr is complex
198 : : !! nx,ny,nz=Number of point along the three directions.
199 : : !! ldx,ldy,ldz=Leading dimensions of the array.
200 : : !! ndat = Number of FFTS
201 : : !! isign= +1 : fofg(G) => fofr(R);
202 : : !! -1 : fofr(R) => fofg(G)
203 : : !! fofg(2,ldx*ldy*ldz*ndat)=The array to be transformed.
204 : : !!
205 : : !! OUTPUT
206 : : !! fofr(cplex,ldx*ldy*ldz*ndat)=The FFT of fofg
207 : : !!
208 : : !! SOURCE
209 : :
210 : 1263123 : subroutine dfti_seqfourdp(cplex,nx,ny,nz,ldx,ldy,ldz,ndat,isign,fofg,fofr)
211 : :
212 : : !Arguments ------------------------------------
213 : : !scalars
214 : : integer,intent(in) :: cplex,nx,ny,nz,ldx,ldy,ldz,ndat,isign
215 : : !arrays
216 : : real(dp),intent(inout) :: fofg(2*ldx*ldy*ldz*ndat)
217 : : real(dp),intent(inout) :: fofr(cplex*ldx*ldy*ldz*ndat)
218 : :
219 : : !Local variables-------------------------------
220 : : !scalars
221 : : integer,parameter :: iscale1 = 1
222 : : integer :: ii,jj
223 : 1263123 : complex(sp), allocatable :: work_sp(:)
224 : : ! *************************************************************************
225 : :
226 : 498031 : select case (cplex)
227 : : case (2)
228 : : ! Complex to Complex.
229 [ + + ]: 498031 : if (fftcore_mixprec == 1) then
230 : : ! Mixed precision: copy in + in-place + copyout
231 [ - + + - : 100752 : ABI_MALLOC(work_sp, (ldx*ldy*ldz*ndat))
+ - + - -
+ - + ]
232 [ + + ]: 33584 : if (isign == +1) then
233 [ + + ]: 84745007 : work_sp(:) = cmplx(fofg(1::2), fofg(2::2), kind=sp)
234 [ + - ]: 9901 : else if (isign == -1) then
235 [ + + ]: 77821821 : work_sp(:) = cmplx(fofr(1::2), fofr(2::2), kind=sp)
236 : : else
237 : 0 : ABI_BUG("Wrong isign")
238 : : end if
239 : :
240 : 33584 : call dfti_c2c_ip_spc(nx, ny, nz, ldx, ldy, ldz, ndat, iscale1, isign, work_sp)
241 : :
242 [ + + ]: 33584 : if (isign == +1) then
243 : : jj = 1
244 [ + + ]: 84745007 : do ii=1,ldx*ldy*ldz*ndat
245 : 84721324 : fofr(jj) = real(work_sp(ii), kind=dp)
246 : 84721324 : fofr(jj+1) = aimag(work_sp(ii))
247 : 84745007 : jj = jj + 2
248 : : end do
249 [ + - ]: 9901 : else if (isign == -1) then
250 : : jj = 1
251 [ + + ]: 77821821 : do ii=1,ldx*ldy*ldz*ndat
252 : 77811920 : fofg(jj) = real(work_sp(ii), kind=dp)
253 : 77811920 : fofg(jj+1) = aimag(work_sp(ii))
254 : 77821821 : jj = jj + 2
255 : : end do
256 : : end if
257 [ - + ]: 33584 : ABI_FREE(work_sp)
258 : :
259 : : else
260 : : ! double precision version.
261 : 699204 : select case (isign)
262 : : case (+1)
263 : 234757 : call dfti_many_dft_op(nx,ny,nz,ldx,ldy,ldz,ndat,isign,fofg,fofr)
264 : : case (-1) ! -1
265 : 229690 : call dfti_many_dft_op(nx,ny,nz,ldx,ldy,ldz,ndat,isign,fofr,fofg)
266 : : case default
267 [ + + - ]: 464447 : ABI_BUG("Wrong isign")
268 : : end select
269 : : end if
270 : :
271 : : case (1)
272 : : ! Real case.
273 : :
274 : : ! MG: June 24. 2025
275 : : ! dfti_seqfourdp does not work as expected when cplex= 1 and ngfft(1:3) != ngfft(4:6)
276 : : ! very likely due to the use of r->c, c->r transforms.
277 : : ! I don't know if it's a bug as the error seems to depend on the mkl version.
278 : : ! To bypass this problem, we change the params on the fly so that ngfft(1:3) == ngfft(4:6)
279 : : ! when FFT_DFTI is used.
280 : : ! Note however that we never call fourdp with ngfft(1:3) != ngftt(4:6) so this is not a serious problem.
281 : : ! An additional check is done inside dfti_seqfourdp
282 : :
283 [ + - + - : 765092 : if (nx /= ldx .or. ny /= ldy .or. nz /= ldz) then
- + ]
284 : 0 : ABI_ERROR("dfti_seqfourdp is buggy/not portable when nx /= ldx .or. ny /= ldy .or. nz /= ldz")
285 : : end if
286 : :
287 : 1162972 : select case (isign)
288 : : case (+1) ! G --> R
289 : 397880 : call dfti_c2r_op(nx,ny,nz,ldx,ldy,ldz,ndat,fofg,fofr)
290 : : case (-1) ! R --> G
291 : 367212 : call dfti_r2c_op(nx,ny,nz,ldx,ldy,ldz,ndat,fofr,fofg)
292 : : case default
293 [ + + - ]: 765092 : ABI_BUG("Wrong isign")
294 : : end select
295 : :
296 : : case default
297 [ + + - ]: 1263123 : ABI_BUG("Wrong value for cplex")
298 : : end select
299 : :
300 : 1263123 : end subroutine dfti_seqfourdp
301 : : !!***
302 : :
303 : : !----------------------------------------------------------------------
304 : :
305 : : !!****f* m_dfti/dfti_seqfourwf
306 : : !! NAME
307 : : !! dfti_seqfourwf
308 : : !!
309 : : !! FUNCTION
310 : : !! Carry out composite Fourier transforms between real and reciprocal (G) space.
311 : : !! Wavefunctions, contained in a sphere in reciprocal space,
312 : : !! can be FFT to real space. They can also be FFT from real space
313 : : !! to a sphere. Also, the density maybe accumulated, and a local potential can be applied.
314 : : !!
315 : : !! The different options are :
316 : : !! - option=0 --> reciprocal to real space and output the result.
317 : : !! - option=1 --> reciprocal to real space and accumulate the density.
318 : : !! - option=2 --> reciprocal to real space, apply the local potential to the wavefunction
319 : : !! in real space and produce the result in reciprocal space.
320 : : !! - option=3 --> real space to reciprocal space.
321 : : !! NOTE that in this case, fftalg=1x1 MUST be used. This may be changed in the future.
322 : : !!
323 : : !! INPUTS
324 : : !! cplex= if 1 , denpot is real, if 2 , denpot is complex
325 : : !! (cplex=2 only allowed for option=2, and istwf_k=1)
326 : : !! not relevant if option=0 or option=3, so cplex=0 can be used to minimize memory
327 : : !! fofgin(2,npwin)=holds input wavefunction in G vector basis sphere.
328 : : !! (intent(in) but the routine sphere can modify it for another iflag)
329 : : !! gboundin(2*mgfft+8,2)=sphere boundary info for reciprocal to real space
330 : : !! gboundout(2*mgfft+8,2)=sphere boundary info for real to reciprocal space
331 : : !! istwf_k=option parameter that describes the storage of wfs
332 : : !! kg_kin(3,npwin)=reduced planewave coordinates, input
333 : : !! kg_kout(3,npwout)=reduced planewave coordinates, output
334 : : !! mgfft=maximum size of 1D FFTs
335 : : !! ndat=number of FFT to do in //
336 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
337 : : !! npwin=number of elements in fofgin array (for option 0, 1 and 2)
338 : : !! npwout=number of elements in fofgout array (for option 2 and 3)
339 : : !! ldx,ldy,ldz=ngfft(4),ngfft(5),ngfft(6), dimensions of fofr.
340 : : !! option= if 0: do direct FFT
341 : : !! if 1: do direct FFT, then sum the density
342 : : !! if 2: do direct FFT, multiply by the potential, then do reverse FFT
343 : : !! if 3: do reverse FFT only
344 : : !! weight_r=weight to be used for the accumulation of the density in real space
345 : : !! (needed only when option=1)
346 : :
347 : : !! OUTPUT
348 : : !! (see side effects)
349 : : !!
350 : : !! SIDE EFFECTS
351 : : !! Input/Output
352 : : !! for option==0, fofgin(2,npwin*ndat)=holds input wavefunction in G sphere;
353 : : !! fofr(2,ldx,ldy,ldz) contains the output Fourier Transform of fofgin;
354 : : !! no use of denpot, fofgout and npwout.
355 : : !! for option==1, fofgin(2,npwin*ndat)=holds input wavefunction in G sphere;
356 : : !! denpot(cplex*ldx,ldy,ldz) contains the input density at input,
357 : : !! and the updated density at output (accumulated);
358 : : !! no use of fofgout and npwout.
359 : : !! for option==2, fofgin(2,npwin*ndat)=holds input wavefunction in G sphere;
360 : : !! denpot(cplex*ldx,ldy,ldz) contains the input local potential;
361 : : !! fofgout(2,npwout*ndat) contains the output function;
362 : : !! for option==3, fofr(2,ldx,ldy,ldz*ndat) contains the input real space wavefunction;
363 : : !! fofgout(2,npwout*ndat) contains its output Fourier transform;
364 : : !! no use of fofgin and npwin.
365 : : !!
366 : : !! SOURCE
367 : :
368 : 59878649 : subroutine dfti_seqfourwf(cplex,denpot,fofgin,fofgout,fofr,gboundin,gboundout,istwf_k, &
369 : 59878649 : kg_kin,kg_kout,mgfft,ndat,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
370 : :
371 : : !Arguments ------------------------------------
372 : : !scalars
373 : : integer,intent(in) :: cplex,istwf_k,ldx,ldy,ldz,ndat,npwin,npwout,option,mgfft
374 : : real(dp),intent(in) :: weight_i,weight_r
375 : : !arrays
376 : : integer,intent(in) :: gboundin(2*mgfft+8,2),gboundout(2*mgfft+8,2)
377 : : integer,intent(in) :: kg_kin(3,npwin),kg_kout(3,npwout),ngfft(18)
378 : : real(dp),intent(inout) :: denpot(cplex*ldx,ldy,ldz),fofgin(2,npwin*ndat)
379 : : real(dp),intent(inout) :: fofr(2,ldx*ldy*ldz*ndat)
380 : : real(dp),intent(out) :: fofgout(2,npwout*ndat)
381 : :
382 : : !Local variables-------------------------------
383 : : !scalars
384 : : integer,parameter :: ndat1=1
385 : : integer :: nx,ny,nz,fftalg,fftalga,fftalgc,fftcache,dat,ptg,ptr,ptgin,ptgout,nthreads
386 : : logical :: use_fftrisc
387 : : character(len=500) :: msg
388 : : ! *************************************************************************
389 : :
390 [ + - + + : 162778840 : if (all(option /= [0, 1, 2, 3])) then
- + ]
391 : 0 : write(msg,'(a,i0,a)')' Option:',option,' is not allowed. Only option=0, 1, 2 or 3 are allowed presently.'
392 : 0 : ABI_ERROR(msg)
393 : : end if
394 : :
395 [ + + - + ]: 59878649 : if (option == 1 .and. cplex /= 1) then
396 : 0 : ABI_ERROR(sjoin("With option number 1, cplex must be 1 but it is cplex:", itoa(cplex)))
397 : : end if
398 : :
399 [ + + - + ]: 59878649 : if (option==2 .and. (cplex/=1 .and. cplex/=2)) then
400 : 0 : ABI_ERROR(sjoin("With the option number 2, cplex must be 1 or 2, but it is cplex:", itoa(cplex)))
401 : : end if
402 : :
403 : 59878649 : nx=ngfft(1); ny=ngfft(2); nz=ngfft(3)
404 : 59878649 : fftalg=ngfft(7); fftalga=fftalg/100; fftalgc=MOD(fftalg,10)
405 : 59878649 : fftcache=ngfft(8)
406 : :
407 : 59878649 : use_fftrisc = (fftalgc==2)
408 [ + + + + ]: 59878649 : if (istwf_k==2.and.option==3) use_fftrisc = .FALSE.
409 [ + + + + : 165061272 : if (istwf_k>2.and.ANY(option==(/0,3/))) use_fftrisc = .FALSE.
+ + + + ]
410 : :
411 : : nthreads = xomp_get_num_threads(open_parallel=.TRUE.)
412 : :
413 [ + + ]: 59878614 : if (use_fftrisc) then
414 : : !call wrtout(std_out, " calls dfti_fftrisc")
415 [ + + ]: 59878598 : if (ndat == 1) then
416 [ + + ]: 59047491 : if (fftcore_mixprec == 0) then
417 : : call dfti_fftrisc_dp(cplex,denpot,fofgin,fofgout,fofr,gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
418 : 58855179 : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
419 : : else
420 : : call dfti_fftrisc_mixprec(cplex,denpot,fofgin,fofgout,fofr,gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
421 : 192312 : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
422 : : end if
423 : : else
424 : : ! All this boilerplate code is needed because the caller might pass zero-sized arrays
425 : : ! for the arguments that are not referenced and we don't want to have problems at run-time.
426 : : ! Moreover option 1 requires a special treatment when threads are started at this level.
427 : :
428 : : SELECT CASE (option)
429 : : CASE (0)
430 : : !
431 : : ! fofgin -> fofr, no use of denpot, fofgout and npwout.
432 : : if (.not.dfti_spawn_threads_here(ndat,nthreads)) then
433 [ + + ]: 1115368 : do dat=1,ndat
434 : 912847 : ptg = 1 + (dat-1)*npwin
435 : 912847 : ptr = 1 + (dat-1)*ldx*ldy*ldz
436 : : call dfti_fftrisc_dp(cplex,denpot,fofgin(1,ptg),fofgout,fofr(1,ptr),gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
437 : 1115368 : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
438 : : end do
439 : : else
440 : : !$OMP PARALLEL DO PRIVATE(ptg,ptr)
441 : : do dat=1,ndat
442 : : ptg = 1 + (dat-1)*npwin
443 : : ptr = 1 + (dat-1)*ldx*ldy*ldz
444 : : call dfti_fftrisc_dp(cplex,denpot,fofgin(1,ptg),fofgout,fofr(1,ptr),gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
445 : : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
446 : : end do
447 : : end if
448 : :
449 : : CASE (1)
450 : : !fofgin -> local ur and accumulate density in denpot
451 : : ! TODO this is delicate part to do in parallel, as one should OMP reduce denpot.
452 [ + + ]: 10 : do dat=1,ndat
453 : 8 : ptg = 1 + (dat-1)*npwin
454 : 8 : ptr = 1 + (dat-1)*ldx*ldy*ldz
455 : : call dfti_fftrisc_dp(cplex,denpot,fofgin(1,ptg),fofgout,fofr,gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
456 : 10 : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
457 : : end do
458 : :
459 : : CASE (2)
460 : : ! <G|vloc(r)|fofgin(r)> in fofgout
461 : : if (.not.dfti_spawn_threads_here(ndat,nthreads)) then
462 [ + + ]: 3921796 : do dat=1,ndat
463 : 3293214 : ptgin = 1 + (dat-1)*npwin
464 : 3293214 : ptgout = 1 + (dat-1)*npwout
465 [ + + ]: 3921796 : if (fftcore_mixprec == 0) then
466 : : call dfti_fftrisc_dp(cplex,denpot,fofgin(1,ptgin),fofgout(1,ptgout),fofr,gboundin,gboundout,istwf_k,&
467 : 3293022 : kg_kin,kg_kout,mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
468 : : else
469 : : call dfti_fftrisc_mixprec(cplex,denpot,fofgin(1,ptgin),fofgout(1,ptgout),fofr,gboundin,gboundout,istwf_k,&
470 : 192 : kg_kin,kg_kout,mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
471 : : end if
472 : : end do
473 : : else
474 : : !$OMP PARALLEL DO PRIVATE(ptgin,ptgout)
475 : : do dat=1,ndat
476 : : ptgin = 1 + (dat-1)*npwin
477 : : ptgout = 1 + (dat-1)*npwout
478 : : call dfti_fftrisc_dp(cplex,denpot,fofgin(1,ptgin),fofgout(1,ptgout),fofr,gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
479 : : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
480 : : end do
481 : : end if
482 : :
483 : : CASE (3)
484 : : ! fofr -> fofgout
485 : 0 : if (.not.dfti_spawn_threads_here(ndat,nthreads)) then
486 [ + + ]: 10 : do dat=1,ndat
487 : 8 : ptr = 1 + (dat-1)*ldx*ldy*ldz
488 : 8 : ptgout = 1 + (dat-1)*npwout
489 : : call dfti_fftrisc_dp(cplex,denpot,fofgin,fofgout(1,ptgout),fofr(1,ptr),gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
490 : 10 : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
491 : : end do
492 : : else
493 : : !$OMP PARALLEL DO PRIVATE(ptr,ptgout)
494 : : do dat=1,ndat
495 : : ptr = 1 + (dat-1)*ldx*ldy*ldz
496 : : ptgout = 1 + (dat-1)*npwout
497 : : call dfti_fftrisc_dp(cplex,denpot,fofgin,fofgout(1,ptgout),fofr(1,ptr),gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
498 : : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,weight_r,weight_i)
499 : : end do
500 : : end if
501 : :
502 : : CASE DEFAULT
503 : 0 : write(msg,'(a,i0,a)')'Option',option,' is not allowed. Only option=0, 1, 2 or 3 are allowed presently.'
504 [ + + + + : 831107 : ABI_ERROR(msg)
- ]
505 : : END SELECT
506 : :
507 : : end if
508 : :
509 : : else
510 : 0 : SELECT CASE (option)
511 : : CASE (0)
512 : : !
513 : : ! FFT u(g) --> u(r)
514 : 0 : if (.not.dfti_spawn_threads_here(ndat,nthreads)) then
515 : 0 : call dfti_fftug_dp(fftalg,fftcache,npwin,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft,kg_kin,gboundin,fofgin,fofr)
516 : : else
517 : : !$OMP PARALLEL DO PRIVATE(ptg, ptr)
518 : : do dat=1,ndat
519 : : ptg = 1 + (dat-1)*npwin
520 : : ptr = 1 + (dat-1)*ldx*ldy*ldz
521 : : call dfti_fftug_dp(fftalg,fftcache,npwin,nx,ny,nz,ldx,ldy,ldz,ndat1,&
522 : : istwf_k,mgfft,kg_kin,gboundin,fofgin(1,ptg),fofr(1,ptr))
523 : : end do
524 : : end if
525 : :
526 : : CASE (1)
527 : : ! TODO this is delicate part to do in parallel, as one should OMP reduce denpot.
528 : 0 : call dfti_fftug_dp(fftalg,fftcache,npwin,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft,kg_kin,gboundin,fofgin,fofr)
529 : 0 : call cg_addtorho(nx,ny,nz,ldx,ldy,ldz,ndat,weight_r,weight_i,fofr,denpot)
530 : :
531 : : CASE (2)
532 : :
533 : 51 : if (.not.dfti_spawn_threads_here(ndat,nthreads)) then
534 : 0 : call dfti_fftug_dp(fftalg,fftcache,npwin,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft,kg_kin,gboundin,fofgin,fofr)
535 : 0 : call cg_vlocpsi(nx,ny,nz,ldx,ldy,ldz,ndat,cplex,denpot,fofr)
536 : :
537 : : ! The data for option==2 is now in fofr.
538 : 0 : call dfti_fftpad_dp(fofr,nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,-1,gboundout)
539 : :
540 : 0 : call cg_box2gsph(nx,ny,nz,ldx,ldy,ldz,ndat,npwout,kg_kout,fofr,fofgout)
541 : : else
542 : :
543 : : !$OMP PARALLEL DO PRIVATE(ptg, ptr)
544 : : do dat=1,ndat
545 : : ptg = 1 + (dat-1)*npwin
546 : : ptr = 1 + (dat-1)*ldx*ldy*ldz
547 : : call dfti_fftug_dp(fftalg,fftcache,npwin,nx,ny,nz,ldx,ldy,ldz,ndat1,&
548 : : istwf_k,mgfft,kg_kin,gboundin,fofgin(1,ptg),fofr(1,ptr))
549 : :
550 : : call cg_vlocpsi(nx,ny,nz,ldx,ldy,ldz,ndat1,cplex,denpot,fofr(1,ptr))
551 : :
552 : : ! The data for option==2 is now in fofr.
553 : : call dfti_fftpad_dp(fofr(1,ptr),nx,ny,nz,ldx,ldy,ldz,ndat1,mgfft,-1,gboundout)
554 : :
555 : : ptg = 1 + (dat-1)*npwout
556 : : call cg_box2gsph(nx,ny,nz,ldx,ldy,ldz,ndat1,npwout,kg_kout,fofr(1,ptr),fofgout(1,ptg))
557 : : end do
558 : : end if
559 : :
560 : : CASE (3)
561 : : ! The data for option==3 is already in fofr.
562 : 0 : if (.not.dfti_spawn_threads_here(ndat,nthreads)) then
563 : 51 : call dfti_fftpad_dp(fofr,nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,-1,gboundout)
564 : 51 : call cg_box2gsph(nx,ny,nz,ldx,ldy,ldz,ndat,npwout,kg_kout,fofr,fofgout)
565 : : else
566 : : !$OMP PARALLEL DO PRIVATE(ptg, ptr)
567 : : do dat=1,ndat
568 : : ptg = 1 + (dat-1)*npwout
569 : : ptr = 1 + (dat-1)*ldx*ldy*ldz
570 : : call dfti_fftpad_dp(fofr(1,ptr),nx,ny,nz,ldx,ldy,ldz,ndat1,mgfft,-1,gboundout)
571 : : call cg_box2gsph(nx,ny,nz,ldx,ldy,ldz,ndat1,npwout,kg_kout,fofr(1,ptr),fofgout(1,ptg))
572 : : end do
573 : : end if
574 : :
575 : : CASE DEFAULT
576 : 0 : write(msg,'(a,i0,a)')'Option',option,' is not allowed. Only option=0, 1, 2 or 3 are allowed presently.'
577 [ - - - + : 51 : ABI_ERROR(msg)
- ]
578 : : END SELECT
579 : : end if
580 : :
581 : 59878649 : end subroutine dfti_seqfourwf
582 : : !!***
583 : :
584 : : !----------------------------------------------------------------------
585 : :
586 : : !!****f* m_dfti/dfti_fftrisc_sp
587 : : !! NAME
588 : : !! dfti_fftrisc_sp
589 : : !!
590 : : !! FUNCTION
591 : : !! Carry out Fourier transforms between real and reciprocal (G) space,
592 : : !! for wavefunctions, contained in a sphere in reciprocal space,
593 : : !! in both directions. Also accomplish some post-processing.
594 : : !! See dfti_fftrisc_dp for API doc.
595 : : !!
596 : : !! SOURCE
597 : :
598 [ + + ]: 15 : subroutine dfti_fftrisc_sp(cplex,denpot,fofgin,fofgout,fofr,gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
599 : : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option, &
600 : : weight_r,weight_i, abi_convention, iscale)
601 : :
602 : : !Arguments ------------------------------------
603 : : !scalars
604 : : integer,intent(in) :: cplex,istwf_k,mgfft,ldx,ldy,ldz,npwin,npwout,option
605 : : real(dp),intent(in) :: weight_i,weight_r
606 : : !arrays
607 : : integer,intent(in) :: gboundin(2*mgfft+8,2),gboundout(2*mgfft+8,2)
608 : : integer,intent(in) :: kg_kin(3,npwin),kg_kout(3,npwout),ngfft(18)
609 : : real(sp),intent(in) :: fofgin(2,npwin)
610 : : real(dp),intent(inout) :: denpot(cplex*ldx,ldy,ldz)
611 : : real(sp),intent(inout) :: fofr(2,ldx*ldy*ldz)
612 : : real(sp),intent(inout) :: fofgout(2,npwout) !vz_i
613 : : logical,optional,intent(in) :: abi_convention
614 : : integer,optional,intent(in) :: iscale
615 : : ! *************************************************************************
616 : :
617 : : #ifdef HAVE_DFTI
618 : :
619 : : #undef FFT_PRECISION
620 : : #undef MYKIND
621 : : #undef MYCZERO
622 : : #undef MYCMPLX
623 : : #undef MYCONJG
624 : :
625 : : #define FFT_PRECISION DFTI_SINGLE
626 : : #define MYKIND SP
627 : : #define MYCZERO (0._sp,0._sp)
628 : : #define MYCMPLX CMPLX
629 : : #define MYCONJG CONJG
630 : :
631 : : #include "dfti_fftrisc.finc"
632 : :
633 : : #else
634 : : ABI_ERROR("DFTI support not activated")
635 : : ABI_UNUSED((/cplex,gboundin(1,1),gboundout(1,1),istwf_k,kg_kin(1,1),kg_kout(1,1),iscale/))
636 : : ABI_UNUSED((/mgfft,ngfft(1),npwin,npwout,ldx,ldy,ldz,option/))
637 : : ABI_UNUSED((/denpot(1,1,1),weight_r,weight_i/))
638 : : ABI_UNUSED((/fofgin(1,1),fofgout(1,1),fofr(1,1)/))
639 : : ABI_UNUSED(abi_convention)
640 : : #endif
641 : :
642 : : end subroutine dfti_fftrisc_sp
643 : : !!***
644 : :
645 : : !----------------------------------------------------------------------
646 : :
647 : : !!****f* m_dfti/dfti_fftrisc_dp
648 : : !! NAME
649 : : !! dfti_fftrisc_dp
650 : : !!
651 : : !! FUNCTION
652 : : !! Carry out Fourier transforms between real and reciprocal (G) space,
653 : : !! for wavefunctions, contained in a sphere in reciprocal space,
654 : : !! in both directions. Also accomplish some post-processing.
655 : : !!
656 : : !! NOTES
657 : : !! Specifically uses rather sophisticated algorithms, based on S Goedecker
658 : : !! routines, specialized for superscalar RISC architecture.
659 : : !! Zero padding : saves 7/12 execution time
660 : : !! Bi-dimensional data locality in most of the routine : cache reuse
661 : : !! For k-point (0 0 0) : takes advantage of symmetry of data.
662 : : !! Note however that no blocking is used, in both 1D z-transform
663 : : !! or subsequent 2D transform. This should be improved.
664 : : !!
665 : : !! INPUTS
666 : : !! cplex= if 1 , denpot is real, if 2 , denpot is complex
667 : : !! (cplex=2 only allowed for option=2 when istwf_k=1)
668 : : !! one can also use cplex=0 if option=0 or option=3
669 : : !! fofgin(2,npwin)=holds input wavefunction in G vector basis sphere.
670 : : !! gboundin(2*mgfft+8,2)=sphere boundary info for reciprocal to real space
671 : : !! gboundout(2*mgfft+8,2)=sphere boundary info for real to reciprocal space
672 : : !! istwf_k=option parameter that describes the storage of wfs
673 : : !! kg_kin(3,npwin)=reduced planewave coordinates, input
674 : : !! kg_kout(3,npwout)=reduced planewave coordinates, output
675 : : !! mgfft=maximum size of 1D FFTs
676 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
677 : : !! npwin=number of elements in fofgin array (for option 0, 1 and 2)
678 : : !! npwout=number of elements in fofgout array (for option 2 and 3)
679 : : !! ldx,ldy,ldz=ngfft(4),ngfft(5),ngfft(6), dimensions of fofr.
680 : : !! option= if 0: do direct FFT
681 : : !! if 1: do direct FFT, then sum the density
682 : : !! if 2: do direct FFT, multiply by the potential, then do reverse FFT
683 : : !! if 3: do reverse FFT only
684 : : !! weight=weight to be used for the accumulation of the density in real space
685 : : !! (needed only when option=1)
686 : : !!
687 : : !! OUTPUT
688 : : !! (see side effects)
689 : : !!
690 : : !! OPTIONS
691 : : !! The different options are:
692 : : !! - reciprocal to real space and output the result (when option=0),
693 : : !! - reciprocal to real space and accumulate the density (when option=1) or
694 : : !! - reciprocal to real space, apply the local potential to the wavefunction
695 : : !! in real space and produce the result in reciprocal space (when option=2)
696 : : !! - real space to reciprocal space (when option=3).
697 : : !! option=0 IS NOT ALLOWED when istwf_k>2
698 : : !! option=3 IS NOT ALLOWED when istwf_k>=2
699 : : !!
700 : : !! SIDE EFFECTS
701 : : !! for option==0, fofgin(2,npwin)=holds input wavefunction in G sphere;
702 : : !! fofr(2,ldx,ldy,ldz) contains the Fourier Transform of fofgin;
703 : : !! no use of denpot, fofgout and npwout.
704 : : !! for option==1, fofgin(2,npwin)=holds input wavefunction in G sphere;
705 : : !! denpot(cplex*ldx,ldy,ldz) contains the input density at input,
706 : : !! and the updated density at output;
707 : : !! no use of fofgout and npwout.
708 : : !! for option==2, fofgin(2,npwin)=holds input wavefunction in G sphere;
709 : : !! denpot(cplex*ldx,ldy,ldz) contains the input local potential;
710 : : !! fofgout(2,npwout) contains the output function;
711 : : !! for option==3, fofr(2,ldx,ldy,ldz) contains the real space wavefunction;
712 : : !! fofgout(2,npwout) contains its Fourier transform;
713 : : !! no use of fofgin and npwin.
714 : : !!
715 : : !! SOURCE
716 : :
717 [ + + ]: 64055066 : subroutine dfti_fftrisc_dp(cplex,denpot,fofgin,fofgout,fofr,gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
718 : : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option, &
719 : : weight_r, weight_i, abi_convention, iscale)
720 : :
721 : : !Arguments ------------------------------------
722 : : !scalars
723 : : integer,intent(in) :: cplex,istwf_k,mgfft,ldx,ldy,ldz,npwin,npwout,option
724 : : real(dp),intent(in) :: weight_i,weight_r
725 : : !arrays
726 : : integer,intent(in) :: gboundin(2*mgfft+8,2),gboundout(2*mgfft+8,2)
727 : : integer,intent(in) :: kg_kin(3,npwin),kg_kout(3,npwout),ngfft(18)
728 : : real(dp),intent(in) :: fofgin(2,npwin)
729 : : real(dp),intent(inout) :: denpot(cplex*ldx,ldy,ldz)
730 : : real(dp),intent(inout) :: fofr(2,ldx*ldy*ldz)
731 : : real(dp),intent(inout) :: fofgout(2,npwout) !vz_i
732 : : logical,optional,intent(in) :: abi_convention
733 : : integer,optional,intent(in) :: iscale
734 : : ! *************************************************************************
735 : :
736 : : #ifdef HAVE_DFTI
737 : :
738 : : #undef FFT_PRECISION
739 : : #undef MYKIND
740 : : #undef MYCZERO
741 : : #undef MYCMPLX
742 : : #undef MYCONJG
743 : :
744 : : #define FFT_PRECISION DFTI_DOUBLE
745 : : #define MYKIND DP
746 : : #define MYCZERO (0._dp,0._dp)
747 : : #define MYCMPLX DCMPLX
748 : : #define MYCONJG DCONJG
749 : :
750 : : #include "dfti_fftrisc.finc"
751 : :
752 : : #else
753 : : ABI_ERROR("DFTI support not activated")
754 : : ABI_UNUSED((/cplex,gboundin(1,1),gboundout(1,1),istwf_k,kg_kin(1,1),kg_kout(1,1),iscale/))
755 : : ABI_UNUSED((/mgfft,ngfft(1),npwin,npwout,ldx,ldy,ldz,option/))
756 : : ABI_UNUSED((/denpot(1,1,1),fofgin(1,1),fofgout(1,1),fofr(1,1),weight_r,weight_i/))
757 : : ABI_UNUSED(abi_convention)
758 : : #endif
759 : :
760 : : end subroutine dfti_fftrisc_dp
761 : : !!***
762 : :
763 : : !----------------------------------------------------------------------
764 : :
765 : : !!****f* m_dfti/dfti_fftrisc_mixprec
766 : : !! NAME
767 : : !! dfti_fftrisc_mixprec
768 : : !!
769 : : !! FUNCTION
770 : : !! Carry out Fourier transforms between real and reciprocal (G) space,
771 : : !! for wavefunctions, contained in a sphere in reciprocal space,
772 : : !! in both directions. Also accomplish some post-processing.
773 : : !! This is the Mixed Precision version (dp in input, FFT done with sp data, output is dp)
774 : : !! See dfti_fftrisc_dp for API doc.
775 : : !!
776 : : !! SOURCE
777 : :
778 [ - + ]: 192504 : subroutine dfti_fftrisc_mixprec(cplex,denpot,fofgin,fofgout,fofr,gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
779 : : mgfft,ngfft,npwin,npwout,ldx,ldy,ldz,option,&
780 : : weight_r,weight_i, abi_convention, iscale) ! optional
781 : :
782 : : !Arguments ------------------------------------
783 : : !scalars
784 : : integer,intent(in) :: cplex,istwf_k,mgfft,ldx,ldy,ldz,npwin,npwout,option
785 : : real(dp),intent(in) :: weight_i,weight_r
786 : : !arrays
787 : : integer,intent(in) :: gboundin(2*mgfft+8,2),gboundout(2*mgfft+8,2)
788 : : integer,intent(in) :: kg_kin(3,npwin),kg_kout(3,npwout),ngfft(18)
789 : : real(dp),intent(in) :: fofgin(2,npwin)
790 : : real(dp),intent(inout) :: denpot(cplex*ldx,ldy,ldz)
791 : : real(dp),intent(inout) :: fofr(2,ldx*ldy*ldz)
792 : : real(dp),intent(inout) :: fofgout(2,npwout) !vz_i
793 : : logical,optional,intent(in) :: abi_convention
794 : : integer,optional,intent(in) :: iscale
795 : : ! *************************************************************************
796 : :
797 : : #ifdef HAVE_DFTI
798 : :
799 : : #undef FFT_PRECISION
800 : : #undef MYKIND
801 : : #undef MYCZERO
802 : : #undef MYCMPLX
803 : : #undef MYCONJG
804 : :
805 : : #define FFT_PRECISION DFTI_SINGLE
806 : : #define MYKIND SP
807 : : #define MYCZERO (0._sp,0._sp)
808 : : #define MYCMPLX CMPLX
809 : : #define MYCONJG CONJG
810 : :
811 : : #define HAVE_DFTI_MIXED_PRECISION 1
812 : :
813 : : #include "dfti_fftrisc.finc"
814 : :
815 : : #undef HAVE_DFTI_MIXED_PRECISION
816 : :
817 : : #else
818 : : ABI_ERROR("DFTI support not activated")
819 : : ABI_UNUSED((/cplex,gboundin(1,1),gboundout(1,1),istwf_k,kg_kin(1,1),kg_kout(1,1),iscale/))
820 : : ABI_UNUSED((/mgfft,ngfft(1),npwin,npwout,ldx,ldy,ldz,option/))
821 : : ABI_UNUSED((/denpot(1,1,1),fofgin(1,1),fofgout(1,1),fofr(1,1),weight_r,weight_i/))
822 : : ABI_UNUSED(abi_convention)
823 : : #endif
824 : :
825 : : end subroutine dfti_fftrisc_mixprec
826 : : !!***
827 : :
828 : : !----------------------------------------------------------------------
829 : :
830 : : !!****f* m_dfti/dfti_fftug_dp
831 : : !! NAME
832 : : !! dfti_fftug_dp
833 : : !!
834 : : !! FUNCTION
835 : : !! Compute ndat zero-padded FFTs from G to R space.
836 : : !! Mainly used for the transform of wavefunctions.
837 : : !! TARGET: dp arrays with real and imaginary part
838 : : !!
839 : : !! INPUTS
840 : : !! fftalg=FFT algorithm (see input variable)
841 : : !! fftcache=size of the cache (kB)
842 : : !! npw_k=number of plane waves for this k-point.
843 : : !! nx,ny,nz=Number of point along the three directions.
844 : : !! ldx,ldy,ldz=Leading dimensions of the array.
845 : : !! ndat=Number of transforms
846 : : !! istwf_k=Option describing the storage of the wavefunction.
847 : : !! mgfft=Max number of FFT divisions (used to dimension gbound)
848 : : !! kg_k(3,npw_k)=G-vectors in reduced coordinates
849 : : !! gbound(2*mgfft+8,2)=Table for padded-FFT. See sphereboundary.
850 : : !! ug(npw_k*ndat)=wavefunctions in reciprocal space.
851 : : !!
852 : : !! OUTPUT
853 : : !! ur(ldx*ldy*ldz*ndat)=wavefunctions in real space.
854 : : !!
855 : : !! SOURCE
856 : :
857 : 0 : subroutine dfti_fftug_dp(fftalg, fftcache, npw_k, nx, ny, nz, ldx, ldy, ldz, ndat, &
858 : 0 : istwf_k, mgfft, kg_k, gbound, ug, ur, &
859 : : isign, iscale) ! optional
860 : :
861 : : !Arguments ------------------------------------
862 : : !scalars
863 : : integer,intent(in) :: fftalg,fftcache
864 : : integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft
865 : : !arrays
866 : : integer,intent(in) :: gbound(2*mgfft+8,2),kg_k(3,npw_k)
867 : : real(dp),target,intent(in) :: ug(2*npw_k*ndat)
868 : : real(dp),target,intent(inout) :: ur(2*ldx*ldy*ldz*ndat)
869 : : integer,optional,intent(in) :: isign, iscale
870 : :
871 : : #ifdef HAVE_DFTI
872 : : !Local variables-------------------------------
873 : : integer,parameter :: dist=2
874 : : integer :: iscale__, isign__
875 : : real(dp) :: fofgout(2,0)
876 : 0 : real(dp),contiguous, pointer :: real_ug(:,:),real_ur(:,:)
877 : : ! *************************************************************************
878 : :
879 [ # # ]: 0 : iscale__ = 0; if (present(iscale)) iscale__ = iscale
880 [ # # ]: 0 : isign__ = +1; if (present(isign)) isign__ = isign
881 : :
882 : : #undef TK_PREF
883 : : #define TK_PREF(name) CONCAT(cg_,name)
884 : : #undef TK_PREF_box2gsph
885 : : #define TK_PREF_box2gsph cg_box2gsph
886 : : #undef TK_PREF_gsph2box
887 : : #define TK_PREF_gsph2box cg_gsph2box
888 : :
889 : : #undef FFT_PRECISION
890 : : #define FFT_PRECISION FFT_DOUBLE
891 : :
892 : : #include "fftug.finc"
893 : :
894 : : #undef FFT_PRECISION
895 : :
896 : : #else
897 : : ! Silence compiler warning
898 : : ABI_ERROR("FFT_DFTI support not activated")
899 : : ABI_UNUSED((/fftalg,fftcache,npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft,kg_k(1,1),gbound(1,1), isign, iscale/))
900 : : ABI_UNUSED((/ug(1),ur(1)/))
901 : : #endif
902 : :
903 : 0 : end subroutine dfti_fftug_dp
904 : : !!***
905 : :
906 : : !----------------------------------------------------------------------
907 : :
908 : : !!****f* m_dfti/dfti_fftug_spc
909 : : !! NAME
910 : : !! dfti_fftug_spc
911 : : !!
912 : : !! FUNCTION
913 : : !! Compute ndat zero-padded FFTs from G- to R-space .
914 : : !! Mainly used for the transform of wavefunctions.
915 : : !! TARGET: spc arrays
916 : : !!
917 : : !! INPUTS
918 : : !! fftalg=FFT algorithm (see input variable)
919 : : !! fftcache=size of the cache (kB)
920 : : !! npw_k=number of plane waves for this k-point.
921 : : !! nx,ny,nz=Number of point along the three directions.
922 : : !! ldx,ldy,ldz=Leading dimensions of the array.
923 : : !! ndat=Number of transforms
924 : : !! istwf_k=Option describing the storage of the wavefunction.
925 : : !! mgfft=Max number of FFT divisions (used to dimension gbound)
926 : : !! kg_k(3,npw_k)=G-vectors in reduced coordinates
927 : : !! gbound(2*mgfft+8,2)=Table for padded-FFT. See sphereboundary.
928 : : !! ug(npw_k*ndat)=wavefunctions in reciprocal space.
929 : : !!
930 : : !! OUTPUT
931 : : !! ur(ldx*ldy*ldz*ndat)=wavefunctions in real space.
932 : : !!
933 : : !! SOURCE
934 : :
935 : 18 : subroutine dfti_fftug_spc(fftalg, fftcache, npw_k, nx, ny, nz, ldx, ldy, ldz, ndat, &
936 : 18 : istwf_k, mgfft, kg_k, gbound, ug, ur, &
937 : : isign, iscale) ! optional
938 : :
939 : : !Arguments ------------------------------------
940 : : !scalars
941 : : integer,intent(in) :: fftalg,fftcache
942 : : integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft
943 : : integer,optional,intent(in) :: isign, iscale
944 : : !arrays
945 : : integer,intent(in) :: gbound(2*mgfft+8,2),kg_k(3,npw_k)
946 : : complex(sp),target,intent(in) :: ug(npw_k*ndat)
947 : : complex(sp),target,intent(inout) :: ur(ldx*ldy*ldz*ndat)
948 : :
949 : : #ifdef HAVE_DFTI
950 : : !Local variables-------------------------------
951 : : !scalars
952 : : integer,parameter :: dist=1
953 : : integer :: iscale__, isign__
954 : : real(sp) :: fofgout(2,0)
955 : 18 : real(sp),contiguous, pointer :: real_ug(:,:),real_ur(:,:)
956 : : ! *************************************************************************
957 : :
958 [ - + ]: 18 : iscale__ = 0; if (present(iscale)) iscale__ = iscale
959 [ - + ]: 18 : isign__ = +1; if (present(isign)) isign__ = isign
960 : :
961 : : #undef TK_PREF
962 : : #define TK_PREF(name) CONCAT(cplx_,name)
963 : : #undef TK_PREF_box2gsph
964 : : #define TK_PREF_box2gsph cplx_box2gsph
965 : : #undef TK_PREF_gsph2box
966 : : #define TK_PREF_gsph2box cplx_gsph2box
967 : :
968 : : #undef FFT_PRECISION
969 : : #define FFT_PRECISION FFT_SINGLE
970 : :
971 : : #include "fftug.finc"
972 : :
973 : : #undef FFT_PRECISION
974 : :
975 : : #else
976 : : ! Silence compiler warning
977 : : ABI_ERROR("FFT_DFTI support not activated")
978 : : ABI_UNUSED((/fftalg,fftcache,npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft,kg_k(1,1),gbound(1,1), iscale, isign/))
979 : : ABI_UNUSED((/ug(1),ur(1)/))
980 : : #endif
981 : :
982 : 18 : end subroutine dfti_fftug_spc
983 : : !!***
984 : :
985 : : !----------------------------------------------------------------------
986 : :
987 : : !!****f* m_dfti/dfti_fftug_dpc
988 : : !! NAME
989 : : !! dfti_fftug_dpc
990 : : !!
991 : : !! FUNCTION
992 : : !! Compute ndat zero-padded FFTs from G ro R.
993 : : !! Mainly used for the transform of wavefunctions.
994 : : !! TARGET: DP arrays
995 : : !!
996 : : !! INPUTS
997 : : !! fftalg=FFT algorithm (see input variable)
998 : : !! fftcache=size of the cache (kB)
999 : : !! npw_k=number of plane waves for this k-point.
1000 : : !! nx,ny,nz=Number of point along the three directions.
1001 : : !! ldx,ldy,ldz=Leading dimensions of the array.
1002 : : !! ndat=Number of transforms
1003 : : !! istwf_k=Option describing the storage of the wavefunction.
1004 : : !! mgfft=Max number of FFT divisions (used to dimension gbound)
1005 : : !! kg_k(3,npw_k)=G-vectors in reduced coordinates
1006 : : !! gbound(2*mgfft+8,2)=Table for padded-FFT. See sphereboundary.
1007 : : !! ug(npw_k*ndat)=wavefunctions in reciprocal space
1008 : : !!
1009 : : !! OUTPUT
1010 : : !! ur(ldx*ldy*ldz*ndat)=wavefunctions in real space.
1011 : : !!
1012 : : !! SOURCE
1013 : :
1014 : 215423 : subroutine dfti_fftug_dpc(fftalg, fftcache, npw_k, nx, ny, nz, ldx, ldy, ldz, ndat, &
1015 : 215423 : istwf_k, mgfft, kg_k, gbound, ug, ur, &
1016 : : isign, iscale) ! optional
1017 : :
1018 : : !Arguments ------------------------------------
1019 : : !scalars
1020 : : integer,intent(in) :: fftalg,fftcache
1021 : : integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft
1022 : : !arrays
1023 : : integer,intent(in) :: gbound(2*mgfft+8,2),kg_k(3,npw_k)
1024 : : complex(dp),target,intent(in) :: ug(npw_k*ndat)
1025 : : complex(dp),target,intent(inout) :: ur(ldx*ldy*ldz*ndat) !vz_i
1026 : : integer,optional,intent(in) :: isign, iscale
1027 : :
1028 : : #ifdef HAVE_DFTI
1029 : : !Local variables-------------------------------
1030 : : !scalars
1031 : : integer,parameter :: dist=1
1032 : : integer :: iscale__, isign__
1033 : : !arrays
1034 : : real(dp) :: fofgout(2,0)
1035 : 215423 : real(dp),contiguous, pointer :: real_ug(:,:),real_ur(:,:)
1036 : : ! *************************************************************************
1037 : :
1038 [ - + ]: 215423 : iscale__ = 0; if (present(iscale)) iscale__ = iscale
1039 [ - + ]: 215423 : isign__ = +1; if (present(isign)) isign__ = isign
1040 : :
1041 : : #undef TK_PREF
1042 : : #define TK_PREF(name) CONCAT(cplx_,name)
1043 : : #undef TK_PREF_box2gsph
1044 : : #define TK_PREF_box2gsph cplx_box2gsph
1045 : : #undef TK_PREF_gsph2box
1046 : : #define TK_PREF_gsph2box cplx_gsph2box
1047 : :
1048 : : #undef FFT_PRECISION
1049 : : #define FFT_PRECISION FFT_DOUBLE
1050 : :
1051 : : #include "fftug.finc"
1052 : :
1053 : : #undef FFT_PRECISION
1054 : :
1055 : : #else
1056 : : ! Silence compiler warning
1057 : : ABI_ERROR("FFT_DFTI support not activated")
1058 : : ABI_UNUSED((/fftalg,fftcache,npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft,kg_k(1,1),gbound(1,1),isign,iscale/))
1059 : : ABI_UNUSED((/ug(1),ur(1)/))
1060 : : #endif
1061 : :
1062 : 215423 : end subroutine dfti_fftug_dpc
1063 : : !!***
1064 : :
1065 : : !----------------------------------------------------------------------
1066 : :
1067 : : !!****f* m_dfti/dfti_fftur_dp
1068 : : !! NAME
1069 : : !! dfti_fftur_dp
1070 : : !!
1071 : : !! FUNCTION
1072 : : !! Compute ndat zero-padded FFTs from R- to G-space .
1073 : : !! Mainly used for the transform of wavefunctions.
1074 : : !! TARGET: dp arrays with real and imaginary part.
1075 : : !!
1076 : : !! INPUTS
1077 : : !! fftalg=FFT algorithm (see input variable)
1078 : : !! fftcache=size of the cache (kB)
1079 : : !! npw_k=number of plane waves for this k-point.
1080 : : !! nx,ny,nz=Number of point along the three directions.
1081 : : !! ldx,ldy,ldz=Leading dimensions of the array.
1082 : : !! ndat=Number of transforms
1083 : : !! istwf_k=Option describing the storage of the wavefunction.
1084 : : !! mgfft=Max number of FFT divisions (used to dimension gbound)
1085 : : !! kg_k(3,npw_k)=G-vectors in reduced coordinates
1086 : : !! gbound(2*mgfft+8,2)=Table for padded-FFT. See sphereboundary.
1087 : : !!
1088 : : !! SIDE EFFECT
1089 : : !! ur(2,ldx*ldy*ldz*ndat)= In input: wavefunctions in real space.
1090 : : !! Destroyed in output. Do not use ur anymore!
1091 : : !! OUTPUT
1092 : : !! ug(2,npw_k*ndat)=wavefunctions in reciprocal space.
1093 : : !!
1094 : : !! SOURCE
1095 : :
1096 : 0 : subroutine dfti_fftur_dp(fftalg, fftcache, npw_k, nx, ny, nz, ldx, ldy, ldz, ndat, &
1097 : 0 : istwf_k, mgfft, kg_k, gbound, ur, ug, &
1098 : : isign, iscale) ! optional
1099 : :
1100 : : !Arguments ------------------------------------
1101 : : !scalars
1102 : : integer,intent(in) :: fftalg,fftcache
1103 : : integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft
1104 : : integer,optional,intent(in) :: isign, iscale
1105 : : !arrays
1106 : : integer,intent(in) :: gbound(2*mgfft+8,2),kg_k(3,npw_k)
1107 : : real(dp),target,intent(inout) :: ur(2*ldx*ldy*ldz*ndat)
1108 : : real(dp),target,intent(inout) :: ug(2*npw_k*ndat) !vz_i
1109 : :
1110 : : #ifdef HAVE_DFTI
1111 : : !Local variables-------------------------------
1112 : : !scalars
1113 : : integer,parameter :: dist=2
1114 : : integer :: iscale__, isign__
1115 : : !arrays
1116 : : real(dp) :: dum_ugin(2,0)
1117 : 0 : real(dp),contiguous, pointer :: real_ug(:,:),real_ur(:,:)
1118 : : ! *************************************************************************
1119 : :
1120 : 0 : iscale__ = 1; if (present(iscale)) iscale__ = iscale
1121 [ # # ]: 0 : isign__ = -1; if (present(isign)) isign__ = isign
1122 : :
1123 : : #undef TK_PREF
1124 : : #define TK_PREF(name) CONCAT(cg_,name)
1125 : : #undef TK_PREF_box2gsph
1126 : : #define TK_PREF_box2gsph cg_box2gsph
1127 : : #undef TK_PREF_gsph2box
1128 : : #define TK_PREF_gsph2box cg_gsph2box
1129 : :
1130 : : #undef FFT_PRECISION
1131 : : #define FFT_PRECISION FFT_DOUBLE
1132 : :
1133 : : #include "fftur.finc"
1134 : :
1135 : : #undef FFT_PRECISION
1136 : :
1137 : : #else
1138 : : ! Silence compiler warning
1139 : : ABI_ERROR("FFT_DFTI support not activated")
1140 : : ABI_UNUSED((/fftalg,fftcache/))
1141 : : ABI_UNUSED((/npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft,kg_k(1,1),gbound(1,1),iscale,isign/))
1142 : : ABI_UNUSED((/ug(1),ur(1)/))
1143 : : #endif
1144 : :
1145 : 0 : end subroutine dfti_fftur_dp
1146 : : !!***
1147 : :
1148 : : !----------------------------------------------------------------------
1149 : :
1150 : : !!****f* m_dfti/dfti_fftur_spc
1151 : : !! NAME
1152 : : !! dfti_fftur_spc
1153 : : !!
1154 : : !! FUNCTION
1155 : : !! Compute ndat zero-padded FFTs from R- to G-space .
1156 : : !! Mainly used for the transform of wavefunctions.
1157 : : !! TARGET: spc arrays
1158 : : !!
1159 : : !! INPUTS
1160 : : !! fftalg=FFT algorithm (see input variable)
1161 : : !! fftcache=size of the cache (kB)
1162 : : !! npw_k=number of plane waves for this k-point.
1163 : : !! nx,ny,nz=Number of point along the three directions.
1164 : : !! ldx,ldy,ldz=Leading dimensions of the array.
1165 : : !! ndat=Number of transforms
1166 : : !! istwf_k=Option describing the storage of the wavefunction.
1167 : : !! mgfft=Max number of FFT divisions (used to dimension gbound)
1168 : : !! kg_k(3,npw_k)=G-vectors in reduced coordinates
1169 : : !! gbound(2*mgfft+8,2)=Table for padded-FFT. See sphereboundary.
1170 : : !!
1171 : : !! SIDE EFFECT
1172 : : !! ur(ldx*ldy*ldz*ndat)= In input: wavefunctions in real space.
1173 : : !! Destroyed in output. Do not use ur anymore!
1174 : : !!
1175 : : !! OUTPUT
1176 : : !! ug(npw_k*ndat)=wavefunctions in reciprocal space.
1177 : : !!
1178 : : !! SOURCE
1179 : :
1180 : 18 : subroutine dfti_fftur_spc(fftalg, fftcache, npw_k, nx, ny, nz, ldx, ldy, ldz, ndat, &
1181 : 18 : istwf_k, mgfft, kg_k, gbound, ur, ug, &
1182 : : isign, iscale) ! optional
1183 : :
1184 : : !Arguments ------------------------------------
1185 : : !scalars
1186 : : integer,intent(in) :: fftalg,fftcache
1187 : : integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft
1188 : : !arrays
1189 : : integer,intent(in) :: gbound(2*mgfft+8,2),kg_k(3,npw_k)
1190 : : complex(sp),target,intent(inout) :: ur(ldx*ldy*ldz*ndat)
1191 : : complex(sp),target,intent(inout) :: ug(npw_k*ndat) !vz_i
1192 : : integer,optional,intent(in) :: isign, iscale
1193 : :
1194 : : #ifdef HAVE_DFTI
1195 : : !Local variables-------------------------------
1196 : : !scalars
1197 : : integer,parameter :: dist=1
1198 : : integer :: iscale__, isign__
1199 : : !arrays
1200 : : real(sp) :: dum_ugin(2,0)
1201 : 18 : real(sp),contiguous, pointer :: real_ug(:,:),real_ur(:,:)
1202 : : ! *************************************************************************
1203 : :
1204 : 18 : iscale__ = 1; if (present(iscale)) iscale__ = iscale
1205 [ + - ]: 18 : isign__ = -1; if (present(isign)) isign__ = isign
1206 : :
1207 : : #undef TK_PREF
1208 : : #define TK_PREF(name) CONCAT(cplx_,name)
1209 : : #undef TK_PREF_box2gsph
1210 : : #define TK_PREF_box2gsph cplx_box2gsph
1211 : : #undef TK_PREF_gsph2box
1212 : : #define TK_PREF_gsph2box cplx_gsph2box
1213 : :
1214 : : #undef FFT_PRECISION
1215 : : #define FFT_PRECISION FFT_SINGLE
1216 : :
1217 : : #include "fftur.finc"
1218 : :
1219 : : #undef FFT_PRECISION
1220 : :
1221 : : #else
1222 : : ! Silence compiler warning
1223 : : ABI_ERROR("FFT_DFTI support not activated")
1224 : : ABI_UNUSED((/fftalg,fftcache/))
1225 : : ABI_UNUSED((/npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft,kg_k(1,1),gbound(1,1),isign,iscale/))
1226 : : ABI_UNUSED((/ug(1),ur(1)/))
1227 : : #endif
1228 : :
1229 : 18 : end subroutine dfti_fftur_spc
1230 : : !!***
1231 : :
1232 : : !----------------------------------------------------------------------
1233 : :
1234 : : !!****f* m_dfti/dfti_fftur_dpc
1235 : : !! NAME
1236 : : !! dfti_fftur_dpc
1237 : : !!
1238 : : !! FUNCTION
1239 : : !! Compute ndat zero-padded FFTs from R ro G.
1240 : : !! Mainly used for the transform of wavefunctions.
1241 : : !! TARGET: DP arrays
1242 : : !!
1243 : : !! INPUTS
1244 : : !! fftalg=FFT algorithm (see input variable)
1245 : : !! fftcache=size of the cache (kB)
1246 : : !! npw_k=number of plane waves for this k-point.
1247 : : !! nx,ny,nz=Number of point along the three directions.
1248 : : !! ldx,ldy,ldz=Leading dimensions of the array.
1249 : : !! ndat=Number of transforms
1250 : : !! istwf_k=Option describing the storage of the wavefunction.
1251 : : !! mgfft=Max number of FFT divisions (used to dimension gbound)
1252 : : !! kg_k(3,npw_k)=G-vectors in reduced coordinates
1253 : : !! gbound(2*mgfft+8,2)=Table for padded-FFT. See sphereboundary.
1254 : : !!
1255 : : !! SIDE EFFECT
1256 : : !! ur(ldx*ldy*ldz*ndat)= In input: wavefunctions in real space.
1257 : : !! Destroyed in output. Do not use ur anymore!
1258 : : !! OUTPUT
1259 : : !! ug(npw_k*ndat)=wavefunctions in reciprocal space
1260 : : !!
1261 : : !! SOURCE
1262 : :
1263 : 779986 : subroutine dfti_fftur_dpc(fftalg, fftcache, npw_k, nx, ny, nz, ldx, ldy, ldz, ndat, &
1264 : 779986 : istwf_k, mgfft, kg_k, gbound, ur, ug, &
1265 : : isign, iscale) ! optional
1266 : :
1267 : : !Arguments ------------------------------------
1268 : : !scalars
1269 : : integer,intent(in) :: fftalg,fftcache
1270 : : integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft
1271 : : integer,optional,intent(in) :: isign, iscale
1272 : : !arrays
1273 : : integer,intent(in) :: gbound(2*mgfft+8,2),kg_k(3,npw_k)
1274 : : complex(dp),target,intent(inout) :: ur(ldx*ldy*ldz*ndat)
1275 : : complex(dp),target,intent(inout) :: ug(npw_k*ndat) !vz_i
1276 : :
1277 : : #ifdef HAVE_DFTI
1278 : : !Local variables-------------------------------
1279 : : !scalars
1280 : : integer,parameter :: dist=1
1281 : : integer :: iscale__, isign__
1282 : : !arrays
1283 : : real(dp) :: dum_ugin(2,0)
1284 : 779986 : real(dp),contiguous, pointer :: real_ug(:,:),real_ur(:,:)
1285 : : ! *************************************************************************
1286 : :
1287 : 779986 : iscale__ = 1; if (present(iscale)) iscale__ = iscale
1288 [ + - ]: 779986 : isign__ = -1; if (present(isign)) isign__ = isign
1289 : :
1290 : : #undef TK_PREF
1291 : : #define TK_PREF(name) CONCAT(cplx_,name)
1292 : : #undef TK_PREF_box2gsph
1293 : : #define TK_PREF_box2gsph cplx_box2gsph
1294 : : #undef TK_PREF_gsph2box
1295 : : #define TK_PREF_gsph2box cplx_gsph2box
1296 : :
1297 : : #undef FFT_PRECISION
1298 : : #define FFT_PRECISION FFT_DOUBLE
1299 : :
1300 : : #include "fftur.finc"
1301 : :
1302 : : #undef FFT_PRECISION
1303 : :
1304 : : #else
1305 : : ! Silence compiler warning
1306 : : ABI_ERROR("FFT_DFTI support not activated")
1307 : : ABI_UNUSED((/fftalg,fftcache/))
1308 : : ABI_UNUSED((/npw_k,nx,ny,nz,ldx,ldy,ldz,ndat,istwf_k,mgfft,kg_k(1,1),gbound(1,1),iscale,isign/))
1309 : : ABI_UNUSED((/ug(1),ur(1)/))
1310 : : #endif
1311 : :
1312 : 779986 : end subroutine dfti_fftur_dpc
1313 : : !!***
1314 : :
1315 : : !----------------------------------------------------------------------
1316 : :
1317 : : !!****f* m_dfti/dfti_c2c_ip_spc
1318 : : !! NAME
1319 : : !! dfti_c2c_ip_spc
1320 : : !!
1321 : : !! FUNCTION
1322 : : !! Driver routine for in-place 3D complex-complex FFT. TARGET: SPC arrays
1323 : : !!
1324 : : !! INPUTS
1325 : : !! nx,ny,nz=Number of points along the three directions.
1326 : : !! ldx,ldy,ldz=Physical dimensions of the array.
1327 : : !! ndat=Number of FFTs to be done.
1328 : : !! isign= +1 : ff(G) => ff(R); -1 : ff(R) => ff(G)
1329 : : !!
1330 : : !! SIDE EFFECTS
1331 : : !! ff(ldx*ldy*ldz*ndat)=
1332 : : !! In input: the complex array to be transformed.
1333 : : !! In output: the Fourier transform in the space specified by isign.
1334 : : !!
1335 : : !! SOURCE
1336 : :
1337 : 33608 : subroutine dfti_c2c_ip_spc(nx, ny, nz, ldx, ldy, ldz, ndat, iscale, isign, ff)
1338 : :
1339 : : !Arguments ------------------------------------
1340 : : !scalars
1341 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat,iscale,isign
1342 : : !arrays
1343 : : complex(sp),intent(inout) :: ff(ldx*ldy*ldz*ndat)
1344 : : ! *************************************************************************
1345 : :
1346 : : ! Include Fortran template
1347 : : #undef DEV_DFTI_PRECISION
1348 : : #define DEV_DFTI_PRECISION DFTI_SINGLE
1349 : :
1350 : : #include "dfti_c2c_ip.finc"
1351 : :
1352 : 33608 : end subroutine dfti_c2c_ip_spc
1353 : : !!***
1354 : :
1355 : : !----------------------------------------------------------------------
1356 : :
1357 : : !!****f* m_dfti/dfti_c2c_ip_dpc
1358 : : !! NAME
1359 : : !! dfti_c2c_ip_dpc
1360 : : !!
1361 : : !! FUNCTION
1362 : : !! Driver routine for in-place 3D complex-complex FFT. TARGET: DP arrays
1363 : : !!
1364 : : !! INPUTS
1365 : : !! nx,ny,nz=Number of points along the three directions.
1366 : : !! ldx,ldy,ldz=Physical dimensions of the array.
1367 : : !! ndat=Number of FFTs to be done.
1368 : : !! iscale=0 if G --> R FFT should not be scaled.
1369 : : !! isign= +1 : ff(G) => ff(R); -1 : ff(R) => ff(G)
1370 : : !!
1371 : : !! SIDE EFFECTS
1372 : : !! ff(ldx*ldy*ldz*ndat)=
1373 : : !! In input: the complex array to be transformed.
1374 : : !! In output: the Fourier transformed in the space specified by isign.
1375 : : !!
1376 : : !! SOURCE
1377 : :
1378 : 12121 : subroutine dfti_c2c_ip_dpc(nx, ny, nz, ldx, ldy, ldz, ndat, iscale, isign, ff)
1379 : :
1380 : : !Arguments ------------------------------------
1381 : : !scalars
1382 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat,iscale,isign
1383 : : !arrays
1384 : : complex(dp),intent(inout) :: ff(ldx*ldy*ldz*ndat)
1385 : : ! *************************************************************************
1386 : :
1387 : : ! Include Fortran template
1388 : : #undef DEV_DFTI_PRECISION
1389 : : #define DEV_DFTI_PRECISION DFTI_DOUBLE
1390 : :
1391 : : #include "dfti_c2c_ip.finc"
1392 : :
1393 : 12121 : end subroutine dfti_c2c_ip_dpc
1394 : : !!***
1395 : :
1396 : : !----------------------------------------------------------------------
1397 : :
1398 : : !!****f* m_dfti/dfti_c2c_op_spc
1399 : : !! NAME
1400 : : !! dfti_c2c_op_spc
1401 : : !!
1402 : : !! FUNCTION
1403 : : !! Driver routine for out-of-place 3D complex-complex FFT of lengths nx, ny, nz.
1404 : : !! TARGET: spc arrays
1405 : : !!
1406 : : !! INPUTS
1407 : : !! nx,ny,nz=Number of points along the three directions.
1408 : : !! ldx,ldy,ldz=Physical dimensions of the array.
1409 : : !! ndat=Number of FFTs to be done.
1410 : : !! iscale=0 if G --> R FFT should not be scaled.
1411 : : !! isign= +1 : ff(G) => gg(R); -1 : ff(R) => gg(G)
1412 : : !! ff(ldx*ldy*ldz*ndat)=The array to be transformed.
1413 : : !!
1414 : : !! OUTPUT
1415 : : !! gg(ldx*ldy*ldz*ndat)=The FFT of ff.
1416 : : !!
1417 : : !! SOURCE
1418 : :
1419 : 24 : subroutine dfti_c2c_op_spc(nx, ny, nz, ldx, ldy, ldz, ndat, iscale, isign, ff, gg)
1420 : :
1421 : : !Arguments ------------------------------------
1422 : : !scalars
1423 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,isign,ndat, iscale
1424 : : !arrays
1425 : : complex(sp),intent(in) :: ff(ldx*ldy*ldz*ndat)
1426 : : complex(sp),intent(out) :: gg(ldx*ldy*ldz*ndat)
1427 : : ! *************************************************************************
1428 : :
1429 : : ! Include Fortran template
1430 : : #undef DEV_DFTI_PRECISION
1431 : : #define DEV_DFTI_PRECISION DFTI_SINGLE
1432 : :
1433 : : #include "dfti_c2c_op.finc"
1434 : :
1435 : 24 : end subroutine dfti_c2c_op_spc
1436 : : !!***
1437 : :
1438 : : !----------------------------------------------------------------------
1439 : :
1440 : : !!****f* m_dfti/dfti_c2c_op_dpc
1441 : : !! NAME
1442 : : !! dfti_c2c_op_dpc
1443 : : !!
1444 : : !! FUNCTION
1445 : : !! Driver routine for out-of-place 3D complex-complex FFT of lengths nx, ny, nz.
1446 : : !! TARGET: DP arrays
1447 : : !!
1448 : : !! INPUTS
1449 : : !! nx,ny,nz=Number of points along the three directions.
1450 : : !! ldx,ldy,ldz=Physical dimensions of the array.
1451 : : !! ndat=Number of FFTs to be done.
1452 : : !! iscale=0 if G --> R FFT should not be scaled.
1453 : : !! isign= +1 : ff(G) => gg(R); -1 : ff(R) => gg(G)
1454 : : !! ff(ldx*ldy*ldz*ndat)=The array to be transformed.
1455 : : !!
1456 : : !! OUTPUT
1457 : : !! gg(ldx*ldy*ldz*ndat)=The FFT of ff.
1458 : : !!
1459 : : !! SOURCE
1460 : :
1461 : 467313 : subroutine dfti_c2c_op_dpc(nx, ny, nz, ldx, ldy, ldz, ndat, iscale, isign, ff, gg)
1462 : :
1463 : : !Arguments ------------------------------------
1464 : : !scalars
1465 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,isign,ndat,iscale
1466 : : !arrays
1467 : : complex(dp),intent(in) :: ff(ldx*ldy*ldz*ndat)
1468 : : complex(dp),intent(out) :: gg(ldx*ldy*ldz*ndat)
1469 : : ! *************************************************************************
1470 : :
1471 : : ! Include Fortran template
1472 : : #undef DEV_DFTI_PRECISION
1473 : : #define DEV_DFTI_PRECISION DFTI_DOUBLE
1474 : :
1475 : : #include "dfti_c2c_op.finc"
1476 : :
1477 : 467313 : end subroutine dfti_c2c_op_dpc
1478 : : !!***
1479 : :
1480 : : !----------------------------------------------------------------------
1481 : :
1482 : : !!****f* m_dfti/dfti_many_dft_op
1483 : : !! NAME
1484 : : !! dfti_many_dft_op
1485 : : !!
1486 : : !! FUNCTION
1487 : : !! Driver routine for many out-of-place 3D complex-to-complex FFTs of lengths nx, ny, nz.
1488 : : !!
1489 : : !! INPUTS
1490 : : !! nx,ny,nz=Number of points along the three directions.
1491 : : !! ldx,ldy,ldz=Physical dimension of the fin and fout arrays (to avoid cache conflicts).
1492 : : !! ndat=Number of FFTs to be done.
1493 : : !! fin(2*ldx*ldy*ldz*ndat)=The complex array to be transformed.
1494 : : !! isign=sign of Fourier transform exponent: current convention uses
1495 : : !! +1 for transforming from G to r,
1496 : : !! -1 for transforming from r to G.
1497 : : !!
1498 : : !! OUTPUT
1499 : : !! fout(2,ldx*ldy*ldz*ndat)=The Fourier transform of fin.
1500 : : !!
1501 : : !! SOURCE
1502 : :
1503 : 464447 : subroutine dfti_many_dft_op(nx,ny,nz,ldx,ldy,ldz,ndat,isign,fin,fout)
1504 : :
1505 : : !Arguments ------------------------------------
1506 : : !scalars
1507 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat,isign
1508 : : !arrays
1509 : : real(dp),target,intent(in) :: fin(2*ldx*ldy*ldz*ndat)
1510 : : real(dp),target,intent(out) :: fout(2*ldx*ldy*ldz*ndat)
1511 : :
1512 : : #ifdef HAVE_DFTI
1513 : : !Local variables-------------------------------
1514 : : !scalars
1515 : : integer,parameter :: iscale1 = 1
1516 : : type(C_ptr) :: fin_cptr, fout_cptr
1517 : :
1518 : : !arrays
1519 : 464447 : complex(dp),contiguous, pointer :: fin_fptr(:),fout_fptr(:)
1520 : : ! *************************************************************************
1521 : :
1522 : : ! Associate complex pointers with real inputs via the C pointers
1523 : 464447 : fin_cptr = C_loc(fin)
1524 : 464447 : call C_F_pointer(fin_cptr,fin_fptr, shape=[ldx*ldy*ldz*ndat])
1525 : :
1526 : 464447 : fout_cptr = C_loc(fout)
1527 : 464447 : call C_F_pointer(fout_cptr,fout_fptr, shape=[ldx*ldy*ldz*ndat])
1528 : :
1529 : : ! Call complex version --> a lot of boilerplate code avoided
1530 : 464447 : call dfti_c2c_op(nx, ny, nz, ldx, ldy, ldz, ndat, iscale1, isign, fin_fptr, fout_fptr)
1531 : :
1532 : : #else
1533 : : ABI_ERROR("FFT_DFTI support not activated")
1534 : : ABI_UNUSED((/nx,ny,nz,ldx,ldy,ldz,ndat,isign/))
1535 : : ABI_UNUSED(fin(1))
1536 : : ABI_UNUSED(fout(1))
1537 : : #endif
1538 : :
1539 : 464447 : end subroutine dfti_many_dft_op
1540 : : !!***
1541 : :
1542 : : !----------------------------------------------------------------------
1543 : :
1544 : : !!****f* m_dfti/dfti_many_dft_ip
1545 : : !! NAME
1546 : : !! dfti_many_dft_ip
1547 : : !!
1548 : : !! FUNCTION
1549 : : !! Driver routine for many in-place 3D complex-to-complex FFTs of lengths nx, ny, nz.
1550 : : !!
1551 : : !! INPUTS
1552 : : !! nx,ny,nz=Number of points along the three directions.
1553 : : !! ldx,ldy,ldz=Physical dimension of the finout array (to avoid cache conflicts).
1554 : : !! ndat=Number of FFTs to be done.
1555 : : !! isign=sign of Fourier transform exponent: current convention uses
1556 : : !! +1 for transforming from G to r,
1557 : : !! -1 for transforming from r to G.
1558 : : !!
1559 : : !! OUTPUT
1560 : : !! finout(2,ldx*ldy*ldz*ndat)=
1561 : : !! In input: The complex array to be transformed.
1562 : : !! In output: The FFT results.
1563 : : !!
1564 : : !! SOURCE
1565 : :
1566 : 0 : subroutine dfti_many_dft_ip(nx,ny,nz,ldx,ldy,ldz,ndat,isign,finout)
1567 : :
1568 : : !Arguments ------------------------------------
1569 : : !scalars
1570 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat,isign
1571 : : !arrays
1572 : : real(dp),target,intent(inout) :: finout(2*ldx*ldy*ldz*ndat)
1573 : :
1574 : : #ifdef HAVE_DFTI
1575 : : !Local variables-------------------------------
1576 : : !scalars
1577 : : integer,parameter :: iscale1 = 1
1578 : : type(C_ptr) :: finout_cptr
1579 : : !arrays
1580 : : complex(dp),contiguous, pointer :: finout_fptr(:)
1581 : : ! *************************************************************************
1582 : :
1583 : : ! Associate complex finout_fptr with real ffinout via the C pointer
1584 : 0 : finout_cptr = C_loc(finout)
1585 : 0 : call C_F_pointer(finout_cptr,finout_fptr, shape=(/ldx*ldy*ldz*ndat/))
1586 : :
1587 : : ! Call complex version --> a lot of boilerplate code avoided
1588 : 0 : call dfti_c2c_ip(nx, ny, nz, ldx, ldy, ldz, ndat, iscale1, isign, finout_fptr)
1589 : :
1590 : : #else
1591 : : ABI_ERROR("FFT_DFTI support not activated")
1592 : : ABI_UNUSED((/nx,ny,nz,ldx,ldy,ldz,ndat,isign/))
1593 : : ABI_UNUSED(finout(1))
1594 : : #endif
1595 : :
1596 : 0 : end subroutine dfti_many_dft_ip
1597 : : !!***
1598 : :
1599 : : !----------------------------------------------------------------------
1600 : :
1601 : : !!****f* m_dfti/dfti_fftpad_dp
1602 : : !! NAME
1603 : : !! dfti_fftpad_dp
1604 : : !!
1605 : : !! FUNCTION
1606 : : !! This routine transforms wavefunctions using 3D zero-padded FFTs with DFTI.
1607 : : !! The 3D ffts are computed only on lines and planes which have non zero elements (see zpad_init)
1608 : : !! FFT transform is in-place.
1609 : : !!
1610 : : !! INPUTS
1611 : : !! nx,ny,nz=Logical dimensions of the FFT mesh.
1612 : : !! ldx,ldy,ldz=Physical dimension of the f array (to avoid cache conflicts).
1613 : : !! ndat=Number of FFTs
1614 : : !! mgfft=MAX(nx,ny,nz), only used to dimension gbound
1615 : : !! isign=The sign of the transform.
1616 : : !! gbound(2*mgfft+8,2)= The boundaries of the basis sphere of G vectors at a given k-point.
1617 : : !! See sphereboundary for more info.
1618 : : !!
1619 : : !! SIDE EFFECTS
1620 : : !! ff(2*ldx*ldy*ldz*ndat)=
1621 : : !! input: The array with the data to be transformed.
1622 : : !! output: The results of the FFT.
1623 : : !!
1624 : : !! SOURCE
1625 : :
1626 : 51 : subroutine dfti_fftpad_dp(ff, nx, ny, nz, ldx, ldy, ldz, ndat, mgfft, isign, gbound, iscale)
1627 : :
1628 : : !Arguments ------------------------------------
1629 : : !scalars
1630 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,isign
1631 : : !arrays
1632 : : integer,intent(in) :: gbound(2*mgfft+8,2)
1633 : : real(dp),target,intent(inout) :: ff(2*ldx*ldy*ldz*ndat)
1634 : : integer,optional,intent(in) :: iscale
1635 : :
1636 : : !Local variables-------------------------------
1637 : : #ifdef HAVE_DFTI
1638 : : !scalars
1639 : : type(C_ptr) :: cptr
1640 : : integer :: iscale__
1641 : : !arrays
1642 : : complex(dp),contiguous, pointer :: fptr(:)
1643 : : ! *************************************************************************
1644 : :
1645 : 51 : iscale__ = merge(1, 0, isign == -1); if (present(iscale)) iscale__ = iscale
1646 : :
1647 : : ! Associate complex fptr with real ff via the C pointer
1648 : 51 : cptr = C_loc(ff)
1649 : 51 : call C_F_pointer(cptr, fptr, SHAPE=[ldx*ldy*ldz*ndat])
1650 : :
1651 : : ! Call complex version --> a lot of boilerplate code avoided
1652 : 51 : call dfti_fftpad_dpc(fptr, nx, ny, nz, ldx, ldy, ldz, ndat, mgfft, isign, gbound)
1653 : :
1654 : : #else
1655 : : ABI_ERROR("FFT_DFTI support not activated")
1656 : : ABI_UNUSED((/nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,isign,iscale/))
1657 : : ABI_UNUSED(gbound(1,1))
1658 : : ABI_UNUSED(ff(1))
1659 : : #endif
1660 : :
1661 : 51 : end subroutine dfti_fftpad_dp
1662 : : !!***
1663 : :
1664 : : !----------------------------------------------------------------------
1665 : :
1666 : : !!****f* m_dfti/dfti_fftpad_dpc
1667 : : !! NAME
1668 : : !! dfti_fftpad_dpc
1669 : : !!
1670 : : !! FUNCTION
1671 : : !! This routine transforms wavefunctions using 3D zero-padded FFTs with DFTI.
1672 : : !! The 3D ffts are computed only on lines and planes which have non zero elements (see zpad_init)
1673 : : !! FFT transform is in-place. Target: complex double-precision arrays.
1674 : : !!
1675 : : !! INPUTS
1676 : : !! nx,ny,nz=Logical dimensions of the FFT mesh.
1677 : : !! ldx,ldy,ldz=Physical dimension of the f array (to avoid cache conflicts).
1678 : : !! ndat=Number of FFTs.
1679 : : !! mgfft=MAX(nx,ny,nz), only used to dimension gbound
1680 : : !! isign=The sign of the transform.
1681 : : !! gbound(2*mgfft+8,2)= The boundaries of the basis sphere of G vectors at a given k-point.
1682 : : !! See sphereboundary for more info.
1683 : : !!
1684 : : !! SIDE EFFECTS
1685 : : !! ff(ldx*ldy*ldz*ndat)=
1686 : : !! input: The array with the data to be transformed.
1687 : : !! output: The results of the FFT.
1688 : : !!
1689 : : !! SOURCE
1690 : :
1691 [ + + ]: 5739087 : subroutine dfti_fftpad_dpc(ff, nx, ny, nz, ldx, ldy, ldz, ndat, mgfft, isign, gbound, iscale)
1692 : :
1693 : : !Arguments ------------------------------------
1694 : : !scalars
1695 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,isign
1696 : : integer,optional,intent(in) :: iscale
1697 : : !arrays
1698 : : integer,intent(in) :: gbound(2*mgfft+8,2)
1699 : : complex(dp),intent(inout) :: ff(ldx*ldy*ldz*ndat)
1700 : :
1701 : : !Local variables-------------------------------
1702 : : #ifdef HAVE_DFTI
1703 : : ! *************************************************************************
1704 : :
1705 : : ! Include Fortran template
1706 : : #undef DEV_DFTI_PRECISION
1707 : : #define DEV_DFTI_PRECISION DFTI_DOUBLE
1708 : :
1709 : : #include "dfti_fftpad.finc"
1710 : :
1711 : : #else
1712 : : ABI_ERROR("FFT_DFTI support not activated")
1713 : : ABI_UNUSED((/nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,isign,iscale/))
1714 : : ABI_UNUSED(gbound(1,1))
1715 : : ABI_UNUSED(ff(1))
1716 : : #endif
1717 : :
1718 [ - + - + ]: 5739087 : end subroutine dfti_fftpad_dpc
1719 : : !!***
1720 : :
1721 : : !----------------------------------------------------------------------
1722 : :
1723 : : !!****f* m_dfti/dfti_fftpad_spc
1724 : : !! NAME
1725 : : !! dfti_fftpad_spc
1726 : : !!
1727 : : !! FUNCTION
1728 : : !! This routine transforms wavefunctions using 3D zero-padded FFTs with DFTI.
1729 : : !! The 3D ffts are computed only on lines and planes which have non zero elements (see zpad_init)
1730 : : !! FFT transform is in-place. Target: complex SPC arrays.
1731 : : !!
1732 : : !! INPUTS
1733 : : !! nx,ny,nz=Logical dimensions of the FFT mesh.
1734 : : !! ldx,ldy,ldz=Physical dimension of the f array (to avoid cache conflicts).
1735 : : !! ndat=Number of FFTs.
1736 : : !! mgfft=MAX(nx,ny,nz), only used to dimension gbound
1737 : : !! isign=The sign of the transform.
1738 : : !! gbound(2*mgfft+8,2)= The boundaries of the basis sphere of G vectors at a given k-point.
1739 : : !! See sphereboundary for more info.
1740 : : !!
1741 : : !! SIDE EFFECTS
1742 : : !! ff(ldx*ldy*ldz*ndat)=
1743 : : !! input: The array with the data to be transformed.
1744 : : !! output: The results of the FFT.
1745 : : !!
1746 : : !! SOURCE
1747 : :
1748 [ + + ]: 30 : subroutine dfti_fftpad_spc(ff, nx, ny, nz, ldx, ldy, ldz, ndat, mgfft, isign, gbound, iscale)
1749 : :
1750 : : !Arguments ------------------------------------
1751 : : !scalars
1752 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,isign
1753 : : integer,optional,intent(in) :: iscale
1754 : : !arrays
1755 : : integer,intent(in) :: gbound(2*mgfft+8,2)
1756 : : complex(sp),intent(inout) :: ff(ldx*ldy*ldz*ndat)
1757 : :
1758 : : #ifdef HAVE_DFTI
1759 : :
1760 : : ! Include Fortran template
1761 : : #undef DEV_DFTI_PRECISION
1762 : : #define DEV_DFTI_PRECISION DFTI_SINGLE
1763 : :
1764 : : #include "dfti_fftpad.finc"
1765 : :
1766 : : #else
1767 : : ABI_ERROR("FFT_DFTI support not activated")
1768 : : ABI_UNUSED((/nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,isign,iscale/))
1769 : : ABI_UNUSED(gbound(1,1))
1770 : : ABI_UNUSED(ff(1))
1771 : : #endif
1772 : :
1773 [ - + - + ]: 30 : end subroutine dfti_fftpad_spc
1774 : : !!***
1775 : :
1776 : : !----------------------------------------------------------------------
1777 : :
1778 : : !!****f* m_dfti/dfti_r2c_op_dpc
1779 : : !! NAME
1780 : : !! dfti_r2c_op_dpc
1781 : : !!
1782 : : !! FUNCTION
1783 : : !! Driver routine for out-of-place 3D real-to-complex FFT of lengths nx, ny, nz.
1784 : : !!
1785 : : !! INPUTS
1786 : : !! nx,ny,nz=Number of points along the three directions.
1787 : : !! ldx,ldy,ldz=Physical dimensions of the f array (to avoid cache conflicts).
1788 : : !! ff(ldx*ldy*ldz*ndat)=The real array to be transformed.
1789 : : !! ndat=Number of FFTs to be done.
1790 : : !!
1791 : : !! OUTPUT
1792 : : !! gg(ldx*ldy*ldz*ndat)=The forward FFT of ff (complex valued)
1793 : : !!
1794 : : !! SOURCE
1795 : :
1796 : 367212 : subroutine dfti_r2c_op_dpc(nx, ny, nz, ldx, ldy, ldz, ndat, ff, gg)
1797 : :
1798 : : !Arguments ------------------------------------
1799 : : !scalars
1800 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
1801 : : !arrays
1802 : : real(dp),intent(in) :: ff(ldx*ldy*ldz*ndat)
1803 : : complex(dp),intent(out) :: gg(ldx*ldy*ldz*ndat)
1804 : :
1805 : : #ifdef HAVE_DFTI
1806 : : !Local variables-------------------------------
1807 : : !scalars
1808 : : integer :: status,nhp,padx,i1,i2,i3,igp,igf,imgf,ii,arr3(3)
1809 : : integer :: i1inv,i2inv,i3inv,idat,padatf,arr(4)
1810 : : type(DFTI_DESCRIPTOR),pointer :: Desc
1811 : : type(C_PTR) :: cptr
1812 : : !arrays
1813 : 367212 : integer,allocatable :: i1inver(:),i2inver(:),i3inver(:)
1814 : : complex(dp),contiguous, pointer :: gg_hp(:)
1815 : : ! *************************************************************************
1816 : :
1817 : 367212 : padx = (nx/2+1)
1818 : 367212 : nhp = (nx/2+1)*ny*nz
1819 : :
1820 : 367212 : call dfti_alloc_complex(nhp*ndat,cptr,gg_hp)
1821 : :
1822 : 367212 : arr3(1) = nx
1823 : 367212 : arr3(2) = ny
1824 : 367212 : arr3(3) = nz
1825 : 367212 : status = DftiCreateDescriptor(Desc, DFTI_DOUBLE, DFTI_REAL, 3, arr3 )
1826 [ - + ]: 367212 : DFTI_CHECK(status)
1827 : :
1828 : 367212 : status = DftiSetValue(Desc, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX)
1829 : 367212 : status = DftiSetValue(Desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE )
1830 : 367212 : arr(1) = 0
1831 : 367212 : arr(2) = 1
1832 : 367212 : arr(3) = ldx
1833 : 367212 : arr(4) = ldx*ldy
1834 : 367212 : status = DftiSetValue(Desc, DFTI_INPUT_STRIDES, arr)
1835 : 367212 : status = DftiSetValue(Desc, DFTI_INPUT_DISTANCE, ldx*ldy*ldz)
1836 : : arr(1) = 0
1837 : : arr(2) = 1
1838 : 367212 : arr(3) = padx
1839 : 367212 : arr(4) = padx*ny
1840 : 367212 : status = DftiSetValue(Desc, DFTI_OUTPUT_STRIDES, arr)
1841 : 367212 : status = DftiSetValue(Desc, DFTI_OUTPUT_DISTANCE, nhp)
1842 : 367212 : status = DftiSetValue(Desc, DFTI_NUMBER_OF_TRANSFORMS, ndat)
1843 : 367212 : status = DftiSetValue(Desc, DFTI_FORWARD_SCALE, one / DBLE(nx*ny*nz) )
1844 [ - + ]: 367212 : DFTI_CHECK(status)
1845 : :
1846 [ - + ]: 367212 : DFTI_CHECK( DftiCommitDescriptor(Desc) )
1847 : :
1848 [ - + ]: 367212 : DFTI_CHECK( DftiComputeForward(Desc, ff, gg_hp) )
1849 : :
1850 [ - + ]: 367212 : DFTI_CHECK( DftiFreeDescriptor(Desc) )
1851 : :
1852 : : ! Reconstruct full FFT: Hermitian redundancy: out[i] is the conjugate of out[n-i]
1853 [ - + + - : 1101636 : ABI_MALLOC(i1inver,(padx))
+ - + - -
+ - + ]
1854 [ - + + - : 1101636 : ABI_MALLOC(i2inver,(ny))
+ - + - -
+ - + ]
1855 [ - + + - : 1101636 : ABI_MALLOC(i3inver,(nz))
+ - + - -
+ - + ]
1856 : :
1857 : 367212 : i1inver(1)=1
1858 [ + + ]: 4544281 : do i1=2,padx
1859 : 4544281 : i1inver(i1)=nx+2-i1
1860 : : end do
1861 : :
1862 : 367212 : i2inver(1)=1
1863 [ + + ]: 8038833 : do i2=2,ny
1864 : 8038833 : i2inver(i2)=ny+2-i2
1865 : : end do
1866 : :
1867 : 367212 : i3inver(1)=1
1868 [ + + ]: 9187351 : do i3=2,nz
1869 : 9187351 : i3inver(i3)=nz+2-i3
1870 : : end do
1871 : :
1872 : : igp=0
1873 [ + + ]: 734442 : do idat=1,ndat
1874 : 367230 : padatf=(idat-1)*ldx*ldy*ldz
1875 [ + + ]: 9922063 : do i3=1,nz
1876 : 9187621 : i3inv = i3inver(i3)
1877 [ + + ]: 249317781 : do i2=1,ny
1878 : 239762930 : i2inv = i2inver(i2)
1879 [ + + ]: 5129917800 : do i1=1,padx
1880 : 4880967249 : igp=igp+1
1881 : 4880967249 : igf = i1 + (i3-1)*ldx*ldy + (i2-1)*ldx + padatf
1882 : 4880967249 : gg(igf) = gg_hp(igp)
1883 : 4880967249 : i1inv = i1inver(i1)
1884 [ + + ]: 5120730179 : if (i1inv/=i1) then
1885 : 4426842889 : imgf = i1inv + (i3inv-1)*ldx*ldy + (i2inv-1)*ldx + padatf
1886 : 4426842889 : gg(imgf) = DCONJG(gg_hp(igp))
1887 : : end if
1888 : : end do
1889 : : end do
1890 : : end do
1891 : : end do
1892 : :
1893 : 367212 : ABI_FREE(i1inver)
1894 : 367212 : ABI_FREE(i2inver)
1895 : 367212 : ABI_FREE(i3inver)
1896 : :
1897 : 367212 : call dfti_free(cptr)
1898 : :
1899 : : #else
1900 : : ABI_ERROR("FFT_DFTI support not activated")
1901 : : ABI_UNUSED((/nx,ny,nz,ldx,ldy,ldz/))
1902 : : ABI_UNUSED(ff)
1903 : : ABI_UNUSED(gg(1))
1904 : : #endif
1905 : :
1906 : 367212 : end subroutine dfti_r2c_op_dpc
1907 : : !!***
1908 : :
1909 : : !----------------------------------------------------------------------
1910 : :
1911 : : !!****f* m_dfti/dfti_r2c_op_dp
1912 : : !! NAME
1913 : : !! dfti_r2c_op_dp
1914 : : !!
1915 : : !! FUNCTION
1916 : : !! Driver routine for out-of-place 3D real-to-complex FFT of lengths nx, ny, nz.
1917 : : !!
1918 : : !! INPUTS
1919 : : !! nx,ny,nz=Number of points along the three directions.
1920 : : !! ldx,ldy,ldz=Physical dimensions of the f array (to avoid cache conflicts).
1921 : : !! ndat=Number of FFTs to be done.
1922 : : !! ff(ldx*ldy*ldz*ndat)=The real array to be transformed.
1923 : : !!
1924 : : !! OUTPUT
1925 : : !! gg(2*ldx*ldy*ldz*ndat)=The forward FFT of ff (real valued)
1926 : : !!
1927 : : !! SOURCE
1928 : :
1929 : 367212 : subroutine dfti_r2c_op_dp(nx, ny, nz, ldx, ldy, ldz, ndat, ff, gg)
1930 : :
1931 : : !Arguments ------------------------------------
1932 : : !scalars
1933 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
1934 : : !arrays
1935 : : real(dp),intent(in) :: ff(ldx*ldy*ldz*ndat)
1936 : : real(dp),target,intent(out) :: gg(2*ldx*ldy*ldz*ndat)
1937 : :
1938 : : #ifdef HAVE_DFTI
1939 : : !Local variables-------------------------------
1940 : : !scalars
1941 : : type(C_ptr) :: gg_cptr
1942 : : !arrays
1943 : : complex(dp),contiguous, pointer :: gg_fptr(:)
1944 : : ! *************************************************************************
1945 : :
1946 : 367212 : gg_cptr = C_loc(gg)
1947 : 367212 : call C_F_pointer(gg_cptr,gg_fptr, shape=(/ldx*ldy*ldz*ndat/))
1948 : :
1949 : 367212 : call dfti_r2c_op_dpc(nx,ny,nz,ldx,ldy,ldz,ndat,ff,gg_fptr)
1950 : :
1951 : : #else
1952 : : ABI_ERROR("FFT_DFTI support not activated")
1953 : : ABI_UNUSED((/nx,ny,nz,ldx,ldy,ldz/))
1954 : : ABI_UNUSED(ff)
1955 : : ABI_UNUSED(gg(1))
1956 : : #endif
1957 : :
1958 : 367212 : end subroutine dfti_r2c_op_dp
1959 : : !!***
1960 : :
1961 : : !----------------------------------------------------------------------
1962 : :
1963 : : !!****f* m_dfti/dfti_c2r_op_dpc
1964 : : !! NAME
1965 : : !! dfti_c2r_op_dpc
1966 : : !!
1967 : : !! FUNCTION
1968 : : !! Driver routine for out-of-place 3D complex-to-real FFT of lengths nx, ny, nz.
1969 : : !!
1970 : : !! INPUTS
1971 : : !! nx,ny,nz=Number of point along the three directions.
1972 : : !! ldx,ldy,ldz=Physical dimension of the f array (to avoid cache conflicts).
1973 : : !! ndat=Number of FFTs to be done.
1974 : : !! ff(2,ldx*ldy*ldz*ndat)=The complex array to be transformed.
1975 : : !!
1976 : : !! OUTPUT
1977 : : !! gg(2,ldx*ldy*ldz*ndat)=The backwards real FFT of ff.
1978 : : !!
1979 : : !! SOURCE
1980 : :
1981 : 397880 : subroutine dfti_c2r_op_dpc(nx, ny, nz, ldx, ldy, ldz, ndat, ff, gg)
1982 : :
1983 : : !Arguments ------------------------------------
1984 : : !scalars
1985 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
1986 : : !arrays
1987 : : complex(dp),intent(in) :: ff(ldx*ldy*ldz*ndat)
1988 : : real(dp),intent(out) :: gg(ldx*ldy*ldz*ndat)
1989 : :
1990 : : #ifdef HAVE_DFTI
1991 : : !Local variables-------------------------------
1992 : : !scalars
1993 : : integer :: status,nhp,padx,i2,i3,igp,igf,idat,padatf,padatp,ii,arr(4),arr3(3)
1994 : : type(DFTI_DESCRIPTOR),pointer :: Desc
1995 : : type(C_PTR) :: cptr
1996 : : !arrays
1997 : : complex(dp),contiguous, pointer :: ff_hp(:)
1998 : : ! *************************************************************************
1999 : :
2000 : : !stride = 1
2001 : : !idist = nhp
2002 : : !odist = nx*ny*nz
2003 : : !n = (/nx,ny,nz/)
2004 : : !inembed = (/(nx/2+1),ny,nz/)
2005 : : !onembed = (/nx,ny,nz/) ! check this
2006 : : !my_plan = retrieve_plan3(n,ndat,inembed,stride,idist,onembed,stride,odist,FFTW_BACKWARD,my_flags,Saved_plans)
2007 : :
2008 : : ! Fill the Hermitian part: Hermitian redundancy: out[i] is the conjugate of out[n-i]
2009 : 397880 : padx = (nx/2+1)
2010 : 397880 : nhp = padx*ny*nz
2011 : :
2012 : 397880 : call dfti_alloc_complex(nhp*ndat,cptr,ff_hp)
2013 : :
2014 : : !!$OMP PARALLEL DO PRIVATE(padatf,padatp,igf,igp)
2015 [ + + ]: 795778 : do idat=1,ndat
2016 : 397898 : padatf=(idat-1)*ldx*ldy*ldz
2017 : 397898 : padatp=(idat-1)*padx*ny*nz
2018 [ + + ]: 10638880 : do i3=1,nz
2019 [ + + ]: 262980735 : do i2=1,ny
2020 : 252739735 : igf = (i3-1)*ldx*ldy + (i2-1)*ldx + padatf
2021 : 252739735 : igp = (i3-1)*padx*ny + (i2-1)*padx + padatp
2022 [ + + ]: 5214745029 : ff_hp(igp+1:igp+padx) = ff(igf+1:igf+padx)
2023 : : end do
2024 : : end do
2025 : : end do
2026 : :
2027 : 397880 : arr3(1) = nx
2028 : 397880 : arr3(2) = ny
2029 : 397880 : arr3(3) = nz
2030 : 397880 : status = DftiCreateDescriptor(Desc, DFTI_DOUBLE, DFTI_REAL, 3, arr3 )
2031 [ - + ]: 397880 : DFTI_CHECK(status)
2032 : :
2033 : 397880 : status = DftiSetValue(Desc, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX)
2034 : 397880 : status = DftiSetValue(Desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE )
2035 : 397880 : arr(1) = 0
2036 : 397880 : arr(2) = 1
2037 : 397880 : arr(3) = padx
2038 : 397880 : arr(4) = padx*ny
2039 : 397880 : status = DftiSetValue(Desc, DFTI_INPUT_STRIDES, arr)
2040 : 397880 : status = DftiSetValue(Desc, DFTI_INPUT_DISTANCE, nhp)
2041 : : arr(1) = 0
2042 : : arr(2) = 1
2043 : 397880 : arr(3) = ldx
2044 : 397880 : arr(4) = ldx*ldy
2045 : 397880 : status = DftiSetValue(Desc, DFTI_OUTPUT_STRIDES, arr)
2046 : 397880 : status = DftiSetValue(Desc, DFTI_OUTPUT_DISTANCE, ldx*ldy*ldz)
2047 : 397880 : status = DftiSetValue(Desc, DFTI_NUMBER_OF_TRANSFORMS, ndat)
2048 [ - + ]: 397880 : DFTI_CHECK(status)
2049 : :
2050 [ - + ]: 397880 : DFTI_CHECK( DftiCommitDescriptor(Desc) )
2051 : :
2052 [ - + ]: 397880 : DFTI_CHECK( DftiComputeBackward(Desc, ff_hp, gg) )
2053 : :
2054 [ - + ]: 397880 : DFTI_CHECK( DftiFreeDescriptor(Desc) )
2055 : :
2056 : 397880 : call dfti_free(cptr)
2057 : :
2058 : : #else
2059 : : ABI_ERROR("FFT_DFTI support not activated")
2060 : : ABI_UNUSED((/nx,ny,nz,ldx,ldy,ldz/))
2061 : : ABI_UNUSED(ff(1))
2062 : : ABI_UNUSED(gg(1))
2063 : : #endif
2064 : :
2065 : 397880 : end subroutine dfti_c2r_op_dpc
2066 : : !!***
2067 : :
2068 : : !----------------------------------------------------------------------
2069 : :
2070 : : !!****f* m_dfti/dfti_c2r_op_dp
2071 : : !! NAME
2072 : : !! dfti_c2r_op_dp
2073 : : !!
2074 : : !! FUNCTION
2075 : : !! Driver routine for out-of-place 3D complex-to-real FFT of lengths nx, ny, nz.
2076 : : !!
2077 : : !! INPUTS
2078 : : !! nx,ny,nz=Number of point along the three directions.
2079 : : !! ldx,ldy,ldz=Physical dimension of the f array (to avoid cache conflicts).
2080 : : !! ndat=Number of FFTs to be done.
2081 : : !! ff(2,ldx*ldy*ldz*ndat)=The complex array to be transformed.
2082 : : !!
2083 : : !! OUTPUT
2084 : : !! gg(ldx*ldy*ldz*ndat)=The backwards real FFT of ff.
2085 : : !!
2086 : : !! SOURCE
2087 : :
2088 : 397880 : subroutine dfti_c2r_op_dp(nx, ny, nz, ldx, ldy, ldz, ndat, ff, gg)
2089 : :
2090 : : !Arguments ------------------------------------
2091 : : !scalars
2092 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
2093 : : !arrays
2094 : : real(dp),target,intent(in) :: ff(2*ldx*ldy*ldz*ndat)
2095 : : real(dp),intent(inout) :: gg(ldx*ldy*ldz*ndat) !vz_i
2096 : :
2097 : : #ifdef HAVE_DFTI
2098 : : !Local variables-------------------------------
2099 : : !scalars
2100 : : type(C_ptr) :: ff_cptr
2101 : : !arrays
2102 : : complex(dp),contiguous, pointer :: ff_fptr(:)
2103 : : ! *************************************************************************
2104 : :
2105 : 397880 : ff_cptr = C_loc(ff)
2106 : 397880 : call C_F_pointer(ff_cptr,ff_fptr, shape=(/ldx*ldy*ldz*ndat/))
2107 : :
2108 : 397880 : call dfti_c2r_op_dpc(nx,ny,nz,ldx,ldy,ldz,ndat,ff_fptr,gg)
2109 : :
2110 : : #else
2111 : : ABI_ERROR("FFT_DFTI support not activated")
2112 : : ABI_UNUSED((/nx,ny,nz,ldx,ldy,ldz/))
2113 : : ABI_UNUSED((/ff(1),gg(1)/))
2114 : : #endif
2115 : :
2116 : 397880 : end subroutine dfti_c2r_op_dp
2117 : : !!***
2118 : :
2119 : : !----------------------------------------------------------------------
2120 : :
2121 : : !!****f* m_dfti/dfti_check_status
2122 : : !! NAME
2123 : : !! dfti_check_status
2124 : : !!
2125 : : !! FUNCTION
2126 : : !! Error handler for DFTI wrappers. Print error message and abort.
2127 : : !!
2128 : : !! INPUTS
2129 : : !! status = status error reported by dfti.
2130 : : !! [file] = file name
2131 : : !! [line] = line number
2132 : : !!
2133 : : !! SOURCE
2134 : :
2135 : : #ifdef HAVE_DFTI
2136 : :
2137 : 0 : subroutine dfti_check_status(status, file, line)
2138 : :
2139 : : !Arguments ------------------------------------
2140 : : !scalars
2141 : : integer,intent(in) :: status
2142 : : integer,optional,intent(in) :: line
2143 : : character(len=*),optional,intent(in) :: file
2144 : :
2145 : : !Local variables-------------------------------
2146 : : integer :: f90line
2147 : : character(len=10) :: lnum
2148 : : character(len=500) :: f90name
2149 : : character(len=500) :: my_msg
2150 : : character(len=DFTI_MAX_MESSAGE_LENGTH+500) :: err_msg
2151 : : ! *************************************************************************
2152 : :
2153 [ # # ]: 0 : if (PRESENT(line)) then
2154 : 0 : f90line=line
2155 : : else
2156 : 0 : f90line=0
2157 : : end if
2158 : 0 : call int2char10(f90line,lnum)
2159 : :
2160 [ # # ]: 0 : if (PRESENT(file)) then
2161 [ # # ]: 0 : f90name = basename(file)
2162 : : else
2163 : 0 : f90name='Subroutine Unknown'
2164 : : end if
2165 : :
2166 [ # # ]: 0 : my_msg=strcat(f90name,":",lnum,":")
2167 : :
2168 [ # # ]: 0 : if (status /= 0) then
2169 [ # # ]: 0 : if (.not. DftiErrorClass(status, DFTI_NO_ERROR)) then
2170 [ # # ]: 0 : err_msg = strcat(my_msg," Error: ",DftiErrorMessage(status))
2171 : 0 : ABI_ERROR(err_msg)
2172 : : end if
2173 : : end if
2174 : :
2175 : 0 : end subroutine dfti_check_status
2176 : : !!***
2177 : :
2178 : : #endif
2179 : :
2180 : : !----------------------------------------------------------------------
2181 : :
2182 : : !!****f* m_dfti/dfti_spawn_threads_here
2183 : : !! NAME
2184 : : !! dfti_spawn_threads_here
2185 : : !!
2186 : : !! FUNCTION
2187 : : !! Helper function that returns true if FFT calls should be OMP
2188 : : !! parallelized in the client code.
2189 : : !!
2190 : : !! INPUTS
2191 : : !! ndat=Number of FFT transforms to do
2192 : : !! nthreads = Number of threads available
2193 : : !!
2194 : : !! SOURCE
2195 : :
2196 : : function dfti_spawn_threads_here(ndat, nthreads) result(ans)
2197 : :
2198 : : !Arguments ------------------------------------
2199 : : !scalars
2200 : : integer,intent(in) :: ndat,nthreads
2201 : : logical :: ans
2202 : : ! *************************************************************************
2203 : :
2204 : 2287 : ans = .FALSE.
2205 : : #ifdef HAVE_OPENMP
2206 : : ans = (nthreads > 1 .and. MOD(ndat,nthreads) == 0 .and. .not. USE_LIB_THREADS)
2207 : : #else
2208 : : ABI_UNUSED((/ndat,nthreads/))
2209 : : #endif
2210 : :
2211 : : end function dfti_spawn_threads_here
2212 : : !!***
2213 : :
2214 : : !----------------------------------------------------------------------
2215 : :
2216 : : !!****f* m_dfti/dfti_use_lib_threads
2217 : : !! NAME
2218 : : !! dfti_use_lib_threads
2219 : : !!
2220 : : !! FUNCTION
2221 : : !!
2222 : : !! INPUTS
2223 : : !!
2224 : : !! SOURCE
2225 : :
2226 : 19 : subroutine dfti_use_lib_threads(logvar)
2227 : :
2228 : : !Arguments ------------------------------------
2229 : : logical,intent(in) :: logvar
2230 : : ! *************************************************************************
2231 : :
2232 : 19 : USE_LIB_THREADS = logvar
2233 : :
2234 : 19 : end subroutine dfti_use_lib_threads
2235 : : !!***
2236 : :
2237 : : !----------------------------------------------------------------------
2238 : :
2239 : : !!****f* m_dfti/dfti_alloc_real_dp
2240 : : !! NAME
2241 : : !! dfti_alloc_real_dp
2242 : : !!
2243 : : !! FUNCTION
2244 : : !!
2245 : : !! INPUTS
2246 : : !!
2247 : : !! OUTPUT
2248 : : !!
2249 : : !! SOURCE
2250 : :
2251 : : #ifdef HAVE_DFTI
2252 : :
2253 : : subroutine dfti_alloc_real_dp(size, cptr, fptr)
2254 : :
2255 : : !Arguments ------------------------------------
2256 : : !scalars
2257 : : integer,intent(in) :: size
2258 : : real(dp),contiguous, pointer :: fptr(:)
2259 : : type(C_PTR),intent(out) :: cptr
2260 : : ! *************************************************************************
2261 : :
2262 : : cptr = mkl_malloc( INT(size*C_DOUBLE, KIND=C_SIZE_T), DFTI_DEFAULT_ALIGNMENT_DP)
2263 : : if (.not. C_ASSOCIATED(cptr)) then
2264 : : ABI_ERROR("mkl_malloc returned NULL!")
2265 : : end if
2266 : :
2267 : : call c_f_pointer(cptr, fptr, [size])
2268 : :
2269 : : end subroutine dfti_alloc_real_dp
2270 : : !!***
2271 : :
2272 : : #endif
2273 : :
2274 : : !----------------------------------------------------------------------
2275 : :
2276 : : !!****f* m_dfti/dfti_alloc_complex_spc
2277 : : !! NAME
2278 : : !! dfti_alloc_complex_spc
2279 : : !!
2280 : : !! FUNCTION
2281 : : !!
2282 : : !! INPUTS
2283 : : !!
2284 : : !! OUTPUT
2285 : : !!
2286 : : !! SOURCE
2287 : :
2288 : : #ifdef HAVE_DFTI
2289 : :
2290 : 577567 : subroutine dfti_alloc_complex_spc(size, cptr, fptr)
2291 : :
2292 : : !Arguments ------------------------------------
2293 : : !scalars
2294 : : integer,intent(in) :: size
2295 : : complex(sp),contiguous, pointer :: fptr(:)
2296 : : type(C_PTR),intent(out) :: cptr
2297 : : ! *************************************************************************
2298 : :
2299 : 577567 : cptr = mkl_malloc( INT(2*size*C_FLOAT, KIND=C_SIZE_T), DFTI_DEFAULT_ALIGNMENT_SP)
2300 [ - + ]: 577567 : if (.not. C_ASSOCIATED(cptr)) then
2301 : 0 : ABI_ERROR("mkl_malloc returned NULL!")
2302 : : end if
2303 : :
2304 [ + + ]: 1155134 : call c_f_pointer(cptr, fptr, [size])
2305 : :
2306 : 577567 : end subroutine dfti_alloc_complex_spc
2307 : : !!***
2308 : :
2309 : : #endif
2310 : :
2311 : : !----------------------------------------------------------------------
2312 : :
2313 : : !!****f* m_dfti/dfti_alloc_complex_dpc
2314 : : !! NAME
2315 : : !! dfti_alloc_complex_dpc
2316 : : !!
2317 : : !! FUNCTION
2318 : : !!
2319 : : !! INPUTS
2320 : : !!
2321 : : !! OUTPUT
2322 : : !!
2323 : : !! SOURCE
2324 : :
2325 : : #ifdef HAVE_DFTI
2326 : :
2327 : 196445706 : subroutine dfti_alloc_complex_dpc(size, cptr, fptr)
2328 : :
2329 : : !Arguments ------------------------------------
2330 : : !scalars
2331 : : integer,intent(in) :: size
2332 : : complex(dp),contiguous, pointer :: fptr(:)
2333 : : type(C_PTR),intent(out) :: cptr
2334 : : ! *************************************************************************
2335 : :
2336 : 196445706 : cptr = mkl_malloc( INT(2*size*C_DOUBLE, KIND=C_SIZE_T), DFTI_DEFAULT_ALIGNMENT_DP)
2337 [ - + ]: 196445706 : if (.not. C_ASSOCIATED(cptr)) then
2338 : 0 : ABI_ERROR("mkl_malloc returned NULL!")
2339 : : end if
2340 : :
2341 [ + + ]: 392891412 : call c_f_pointer(cptr, fptr, [size])
2342 : :
2343 : 196445706 : end subroutine dfti_alloc_complex_dpc
2344 : : !!***
2345 : :
2346 : : #endif
2347 : :
2348 : : !----------------------------------------------------------------------
2349 : :
2350 : : END MODULE m_dfti
2351 : : !!***
|
