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 1248302 : 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 1248302 : 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 750271 : 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 1139992 : select case (isign)
288 : case (+1) ! G --> R
289 389721 : call dfti_c2r_op(nx,ny,nz,ldx,ldy,ldz,ndat,fofg,fofr)
290 : case (-1) ! R --> G
291 360550 : call dfti_r2c_op(nx,ny,nz,ldx,ldy,ldz,ndat,fofr,fofg)
292 : case default
293 750271 : ABI_BUG("Wrong isign")
294 : end select
295 :
296 : case default
297 1248302 : ABI_BUG("Wrong value for cplex")
298 : end select
299 :
300 1248302 : 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 59006717 : subroutine dfti_seqfourwf(cplex,denpot,fofgin,fofgout,fofr,gboundin,gboundout,istwf_k, &
369 59006717 : 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 160315703 : 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 59006717 : 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 59006717 : 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 59006717 : nx=ngfft(1); ny=ngfft(2); nz=ngfft(3)
404 59006717 : fftalg=ngfft(7); fftalga=fftalg/100; fftalgc=MOD(fftalg,10)
405 59006717 : fftcache=ngfft(8)
406 :
407 59006717 : use_fftrisc = (fftalgc==2)
408 59006717 : if (istwf_k==2.and.option==3) use_fftrisc = .FALSE.
409 162542388 : if (istwf_k>2.and.ANY(option==(/0,3/))) use_fftrisc = .FALSE.
410 :
411 : nthreads = xomp_get_num_threads(open_parallel=.TRUE.)
412 :
413 59006682 : if (use_fftrisc) then
414 : !call wrtout(std_out, " calls dfti_fftrisc")
415 59006666 : if (ndat == 1) then
416 58239625 : if (fftcore_mixprec == 0) then
417 : call dfti_fftrisc_dp(cplex,denpot,fofgin,fofgout,fofr,gboundin,gboundout,istwf_k,kg_kin,kg_kout,&
418 58047313 : 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 3342418 : do dat=1,ndat
463 2777902 : ptgin = 1 + (dat-1)*npwin
464 2777902 : ptgout = 1 + (dat-1)*npwout
465 3342418 : if (fftcore_mixprec == 0) then
466 : call dfti_fftrisc_dp(cplex,denpot,fofgin(1,ptgin),fofgout(1,ptgout),fofr,gboundin,gboundout,istwf_k,&
467 2777710 : 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 767041 : 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 59006717 : 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 62731627 : 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 214988 : subroutine dfti_fftug_dpc(fftalg, fftcache, npw_k, nx, ny, nz, ldx, ldy, ldz, ndat, &
1015 214988 : 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 214988 : real(dp),contiguous, pointer :: real_ug(:,:),real_ur(:,:)
1036 : ! *************************************************************************
1037 :
1038 214988 : iscale__ = 0; if (present(iscale)) iscale__ = iscale
1039 214988 : 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 214988 : 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 5735994 : 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 5735994 : 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 360550 : 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 360550 : integer,allocatable :: i1inver(:),i2inver(:),i3inver(:)
1814 : complex(dp),contiguous, pointer :: gg_hp(:)
1815 : ! *************************************************************************
1816 :
1817 360550 : padx = (nx/2+1)
1818 360550 : nhp = (nx/2+1)*ny*nz
1819 :
1820 360550 : call dfti_alloc_complex(nhp*ndat,cptr,gg_hp)
1821 :
1822 360550 : arr3(1) = nx
1823 360550 : arr3(2) = ny
1824 360550 : arr3(3) = nz
1825 360550 : status = DftiCreateDescriptor(Desc, DFTI_DOUBLE, DFTI_REAL, 3, arr3 )
1826 360550 : DFTI_CHECK(status)
1827 :
1828 360550 : status = DftiSetValue(Desc, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX)
1829 360550 : status = DftiSetValue(Desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE )
1830 360550 : arr(1) = 0
1831 360550 : arr(2) = 1
1832 360550 : arr(3) = ldx
1833 360550 : arr(4) = ldx*ldy
1834 360550 : status = DftiSetValue(Desc, DFTI_INPUT_STRIDES, arr)
1835 360550 : status = DftiSetValue(Desc, DFTI_INPUT_DISTANCE, ldx*ldy*ldz)
1836 : arr(1) = 0
1837 : arr(2) = 1
1838 360550 : arr(3) = padx
1839 360550 : arr(4) = padx*ny
1840 360550 : status = DftiSetValue(Desc, DFTI_OUTPUT_STRIDES, arr)
1841 360550 : status = DftiSetValue(Desc, DFTI_OUTPUT_DISTANCE, nhp)
1842 360550 : status = DftiSetValue(Desc, DFTI_NUMBER_OF_TRANSFORMS, ndat)
1843 360550 : status = DftiSetValue(Desc, DFTI_FORWARD_SCALE, one / DBLE(nx*ny*nz) )
1844 360550 : DFTI_CHECK(status)
1845 :
1846 360550 : DFTI_CHECK( DftiCommitDescriptor(Desc) )
1847 :
1848 360550 : DFTI_CHECK( DftiComputeForward(Desc, ff, gg_hp) )
1849 :
1850 360550 : DFTI_CHECK( DftiFreeDescriptor(Desc) )
1851 :
1852 : ! Reconstruct full FFT: Hermitian redundancy: out[i] is the conjugate of out[n-i]
1853 1081650 : ABI_MALLOC(i1inver,(padx))
1854 1081650 : ABI_MALLOC(i2inver,(ny))
1855 1081650 : ABI_MALLOC(i3inver,(nz))
1856 :
1857 360550 : i1inver(1)=1
1858 4407639 : do i1=2,padx
1859 4407639 : i1inver(i1)=nx+2-i1
1860 : end do
1861 :
1862 360550 : i2inver(1)=1
1863 7777186 : do i2=2,ny
1864 7777186 : i2inver(i2)=ny+2-i2
1865 : end do
1866 :
1867 360550 : i3inver(1)=1
1868 8925704 : do i3=2,nz
1869 8925704 : i3inver(i3)=nz+2-i3
1870 : end do
1871 :
1872 : igp=0
1873 721118 : do idat=1,ndat
1874 360568 : padatf=(idat-1)*ldx*ldy*ldz
1875 9647092 : do i3=1,nz
1876 8925974 : i3inv = i3inver(i3)
1877 233635581 : do i2=1,ny
1878 224349039 : i2inv = i2inver(i2)
1879 4486547570 : do i1=1,padx
1880 4253272557 : igp=igp+1
1881 4253272557 : igf = i1 + (i3-1)*ldx*ldy + (i2-1)*ldx + padatf
1882 4253272557 : gg(igf) = gg_hp(igp)
1883 4253272557 : i1inv = i1inver(i1)
1884 4477621596 : if (i1inv/=i1) then
1885 3828821252 : imgf = i1inv + (i3inv-1)*ldx*ldy + (i2inv-1)*ldx + padatf
1886 3828821252 : gg(imgf) = DCONJG(gg_hp(igp))
1887 : end if
1888 : end do
1889 : end do
1890 : end do
1891 : end do
1892 :
1893 360550 : ABI_FREE(i1inver)
1894 360550 : ABI_FREE(i2inver)
1895 360550 : ABI_FREE(i3inver)
1896 :
1897 360550 : 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 360550 : 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 360550 : 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 360550 : gg_cptr = C_loc(gg)
1947 360550 : call C_F_pointer(gg_cptr,gg_fptr, shape=(/ldx*ldy*ldz*ndat/))
1948 :
1949 360550 : 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 360550 : 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 389721 : 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 389721 : padx = (nx/2+1)
2010 389721 : nhp = padx*ny*nz
2011 :
2012 389721 : call dfti_alloc_complex(nhp*ndat,cptr,ff_hp)
2013 :
2014 : !!$OMP PARALLEL DO PRIVATE(padatf,padatp,igf,igp)
2015 779460 : do idat=1,ndat
2016 389739 : padatf=(idat-1)*ldx*ldy*ldz
2017 389739 : padatp=(idat-1)*padx*ny*nz
2018 10328073 : do i3=1,nz
2019 247273504 : do i2=1,ny
2020 237335152 : igf = (i3-1)*ldx*ldy + (i2-1)*ldx + padatf
2021 237335152 : igp = (i3-1)*padx*ny + (i2-1)*padx + padatp
2022 4628680725 : ff_hp(igp+1:igp+padx) = ff(igf+1:igf+padx)
2023 : end do
2024 : end do
2025 : end do
2026 :
2027 389721 : arr3(1) = nx
2028 389721 : arr3(2) = ny
2029 389721 : arr3(3) = nz
2030 389721 : status = DftiCreateDescriptor(Desc, DFTI_DOUBLE, DFTI_REAL, 3, arr3 )
2031 389721 : DFTI_CHECK(status)
2032 :
2033 389721 : status = DftiSetValue(Desc, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX)
2034 389721 : status = DftiSetValue(Desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE )
2035 389721 : arr(1) = 0
2036 389721 : arr(2) = 1
2037 389721 : arr(3) = padx
2038 389721 : arr(4) = padx*ny
2039 389721 : status = DftiSetValue(Desc, DFTI_INPUT_STRIDES, arr)
2040 389721 : status = DftiSetValue(Desc, DFTI_INPUT_DISTANCE, nhp)
2041 : arr(1) = 0
2042 : arr(2) = 1
2043 389721 : arr(3) = ldx
2044 389721 : arr(4) = ldx*ldy
2045 389721 : status = DftiSetValue(Desc, DFTI_OUTPUT_STRIDES, arr)
2046 389721 : status = DftiSetValue(Desc, DFTI_OUTPUT_DISTANCE, ldx*ldy*ldz)
2047 389721 : status = DftiSetValue(Desc, DFTI_NUMBER_OF_TRANSFORMS, ndat)
2048 389721 : DFTI_CHECK(status)
2049 :
2050 389721 : DFTI_CHECK( DftiCommitDescriptor(Desc) )
2051 :
2052 389721 : DFTI_CHECK( DftiComputeBackward(Desc, ff_hp, gg) )
2053 :
2054 389721 : DFTI_CHECK( DftiFreeDescriptor(Desc) )
2055 :
2056 389721 : 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 389721 : 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 389721 : 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 389721 : ff_cptr = C_loc(ff)
2106 389721 : call C_F_pointer(ff_cptr,ff_fptr, shape=(/ldx*ldy*ldz*ndat/))
2107 :
2108 389721 : 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 389721 : 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 2113 : 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 192320084 : 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 192320084 : cptr = mkl_malloc( INT(2*size*C_DOUBLE, KIND=C_SIZE_T), DFTI_DEFAULT_ALIGNMENT_DP)
2337 192320084 : if (.not. C_ASSOCIATED(cptr)) then
2338 0 : ABI_ERROR("mkl_malloc returned NULL!")
2339 : end if
2340 :
2341 384640168 : call c_f_pointer(cptr, fptr, [size])
2342 :
2343 192320084 : end subroutine dfti_alloc_complex_dpc
2344 : !!***
2345 :
2346 : #endif
2347 :
2348 : !----------------------------------------------------------------------
2349 :
2350 : END MODULE m_dfti
2351 : !!***
|
