(git:ccc2433)
ps_wavelet_kernel.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Creates the wavelet kernel for the wavelet based poisson solver.
10 !> \author Florian Schiffmann (09.2007,fschiff)
11 ! **************************************************************************************************
12 MODULE ps_wavelet_kernel
13 
14  USE kinds, ONLY: dp
15  USE message_passing, ONLY: mp_comm_type
16  USE ps_wavelet_base, ONLY: scramble_unpack
17  USE ps_wavelet_fft3d, ONLY: ctrig,&
18  ctrig_length,&
19  fftstp
20  USE ps_wavelet_scaling_function, ONLY: scaling_function,&
21  scf_recursion
22  USE ps_wavelet_util, ONLY: f_fft_dimensions,&
23  s_fft_dimensions
24 #include "../base/base_uses.f90"
25 
26  IMPLICIT NONE
27 
28  PRIVATE
29 
30  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ps_wavelet_kernel'
31 
32 ! *** Public data types ***
33 
34  PUBLIC :: createkernel
35 
36 CONTAINS
37 
38 ! **************************************************************************************************
39 !> \brief Allocate a pointer which corresponds to the zero-padded FFT slice needed for
40 !> calculating the convolution with the kernel expressed in the interpolating scaling
41 !> function basis. The kernel pointer is unallocated on input, allocated on output.
42 !> \param geocode Indicates the boundary conditions (BC) of the problem:
43 !> 'F' free BC, isolated systems.
44 !> The program calculates the solution as if the given density is
45 !> "alone" in R^3 space.
46 !> 'S' surface BC, isolated in y direction, periodic in xz plane
47 !> The given density is supposed to be periodic in the xz plane,
48 !> so the dimensions in these direction mus be compatible with the FFT
49 !> Beware of the fact that the isolated direction is y!
50 !> 'P' periodic BC.
51 !> The density is supposed to be periodic in all the three directions,
52 !> then all the dimensions must be compatible with the FFT.
53 !> No need for setting up the kernel.
54 !> \param n01 dimensions of the real space grid to be hit with the Poisson Solver
55 !> \param n02 dimensions of the real space grid to be hit with the Poisson Solver
56 !> \param n03 dimensions of the real space grid to be hit with the Poisson Solver
57 !> \param hx grid spacings. For the isolated BC case for the moment they are supposed to
58 !> be equal in the three directions
59 !> \param hy grid spacings. For the isolated BC case for the moment they are supposed to
60 !> be equal in the three directions
61 !> \param hz grid spacings. For the isolated BC case for the moment they are supposed to
62 !> be equal in the three directions
63 !> \param itype_scf order of the interpolating scaling functions used in the decomposition
64 !> \param iproc ,nproc: number of process, number of processes
65 !> \param nproc ...
66 !> \param kernel pointer for the kernel FFT. Unallocated on input, allocated on output.
67 !> Its dimensions are equivalent to the region of the FFT space for which the
68 !> kernel is injective. This will divide by two each direction,
69 !> since the kernel for the zero-padded convolution is real and symmetric.
70 !> \param mpi_group ...
71 !> \date February 2007
72 !> \author Luigi Genovese
73 !> \note Due to the fact that the kernel dimensions are unknown before the calling, the kernel
74 !> must be declared as pointer in input of this routine.
75 !> To avoid that, one can properly define the kernel dimensions by adding
76 !> the nd1,nd2,nd3 arguments to the PS_dim4allocation routine, then eliminating the pointer
77 !> declaration.
78 ! **************************************************************************************************
79  SUBROUTINE createkernel(geocode, n01, n02, n03, hx, hy, hz, itype_scf, iproc, nproc, kernel, mpi_group)
80 
81  CHARACTER(len=1), INTENT(in) :: geocode
82  INTEGER, INTENT(in) :: n01, n02, n03
83  REAL(kind=dp), INTENT(in) :: hx, hy, hz
84  INTEGER, INTENT(in) :: itype_scf, iproc, nproc
85  REAL(kind=dp), POINTER :: kernel(:)
86 
87  CLASS(mp_comm_type), INTENT(in) :: mpi_group
88 
89  INTEGER :: m1, m2, m3, md1, md2, md3, n1, n2, n3, &
90  nd1, nd2, nd3, nlimd, nlimk
91  REAL(kind=dp) :: hgrid
92 
93  hgrid = max(hx, hy, hz)
94 
95  IF (geocode == 'P') THEN
96 
97  CALL f_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, &
98  md1, md2, md3, nd1, nd2, nd3, nproc)
99 
100  ALLOCATE (kernel(1))
101  nlimd = n2
102  nlimk = 0
103 
104  ELSE IF (geocode == 'S') THEN
105 
106  CALL s_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, &
107  md1, md2, md3, nd1, nd2, nd3, nproc)
108 
109  ALLOCATE (kernel(nd1*nd2*nd3/nproc))
110 
111  !the kernel must be built and scattered to all the processes
112 
113  CALL surfaces_kernel(n1, n2, n3, m3, nd1, nd2, nd3, hx, hz, hy, &
114  itype_scf, kernel, iproc, nproc, mpi_group)
115  !last plane calculated for the density and the kernel
116 
117  nlimd = n2
118  nlimk = n3/2 + 1
119  ELSE IF (geocode == 'F') THEN
120 
121  !Build the Kernel
122 
123  CALL f_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, &
124  md1, md2, md3, nd1, nd2, nd3, nproc)
125  ALLOCATE (kernel(nd1*nd2*nd3/nproc))
126 
127  !the kernel must be built and scattered to all the processes
128  CALL free_kernel(n01, n02, n03, n1, n2, n3, nd1, nd2, nd3, hgrid, &
129  itype_scf, iproc, nproc, kernel, mpi_group)
130 
131  !last plane calculated for the density and the kernel
132  nlimd = n2/2
133  nlimk = n3/2 + 1
134 
135  ELSE
136 
137  cpabort("No wavelet based poisson solver for given geometry")
138 
139  END IF
140 !!! IF (iproc==0) THEN
141 !!! write(*,*)'done.'
142 !!! write(*,'(1x,a,i0)') 'Allocate words for kernel ',nd1*nd2*nd3/nproc
143 !!! !print the load balancing of the different dimensions on screen
144 !!! write(*,'(1x,a,3(i5))')'Grid Dimensions:',n01,n02,n03
145 !!! if (nproc > 1) then
146 !!! write(*,'(1x,a,3(i5),a,3(i5),a,3(i5))')&
147 !!! 'Memory occ. per proc. Density',md1,md3,md2/nproc,' Kernel',nd1,nd2,nd3/nproc
148 !!! write(*,'(1x,a)')'Load Balancing--------------------------------------------'
149 !!! jhd=1000
150 !!! jzd=1000
151 !!! npd=0
152 !!! load_balancing: do jproc=0,nproc-1
153 !!! !print *,'jproc,jfull=',jproc,jproc*md2/nproc,(jproc+1)*md2/nproc
154 !!! if ((jproc+1)*md2/nproc <= nlimd) then
155 !!! jfd=jproc
156 !!! else if (jproc*md2/nproc <= nlimd) then
157 !!! jhd=jproc
158 !!! npd=nint(real(nlimd-(jproc)*md2/nproc,KIND=dp)/real(md2/nproc,KIND=dp)*100._dp)
159 !!! else
160 !!! jzd=jproc
161 !!! exit load_balancing
162 !!! end if
163 !!! end do load_balancing
164 !!! write(*,'(1x,a,i3,a)')'LB_den : processors 0 -',jfd,' work at 100%'
165 !!! if (jfd < nproc-1) write(*,'(1x,a,i3,a,i3,1a)')' processor ',jhd,&
166 !!! ' work at ',npd,'%'
167 !!! if (jhd < nproc-1) write(*,'(1x,a,i3,1a,i3,a)')' processors ',&
168 !!! jzd,' -',nproc-1,' work at 0%'
169 !!! jhk=1000
170 !!! jzk=1000
171 !!! npk=0
172 !!! if (geocode /= 'P') then
173 !!! load_balancingk: do jproc=0,nproc-1
174 !!! !print *,'jproc,jfull=',jproc,jproc*nd3/nproc,(jproc+1)*nd3/nproc
175 !!! if ((jproc+1)*nd3/nproc <= nlimk) then
176 !!! jfk=jproc
177 !!! else if (jproc*nd3/nproc <= nlimk) then
178 !!! jhk=jproc
179 !!! npk=nint(real(nlimk-(jproc)*nd3/nproc,KIND=dp)/real(nd3/nproc,KIND=dp)*100._dp)
180 !!! else
181 !!! jzk=jproc
182 !!! exit load_balancingk
183 !!! end if
184 !!! end do load_balancingk
185 !!! write(*,'(1x,a,i3,a)')'LB_ker : processors 0 -',jfk,' work at 100%'
186 !!! if (jfk < nproc-1) write(*,'(1x,a,i3,a,i3,1a)')' processor ',jhk,&
187 !!! ' work at ',npk,'%'
188 !!! if (jhk < nproc-1) write(*,'(1x,a,i3,1a,i3,a)')' processors ',jzk,' -',nproc-1,&
189 !!! ' work at 0%'
190 !!! end if
191 !!! write(*,'(1x,a)')'The LB per processor is 1/3 LB_den + 2/3 LB_ker-----------'
192 !!! end if
193 !!!
194 !!! END IF
195  END SUBROUTINE createkernel
196 
197 ! **************************************************************************************************
198 !> \brief Build the kernel of the Poisson operator with
199 !> surfaces Boundary conditions
200 !> in an interpolating scaling functions basis.
201 !> Beware of the fact that the nonperiodic direction is y!
202 !> \param n1 Dimensions for the FFT
203 !> \param n2 Dimensions for the FFT
204 !> \param n3 Dimensions for the FFT
205 !> \param m3 Actual dimension in non-periodic direction
206 !> \param nker1 Dimensions of the kernel (nker3=n3/2+1) nker(1,2)=n(1,2)/2+1
207 !> \param nker2 Dimensions of the kernel (nker3=n3/2+1) nker(1,2)=n(1,2)/2+1
208 !> \param nker3 Dimensions of the kernel (nker3=n3/2+1) nker(1,2)=n(1,2)/2+1
209 !> \param h1 Mesh steps in the three dimensions
210 !> \param h2 Mesh steps in the three dimensions
211 !> \param h3 Mesh steps in the three dimensions
212 !> \param itype_scf Order of the scaling function
213 !> \param karray output array
214 !> \param iproc Number of process
215 !> \param nproc number of processes
216 !> \param mpi_group ...
217 !> \date October 2006
218 !> \author L. Genovese
219 ! **************************************************************************************************
220  SUBROUTINE surfaces_kernel(n1, n2, n3, m3, nker1, nker2, nker3, h1, h2, h3, &
221  itype_scf, karray, iproc, nproc, mpi_group)
222 
223  INTEGER, INTENT(in) :: n1, n2, n3, m3, nker1, nker2, nker3
224  REAL(kind=dp), INTENT(in) :: h1, h2, h3
225  INTEGER, INTENT(in) :: itype_scf, nproc, iproc
226  REAL(kind=dp), &
227  DIMENSION(nker1, nker2, nker3/nproc), &
228  INTENT(out) :: karray
229  TYPE(mp_comm_type), INTENT(in) :: mpi_group
230 
231  INTEGER, PARAMETER :: n_points = 2**6, ncache_optimal = 8*1024
232 
233  INTEGER :: i, i1, i2, i3, ic, iend, imu, ind1, ind2, inzee, ipolyord, ireim, istart, j2, &
234  j2nd, j2st, jnd1, jp2, jreim, n_cell, n_range, n_scf, nact2, ncache, nfft, num_of_mus, &
235  shift
236  INTEGER, ALLOCATABLE, DIMENSION(:) :: after, before, now
237  REAL(kind=dp) :: a, b, c, cp, d, diff, dx, fei, fer, foi, &
238  for, fr, mu1, pi, pion, ponx, pony, &
239  sp, value, x
240  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: kernel_scf, x_scf, y_scf
241  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: btrig, cossinarr
242  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: halfft_cache, kernel
243  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kernel_mpi
244  REAL(kind=dp), DIMENSION(9, 8) :: cpol
245 
246 !Better if higher (1024 points are enough 10^{-14}: 2*itype_scf*n_points)
247 ! include "perfdata.inc"
248 !FFT arrays
249 !coefficients for the polynomial interpolation
250 !assign the values of the coefficients
251 
252  karray = 0.0_dp
253  cpol(:, :) = 1._dp
254 
255  cpol(1, 2) = .25_dp
256 
257  cpol(1, 3) = 1._dp/3._dp
258 
259  cpol(1, 4) = 7._dp/12._dp
260  cpol(2, 4) = 8._dp/3._dp
261 
262  cpol(1, 5) = 19._dp/50._dp
263  cpol(2, 5) = 3._dp/2._dp
264 
265  cpol(1, 6) = 41._dp/272._dp
266  cpol(2, 6) = 27._dp/34._dp
267  cpol(3, 6) = 27._dp/272._dp
268 
269  cpol(1, 7) = 751._dp/2989._dp
270  cpol(2, 7) = 73._dp/61._dp
271  cpol(3, 7) = 27._dp/61._dp
272 
273  cpol(1, 8) = -989._dp/4540._dp
274  cpol(2, 8) = -1472._dp/1135._dp
275  cpol(3, 8) = 232._dp/1135._dp
276  cpol(4, 8) = -2624._dp/1135._dp
277 
278  !renormalize values
279  cpol(1, 1) = .5_dp*cpol(1, 1)
280  cpol(1:2, 2) = 2._dp/3._dp*cpol(1:2, 2)
281  cpol(1:2, 3) = 3._dp/8._dp*cpol(1:2, 3)
282  cpol(1:3, 4) = 2._dp/15._dp*cpol(1:3, 4)
283  cpol(1:3, 5) = 25._dp/144._dp*cpol(1:3, 5)
284  cpol(1:4, 6) = 34._dp/105._dp*cpol(1:4, 6)
285  cpol(1:4, 7) = 2989._dp/17280._dp*cpol(1:4, 7)
286  cpol(1:5, 8) = -454._dp/2835._dp*cpol(1:5, 8)
287 
288  !assign the complete values
289  cpol(2, 1) = cpol(1, 1)
290 
291  cpol(3, 2) = cpol(1, 2)
292 
293  cpol(3, 3) = cpol(2, 3)
294  cpol(4, 3) = cpol(1, 3)
295 
296  cpol(4, 4) = cpol(2, 4)
297  cpol(5, 4) = cpol(1, 4)
298 
299  cpol(4, 5) = cpol(3, 5)
300  cpol(5, 5) = cpol(2, 5)
301  cpol(6, 5) = cpol(1, 5)
302 
303  cpol(5, 6) = cpol(3, 6)
304  cpol(6, 6) = cpol(2, 6)
305  cpol(7, 6) = cpol(1, 6)
306 
307  cpol(5, 7) = cpol(4, 7)
308  cpol(6, 7) = cpol(3, 7)
309  cpol(7, 7) = cpol(2, 7)
310  cpol(8, 7) = cpol(1, 7)
311 
312  cpol(6, 8) = cpol(4, 8)
313  cpol(7, 8) = cpol(3, 8)
314  cpol(8, 8) = cpol(2, 8)
315  cpol(9, 8) = cpol(1, 8)
316 
317  !Number of integration points : 2*itype_scf*n_points
318  n_scf = 2*itype_scf*n_points
319  !Allocations
320  ALLOCATE (x_scf(0:n_scf))
321  ALLOCATE (y_scf(0:n_scf))
322 
323  !Build the scaling function
324  CALL scaling_function(itype_scf, n_scf, n_range, x_scf, y_scf)
325  !Step grid for the integration
326  dx = real(n_range, kind=dp)/real(n_scf, kind=dp)
327  !Extend the range (no more calculations because fill in by 0._dp)
328  n_cell = m3
329  n_range = max(n_cell, n_range)
330 
331  !Allocations
332  ncache = ncache_optimal
333  !the HalFFT must be performed only in the third dimension,
334  !and nker3=n3/2+1, hence
335  IF (ncache <= (nker3 - 1)*4) ncache = nker3 - 1*4
336 
337  !enlarge the second dimension of the kernel to be compatible with nproc
338  nact2 = nker2
339  enlarge_ydim: DO
340  IF (nproc*(nact2/nproc) /= nact2) THEN
341  nact2 = nact2 + 1
342  ELSE
343  EXIT enlarge_ydim
344  END IF
345  END DO enlarge_ydim
346 
347  !array for the MPI procedure
348  ALLOCATE (kernel(nker1, nact2/nproc, nker3))
349  ALLOCATE (kernel_mpi(nker1, nact2/nproc, nker3/nproc, nproc))
350  ALLOCATE (kernel_scf(n_range))
351  ALLOCATE (halfft_cache(2, ncache/4, 2))
352  ALLOCATE (cossinarr(2, n3/2 - 1))
353  ALLOCATE (btrig(2, ctrig_length))
354  ALLOCATE (after(7))
355  ALLOCATE (now(7))
356  ALLOCATE (before(7))
357 
358  !constants
359  pi = 4._dp*atan(1._dp)
360 
361  !arrays for the halFFT
362  CALL ctrig(n3/2, btrig, after, before, now, 1, ic)
363 
364  !build the phases for the HalFFT reconstruction
365  pion = 2._dp*pi/real(n3, kind=dp)
366  DO i3 = 2, n3/2
367  x = real(i3 - 1, kind=dp)*pion
368  cossinarr(1, i3 - 1) = cos(x)
369  cossinarr(2, i3 - 1) = -sin(x)
370  END DO
371 
372  ! satisfy valgrind, init arrays to large value, even if the offending bit is (likely?) padding
373  kernel = huge(0._dp)
374  kernel_mpi = huge(0._dp)
375 
376  !calculate the limits of the FFT calculations
377  !that can be performed in a row remaining inside the cache
378  num_of_mus = ncache/(2*n3)
379 
380  diff = 0._dp
381  !order of the polynomial to be used for integration (must be a power of two)
382  ipolyord = 8 !this part should be incorporated inside the numerical integration
383  !here we have to choice the piece of the x-y grid to cover
384 
385  !let us now calculate the fraction of mu that will be considered
386  j2st = iproc*(nact2/nproc)
387  j2nd = min((iproc + 1)*(nact2/nproc), n2/2 + 1)
388 
389  DO ind2 = (n1/2 + 1)*j2st + 1, (n1/2 + 1)*j2nd, num_of_mus
390  istart = ind2
391  iend = min(ind2 + (num_of_mus - 1), (n1/2 + 1)*j2nd)
392  nfft = iend - istart + 1
393  shift = 0
394 
395  !initialization of the interesting part of the cache array
396  halfft_cache(:, :, :) = 0._dp
397 
398  IF (istart == 1) THEN
399  !i2=1
400  shift = 1
401 
402  CALL calculates_green_opt_muzero(n_range, n_scf, ipolyord, x_scf, y_scf, &
403  cpol(1, ipolyord), dx, kernel_scf)
404 
405  !copy of the first zero value
406  halfft_cache(1, 1, 1) = 0._dp
407 
408  DO i3 = 1, m3
409 
410  value = 0.5_dp*h3*kernel_scf(i3)
411  !index in where to copy the value of the kernel
412  CALL indices(ireim, num_of_mus, n3/2 + i3, 1, ind1)
413  !index in where to copy the symmetric value
414  CALL indices(jreim, num_of_mus, n3/2 + 2 - i3, 1, jnd1)
415  halfft_cache(ireim, ind1, 1) = value
416  halfft_cache(jreim, jnd1, 1) = value
417 
418  END DO
419 
420  END IF
421 
422  loopimpulses: DO imu = istart + shift, iend
423 
424  !here there is the value of mu associated to hgrid
425  !note that we have multiplicated mu for hgrid to be comparable
426  !with mu0ref
427 
428  !calculate the proper value of mu taking into account the periodic dimensions
429  !corresponding value of i1 and i2
430  i1 = mod(imu, n1/2 + 1)
431  IF (i1 == 0) i1 = n1/2 + 1
432  i2 = (imu - i1)/(n1/2 + 1) + 1
433  ponx = real(i1 - 1, kind=dp)/real(n1, kind=dp)
434  pony = real(i2 - 1, kind=dp)/real(n2, kind=dp)
435 
436  mu1 = 2._dp*pi*sqrt((ponx/h1)**2 + (pony/h2)**2)*h3
437 
438  CALL calculates_green_opt(n_range, n_scf, itype_scf, ipolyord, x_scf, y_scf, &
439  cpol(1, ipolyord), mu1, dx, kernel_scf)
440 
441  !readjust the coefficient and define the final kernel
442 
443  !copy of the first zero value
444  halfft_cache(1, imu - istart + 1, 1) = 0._dp
445  DO i3 = 1, m3
446  value = -0.5_dp*h3/mu1*kernel_scf(i3)
447  !write(80,*)mu1,i3,kernel_scf(i03)
448  !index in where to copy the value of the kernel
449  CALL indices(ireim, num_of_mus, n3/2 + i3, imu - istart + 1, ind1)
450  !index in where to copy the symmetric value
451  CALL indices(jreim, num_of_mus, n3/2 + 2 - i3, imu - istart + 1, jnd1)
452  halfft_cache(ireim, ind1, 1) = value
453  halfft_cache(jreim, jnd1, 1) = value
454  END DO
455 
456  END DO loopimpulses
457 
458  !now perform the FFT of the array in cache
459  inzee = 1
460  DO i = 1, ic
461  CALL fftstp(num_of_mus, nfft, n3/2, num_of_mus, n3/2, &
462  halfft_cache(1, 1, inzee), halfft_cache(1, 1, 3 - inzee), &
463  btrig, after(i), now(i), before(i), 1)
464  inzee = 3 - inzee
465  END DO
466  !assign the values of the FFT array
467  !and compare with the good results
468  DO imu = istart, iend
469 
470  !corresponding value of i1 and i2
471  i1 = mod(imu, n1/2 + 1)
472  IF (i1 == 0) i1 = n1/2 + 1
473  i2 = (imu - i1)/(n1/2 + 1) + 1
474 
475  j2 = i2 - j2st
476 
477  a = halfft_cache(1, imu - istart + 1, inzee)
478  b = halfft_cache(2, imu - istart + 1, inzee)
479  kernel(i1, j2, 1) = a + b
480  kernel(i1, j2, n3/2 + 1) = a - b
481 
482  DO i3 = 2, n3/2
483  ind1 = imu - istart + 1 + num_of_mus*(i3 - 1)
484  jnd1 = imu - istart + 1 + num_of_mus*(n3/2 + 2 - i3 - 1)
485  cp = cossinarr(1, i3 - 1)
486  sp = cossinarr(2, i3 - 1)
487  a = halfft_cache(1, ind1, inzee)
488  b = halfft_cache(2, ind1, inzee)
489  c = halfft_cache(1, jnd1, inzee)
490  d = halfft_cache(2, jnd1, inzee)
491  fer = .5_dp*(a + c)
492  fei = .5_dp*(b - d)
493  for = .5_dp*(a - c)
494  foi = .5_dp*(b + d)
495  fr = fer + cp*foi - sp*for
496  kernel(i1, j2, i3) = fr
497  END DO
498  END DO
499 
500  END DO
501 
502  !give to each processor a slice of the third dimension
503  IF (nproc > 1) THEN
504  CALL mpi_group%alltoall(kernel, &!nker1*(nact2/nproc)*(nker3/nproc), &
505  kernel_mpi, nker1*(nact2/nproc)*(nker3/nproc))
506 
507  DO jp2 = 1, nproc
508  DO i3 = 1, nker3/nproc
509  DO i2 = 1, nact2/nproc
510  j2 = i2 + (jp2 - 1)*(nact2/nproc)
511  IF (j2 <= nker2) THEN
512  DO i1 = 1, nker1
513  karray(i1, j2, i3) = &
514  kernel_mpi(i1, i2, i3, jp2)
515  END DO
516  END IF
517  END DO
518  END DO
519  END DO
520 
521  ELSE
522  karray(1:nker1, 1:nker2, 1:nker3) = kernel(1:nker1, 1:nker2, 1:nker3)
523  END IF
524 
525  !De-allocations
526  DEALLOCATE (kernel)
527  DEALLOCATE (kernel_mpi)
528  DEALLOCATE (btrig)
529  DEALLOCATE (after)
530  DEALLOCATE (now)
531  DEALLOCATE (before)
532  DEALLOCATE (halfft_cache)
533  DEALLOCATE (kernel_scf)
534  DEALLOCATE (x_scf)
535  DEALLOCATE (y_scf)
536 
537  END SUBROUTINE surfaces_kernel
538 
539 ! **************************************************************************************************
540 !> \brief ...
541 !> \param n ...
542 !> \param n_scf ...
543 !> \param itype_scf ...
544 !> \param intorder ...
545 !> \param xval ...
546 !> \param yval ...
547 !> \param c ...
548 !> \param mu ...
549 !> \param hres ...
550 !> \param g_mu ...
551 ! **************************************************************************************************
552  SUBROUTINE calculates_green_opt(n, n_scf, itype_scf, intorder, xval, yval, c, mu, hres, g_mu)
553  INTEGER, INTENT(in) :: n, n_scf, itype_scf, intorder
554  REAL(kind=dp), DIMENSION(0:n_scf), INTENT(in) :: xval, yval
555  REAL(kind=dp), DIMENSION(intorder+1), INTENT(in) :: c
556  REAL(kind=dp), INTENT(in) :: mu, hres
557  REAL(kind=dp), DIMENSION(n), INTENT(out) :: g_mu
558 
559  REAL(kind=dp), PARAMETER :: mu_max = 0.2_dp
560 
561  INTEGER :: i, iend, ikern, ivalue, izero, n_iter, &
562  nrec
563  REAL(kind=dp) :: f, filter, fl, fr, gleft, gltmp, gright, &
564  grtmp, mu0, ratio, x, x0, x1
565  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: green, green1
566 
567  g_mu = 0.0_dp
568  !We calculate the number of iterations to go from mu0 to mu0_ref
569  IF (mu <= mu_max) THEN
570  n_iter = 0
571  mu0 = mu
572  ELSE
573  n_iter = 1
574  loop_iter: DO
575  ratio = real(2**n_iter, kind=dp)
576  mu0 = mu/ratio
577  IF (mu0 <= mu_max) THEN
578  EXIT loop_iter
579  END IF
580  n_iter = n_iter + 1
581  END DO loop_iter
582  END IF
583 
584  !dimension needed for the correct calculation of the recursion
585  nrec = 2**n_iter*n
586 
587  ALLOCATE (green(-nrec:nrec))
588 
589  !initialization of the branching value
590  ikern = 0
591  izero = 0
592  initialization: DO
593  IF (xval(izero) >= real(ikern, kind=dp) .OR. izero == n_scf) EXIT initialization
594  izero = izero + 1
595  END DO initialization
596  green = 0._dp
597  !now perform the interpolation in right direction
598  ivalue = izero
599  gright = 0._dp
600  loop_right: DO
601  IF (ivalue >= n_scf - intorder - 1) EXIT loop_right
602  DO i = 1, intorder + 1
603  x = xval(ivalue) - real(ikern, kind=dp)
604  f = yval(ivalue)*exp(-mu0*x)
605  filter = intorder*c(i)
606  gright = gright + filter*f
607  ivalue = ivalue + 1
608  END DO
609  ivalue = ivalue - 1
610  END DO loop_right
611  iend = n_scf - ivalue
612  DO i = 1, iend
613  x = xval(ivalue) - real(ikern, kind=dp)
614  f = yval(ivalue)*exp(-mu0*x)
615  filter = intorder*c(i)
616  gright = gright + filter*f
617  ivalue = ivalue + 1
618  END DO
619  gright = hres*gright
620 
621  !the scaling function is symmetric, so the same for the other direction
622  gleft = gright
623 
624  green(ikern) = gleft + gright
625 
626  !now the loop until the last value
627  DO ikern = 1, nrec
628  gltmp = 0._dp
629  grtmp = 0._dp
630  ivalue = izero
631  x0 = xval(izero)
632  loop_integration: DO
633  IF (izero == n_scf) EXIT loop_integration
634  DO i = 1, intorder + 1
635  x = xval(ivalue)
636  fl = yval(ivalue)*exp(mu0*x)
637  fr = yval(ivalue)*exp(-mu0*x)
638  filter = intorder*c(i)
639  gltmp = gltmp + filter*fl
640  grtmp = grtmp + filter*fr
641  ivalue = ivalue + 1
642  IF (xval(izero) >= real(ikern, kind=dp) .OR. izero == n_scf) THEN
643  x1 = xval(izero)
644  EXIT loop_integration
645  END IF
646  izero = izero + 1
647  END DO
648  ivalue = ivalue - 1
649  izero = izero - 1
650  END DO loop_integration
651  gleft = exp(-mu0)*(gleft + hres*exp(-mu0*real(ikern - 1, kind=dp))*gltmp)
652  IF (izero == n_scf) THEN
653  gright = 0._dp
654  ELSE
655  gright = exp(mu0)*(gright - hres*exp(mu0*real(ikern - 1, kind=dp))*grtmp)
656  END IF
657  green(ikern) = gleft + gright
658  green(-ikern) = gleft + gright
659  IF (abs(green(ikern)) <= 1.e-20_dp) THEN
660  nrec = ikern
661  EXIT
662  END IF
663  !print *,ikern,izero,n_scf,gltmp,grtmp,gleft,gright,x0,x1,green(ikern)
664  END DO
665  !now we must calculate the recursion
666  ALLOCATE (green1(-nrec:nrec))
667 
668  !Start the iteration to go from mu0 to mu
669  CALL scf_recursion(itype_scf, n_iter, nrec, green(-nrec), green1(-nrec))
670 
671  DO i = 1, min(n, nrec)
672  g_mu(i) = green(i - 1)
673  END DO
674  DO i = min(n, nrec) + 1, n
675  g_mu(i) = 0._dp
676  END DO
677 
678  DEALLOCATE (green, green1)
679 
680  END SUBROUTINE calculates_green_opt
681 
682 ! **************************************************************************************************
683 !> \brief ...
684 !> \param n ...
685 !> \param n_scf ...
686 !> \param intorder ...
687 !> \param xval ...
688 !> \param yval ...
689 !> \param c ...
690 !> \param hres ...
691 !> \param green ...
692 ! **************************************************************************************************
693  SUBROUTINE calculates_green_opt_muzero(n, n_scf, intorder, xval, yval, c, hres, green)
694  INTEGER, INTENT(in) :: n, n_scf, intorder
695  REAL(kind=dp), DIMENSION(0:n_scf), INTENT(in) :: xval, yval
696  REAL(kind=dp), DIMENSION(intorder+1), INTENT(in) :: c
697  REAL(kind=dp), INTENT(in) :: hres
698  REAL(kind=dp), DIMENSION(n), INTENT(out) :: green
699 
700  INTEGER :: i, iend, ikern, ivalue, izero
701  REAL(kind=dp) :: c0, c1, filter, gl0, gl1, gr0, gr1, x, y
702 
703 !initialization of the branching value
704 
705  ikern = 0
706  izero = 0
707  initialization: DO
708  IF (xval(izero) >= real(ikern, kind=dp) .OR. izero == n_scf) EXIT initialization
709  izero = izero + 1
710  END DO initialization
711  green = 0._dp
712  !first case, ikern=0
713  !now perform the interpolation in right direction
714  ivalue = izero
715  gr1 = 0._dp
716  loop_right: DO
717  IF (ivalue >= n_scf - intorder - 1) EXIT loop_right
718  DO i = 1, intorder + 1
719  x = xval(ivalue)
720  y = yval(ivalue)
721  filter = intorder*c(i)
722  gr1 = gr1 + filter*x*y
723  ivalue = ivalue + 1
724  END DO
725  ivalue = ivalue - 1
726  END DO loop_right
727  iend = n_scf - ivalue
728  DO i = 1, iend
729  x = xval(ivalue)
730  y = yval(ivalue)
731  filter = intorder*c(i)
732  gr1 = gr1 + filter*x*y
733  ivalue = ivalue + 1
734  END DO
735  gr1 = hres*gr1
736  !the scaling function is symmetric
737  gl1 = -gr1
738  gl0 = 0.5_dp
739  gr0 = 0.5_dp
740 
741  green(1) = 2._dp*gr1
742 
743  !now the loop until the last value
744  DO ikern = 1, n - 1
745  c0 = 0._dp
746  c1 = 0._dp
747  ivalue = izero
748  loop_integration: DO
749  IF (izero == n_scf) EXIT loop_integration
750  DO i = 1, intorder + 1
751  x = xval(ivalue)
752  y = yval(ivalue)
753  filter = intorder*c(i)
754  c0 = c0 + filter*y
755  c1 = c1 + filter*y*x
756  ivalue = ivalue + 1
757  IF (xval(izero) >= real(ikern, kind=dp) .OR. izero == n_scf) THEN
758  EXIT loop_integration
759  END IF
760  izero = izero + 1
761  END DO
762  ivalue = ivalue - 1
763  izero = izero - 1
764  END DO loop_integration
765  c0 = hres*c0
766  c1 = hres*c1
767 
768  gl0 = gl0 + c0
769  gl1 = gl1 + c1
770  gr0 = gr0 - c0
771  gr1 = gr1 - c1
772  !general case
773  green(ikern + 1) = real(ikern, kind=dp)*(gl0 - gr0) + gr1 - gl1
774  !print *,ikern,izero,n_scf,gltmp,grtmp,gleft,gright,x0,x1,green(ikern)
775  END DO
776 
777  END SUBROUTINE calculates_green_opt_muzero
778 
779 ! **************************************************************************************************
780 !> \brief ...
781 !> \param var_realimag ...
782 !> \param nelem ...
783 !> \param intrn ...
784 !> \param extrn ...
785 !> \param index ...
786 ! **************************************************************************************************
787  SUBROUTINE indices(var_realimag, nelem, intrn, extrn, index)
788 
789  INTEGER, INTENT(out) :: var_realimag
790  INTEGER, INTENT(in) :: nelem, intrn, extrn
791  INTEGER, INTENT(out) :: index
792 
793  INTEGER :: i
794 
795 !real or imaginary part
796 
797  var_realimag = 2 - mod(intrn, 2)
798 !actual index over half the length
799 
800  i = (intrn + 1)/2
801  !check
802  IF (2*(i - 1) + var_realimag /= intrn) THEN
803  print *, 'error, index=', intrn, 'var_realimag=', var_realimag, 'i=', i
804  END IF
805  !complete index to be assigned
806  index = extrn + nelem*(i - 1)
807 
808  END SUBROUTINE indices
809 
810 ! **************************************************************************************************
811 !> \brief Build the kernel of a gaussian function
812 !> for interpolating scaling functions.
813 !> Do the parallel HalFFT of the symmetrized function and stores into
814 !> memory only 1/8 of the grid divided by the number of processes nproc
815 !>
816 !> Build the kernel (karray) of a gaussian function
817 !> for interpolating scaling functions
818 !> $$ K(j) = \sum_k \omega_k \int \int \phi(x) g_k(x'-x) \delta(x'- j) dx dx' $$
819 !> \param n01 Mesh dimensions of the density
820 !> \param n02 Mesh dimensions of the density
821 !> \param n03 Mesh dimensions of the density
822 !> \param nfft1 Dimensions of the FFT grid (HalFFT in the third direction)
823 !> \param nfft2 Dimensions of the FFT grid (HalFFT in the third direction)
824 !> \param nfft3 Dimensions of the FFT grid (HalFFT in the third direction)
825 !> \param n1k Dimensions of the kernel FFT
826 !> \param n2k Dimensions of the kernel FFT
827 !> \param n3k Dimensions of the kernel FFT
828 !> \param hgrid Mesh step
829 !> \param itype_scf Order of the scaling function (8,14,16)
830 !> \param iproc ...
831 !> \param nproc ...
832 !> \param karray ...
833 !> \param mpi_group ...
834 !> \date February 2006
835 !> \author T. Deutsch, L. Genovese
836 ! **************************************************************************************************
837  SUBROUTINE free_kernel(n01, n02, n03, nfft1, nfft2, nfft3, n1k, n2k, n3k, &
838  hgrid, itype_scf, iproc, nproc, karray, mpi_group)
839 
840  INTEGER, INTENT(in) :: n01, n02, n03, nfft1, nfft2, nfft3, n1k, &
841  n2k, n3k
842  REAL(kind=dp), INTENT(in) :: hgrid
843  INTEGER, INTENT(in) :: itype_scf, iproc, nproc
844  REAL(kind=dp), DIMENSION(n1k, n2k, n3k/nproc), &
845  INTENT(out) :: karray
846  TYPE(mp_comm_type), INTENT(in) :: mpi_group
847 
848  INTEGER, PARAMETER :: n_gauss = 89, n_points = 2**6
849  REAL(kind=dp), PARAMETER :: p0_ref = 1._dp
850 
851  INTEGER :: i, i01, i02, i03, i1, i2, i3, i_gauss, &
852  i_kern, iend, istart, istart1, n1h, &
853  n2h, n3h, n_cell, n_iter, n_range, &
854  n_scf, nker1, nker2, nker3
855  REAL(kind=dp) :: a1, a2, a3, a_range, absci, acc_gauss, &
856  dr_gauss, dx, factor, factor2, kern, &
857  p0_cell, p0gauss, pgauss, ur_gauss
858  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: kern_1_scf, kernel_scf, x_scf, y_scf
859  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: kp
860  REAL(kind=dp), DIMENSION(n_gauss) :: p_gauss, w_gauss
861 
862 !Do not touch !!!!
863 !Better if higher (1024 points are enough 10^{-14}: 2*itype_scf*n_points)
864 !Better p_gauss for calculation
865 !(the support of the exponential should be inside [-n_range/2,n_range/2])
866 !Number of integration points : 2*itype_scf*n_points
867 
868  n_scf = 2*itype_scf*n_points
869  !Set karray
870  karray = 0.0_dp
871  !here we must set the dimensions for the fft part, starting from the nfft
872  !remember that actually nfft2 is associated to n03 and viceversa
873 
874  !dimensions that define the center of symmetry
875  n1h = nfft1/2
876  n2h = nfft2/2
877  n3h = nfft3/2
878 
879  !Auxiliary dimensions only for building the FFT part
880  nker1 = nfft1
881  nker2 = nfft2
882  nker3 = nfft3/2 + 1
883 
884  !adjusting the last two dimensions to be multiples of nproc
885  DO
886  IF (modulo(nker2, nproc) == 0) EXIT
887  nker2 = nker2 + 1
888  END DO
889  DO
890  IF (modulo(nker3, nproc) == 0) EXIT
891  nker3 = nker3 + 1
892  END DO
893 
894  !this will be the array of the kernel in the real space
895  ALLOCATE (kp(n1h + 1, n3h + 1, nker2/nproc))
896 
897  !defining proper extremes for the calculation of the
898  !local part of the kernel
899 
900  istart = iproc*nker2/nproc + 1
901  iend = min((iproc + 1)*nker2/nproc, n2h + n03)
902 
903  istart1 = istart
904  IF (iproc .EQ. 0) istart1 = n2h - n03 + 2
905 
906  !Allocations
907  ALLOCATE (x_scf(0:n_scf))
908  ALLOCATE (y_scf(0:n_scf))
909 
910  !Build the scaling function
911  CALL scaling_function(itype_scf, n_scf, n_range, x_scf, y_scf)
912  !Step grid for the integration
913  dx = real(n_range, kind=dp)/real(n_scf, kind=dp)
914  !Extend the range (no more calculations because fill in by 0._dp)
915  n_cell = max(n01, n02, n03)
916  n_range = max(n_cell, n_range)
917 
918  !Allocations
919  ALLOCATE (kernel_scf(-n_range:n_range))
920  ALLOCATE (kern_1_scf(-n_range:n_range))
921 
922  !Lengthes of the box (use FFT dimension)
923  a1 = hgrid*real(n01, kind=dp)
924  a2 = hgrid*real(n02, kind=dp)
925  a3 = hgrid*real(n03, kind=dp)
926 
927  x_scf(:) = hgrid*x_scf(:)
928  y_scf(:) = 1._dp/hgrid*y_scf(:)
929  dx = hgrid*dx
930  !To have a correct integration
931  p0_cell = p0_ref/(hgrid*hgrid)
932 
933  !Initialization of the gaussian (Beylkin)
934  CALL gequad(p_gauss, w_gauss, ur_gauss, dr_gauss, acc_gauss)
935  !In order to have a range from a_range=sqrt(a1*a1+a2*a2+a3*a3)
936  !(biggest length in the cube)
937  !We divide the p_gauss by a_range**2 and a_gauss by a_range
938  a_range = sqrt(a1*a1 + a2*a2 + a3*a3)
939  factor = 1._dp/a_range
940  !factor2 = factor*factor
941  factor2 = 1._dp/(a1*a1 + a2*a2 + a3*a3)
942  DO i_gauss = 1, n_gauss
943  p_gauss(i_gauss) = factor2*p_gauss(i_gauss)
944  END DO
945  DO i_gauss = 1, n_gauss
946  w_gauss(i_gauss) = factor*w_gauss(i_gauss)
947  END DO
948 
949  kp(:, :, :) = 0._dp
950  !Use in this order (better for accuracy).
951  loop_gauss: DO i_gauss = n_gauss, 1, -1
952  !Gaussian
953  pgauss = p_gauss(i_gauss)
954 
955  !We calculate the number of iterations to go from pgauss to p0_ref
956  n_iter = nint((log(pgauss) - log(p0_cell))/log(4._dp))
957  IF (n_iter <= 0) THEN
958  n_iter = 0
959  p0gauss = pgauss
960  ELSE
961  p0gauss = pgauss/4._dp**n_iter
962  END IF
963 
964  !Stupid integration
965  !Do the integration with the exponential centered in i_kern
966  kernel_scf(:) = 0._dp
967  DO i_kern = 0, n_range
968  kern = 0._dp
969  DO i = 0, n_scf
970  absci = x_scf(i) - real(i_kern, kind=dp)*hgrid
971  absci = absci*absci
972  kern = kern + y_scf(i)*exp(-p0gauss*absci)*dx
973  END DO
974  kernel_scf(i_kern) = kern
975  kernel_scf(-i_kern) = kern
976  IF (abs(kern) < 1.e-18_dp) THEN
977  !Too small not useful to calculate
978  EXIT
979  END IF
980  END DO
981 
982  !Start the iteration to go from p0gauss to pgauss
983  CALL scf_recursion(itype_scf, n_iter, n_range, kernel_scf, kern_1_scf)
984 
985  !Add to the kernel (only the local part)
986 
987  DO i3 = istart1, iend
988  i03 = i3 - n2h - 1
989  ! Crash if index out of range
990  ! Without compiler bounds checking, the calculation might run successfully but
991  ! it is also possible that the Hartree energy will blow up
992  ! This seems to happen with large MPI processor counts if the size of the
993  ! RS grid is not directly compatible with the allowed FFT dimensions in
994  ! subroutine fourier_dim (ps_wavelet_fft3d.F)
995  IF (i03 .LT. -n_range .OR. i03 .GT. n_range) THEN
996  CALL cp_abort(__location__, "Index out of range in wavelet solver. "// &
997  "Try decreasing the number of MPI processors, or adjust the PW_CUTOFF or cell size "// &
998  "so that 2*(number of RS grid points) matches the allowed FFT dimensions "// &
999  "(see ps_wavelet_fft3d.F) exactly.")
1000  END IF
1001  DO i2 = 1, n02
1002  i02 = i2 - 1
1003  DO i1 = 1, n01
1004  i01 = i1 - 1
1005  kp(i1, i2, i3 - istart + 1) = kp(i1, i2, i3 - istart + 1) + w_gauss(i_gauss)* &
1006  kernel_scf(i01)*kernel_scf(i02)*kernel_scf(i03)
1007  END DO
1008  END DO
1009  END DO
1010 
1011  END DO loop_gauss
1012 
1013  !De-allocations
1014  DEALLOCATE (kernel_scf)
1015  DEALLOCATE (kern_1_scf)
1016  DEALLOCATE (x_scf)
1017  DEALLOCATE (y_scf)
1018 
1019 !!!!END KERNEL CONSTRUCTION
1020 
1021 !!$ if(iproc .eq. 0) print *,"Do a 3D PHalFFT for the kernel"
1022 
1023  CALL kernelfft(nfft1, nfft2, nfft3, nker1, nker2, nker3, n1k, n2k, n3k, nproc, iproc, &
1024  kp, karray, mpi_group)
1025 
1026  !De-allocations
1027  DEALLOCATE (kp)
1028 
1029  END SUBROUTINE free_kernel
1030 
1031 ! **************************************************************************************************
1032 !> \brief ...
1033 !> \param n1 ...
1034 !> \param n3 ...
1035 !> \param lot ...
1036 !> \param nfft ...
1037 !> \param i1 ...
1038 !> \param zf ...
1039 !> \param zw ...
1040 ! **************************************************************************************************
1041  SUBROUTINE inserthalf(n1, n3, lot, nfft, i1, zf, zw)
1042  INTEGER, INTENT(in) :: n1, n3, lot, nfft, i1
1043  REAL(kind=dp), DIMENSION(n1/2+1, n3/2+1), &
1044  INTENT(in) :: zf
1045  REAL(kind=dp), DIMENSION(2, lot, n3/2), &
1046  INTENT(out) :: zw
1047 
1048  INTEGER :: i01, i03i, i03r, i3, l1, l3
1049 
1050  zw = 0.0_dp
1051  i3 = 0
1052  DO l3 = 1, n3, 2
1053  i3 = i3 + 1
1054  i03r = abs(l3 - n3/2 - 1) + 1
1055  i03i = abs(l3 - n3/2) + 1
1056  DO l1 = 1, nfft
1057  i01 = abs(l1 - 1 + i1 - n1/2 - 1) + 1
1058  zw(1, l1, i3) = zf(i01, i03r)
1059  zw(2, l1, i3) = zf(i01, i03i)
1060  END DO
1061  END DO
1062 
1063  END SUBROUTINE inserthalf
1064 
1065 ! **************************************************************************************************
1066 !> \brief (Based on suitable modifications of S.Goedecker routines)
1067 !> Calculates the FFT of the distributed kernel
1068 !> \param n1 logical dimension of the transform.
1069 !> \param n2 logical dimension of the transform.
1070 !> \param n3 logical dimension of the transform.
1071 !> \param nd1 Dimensions of work arrays
1072 !> \param nd2 Dimensions of work arrays
1073 !> \param nd3 Dimensions of work arrays
1074 !> \param nk1 ...
1075 !> \param nk2 ...
1076 !> \param nk3 ...
1077 !> \param nproc number of processors used as returned by MPI_COMM_SIZE
1078 !> \param iproc [0:nproc-1] number of processor as returned by MPI_COMM_RANK
1079 !> \param zf Real kernel (input)
1080 !> zf(i1,i2,i3)
1081 !> \param zr Distributed Kernel FFT
1082 !> zr(2,i1,i2,i3)
1083 !> \param mpi_group ...
1084 !> \date February 2006
1085 !> \par Restrictions
1086 !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1087 !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1088 !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
1089 !> This file is distributed under the terms of the
1090 !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
1091 !> \author S. Goedecker, L. Genovese
1092 !> \note As transform lengths
1093 !> most products of the prime factors 2,3,5 are allowed.
1094 !> The detailed table with allowed transform lengths can
1095 !> be found in subroutine CTRIG
1096 ! **************************************************************************************************
1097  SUBROUTINE kernelfft(n1, n2, n3, nd1, nd2, nd3, nk1, nk2, nk3, nproc, iproc, zf, zr, mpi_group)
1098 
1099  INTEGER, INTENT(in) :: n1, n2, n3, nd1, nd2, nd3, nk1, nk2, &
1100  nk3, nproc, iproc
1101  REAL(kind=dp), &
1102  DIMENSION(n1/2+1, n3/2+1, nd2/nproc), &
1103  INTENT(in) :: zf
1104  REAL(kind=dp), DIMENSION(nk1, nk2, nk3/nproc), &
1105  INTENT(inout) :: zr
1106  TYPE(mp_comm_type), INTENT(in) :: mpi_group
1107 
1108  INTEGER, PARAMETER :: ncache_optimal = 8*1024
1109 
1110  INTEGER :: i, i1, i3, ic1, ic2, ic3, inzee, j, j2, &
1111  j2st, j3, jp2st, lot, lzt, ma, mb, &
1112  ncache, nfft
1113  INTEGER, ALLOCATABLE, DIMENSION(:) :: after1, after2, after3, before1, &
1114  before2, before3, now1, now2, now3
1115  REAL(kind=dp) :: twopion
1116  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cosinarr, trig1, trig2, trig3
1117  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zt, zw
1118  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: zmpi2
1119  REAL(kind=dp), ALLOCATABLE, &
1120  DIMENSION(:, :, :, :, :) :: zmpi1
1121 
1122 ! include "perfdata.inc"
1123 !work arrays for transpositions
1124 !work arrays for MPI
1125 !cache work array
1126 !FFT work arrays
1127 !Body
1128 !check input
1129 
1130  cpassert(nd1 .GE. n1)
1131  cpassert(nd2 .GE. n2)
1132  cpassert(nd3 .GE. n3/2 + 1)
1133  cpassert(mod(nd3, nproc) .EQ. 0)
1134  cpassert(mod(nd2, nproc) .EQ. 0)
1135  mark_used(nd1)
1136 
1137  !defining work arrays dimensions
1138  ncache = ncache_optimal
1139  IF (ncache <= max(n1, n2, n3/2)*4) ncache = max(n1, n2, n3/2)*4
1140  lzt = n2
1141  IF (mod(n2, 2) .EQ. 0) lzt = lzt + 1
1142  IF (mod(n2, 4) .EQ. 0) lzt = lzt + 1
1143 
1144  !Allocations
1145  ALLOCATE (trig1(2, ctrig_length))
1146  ALLOCATE (after1(7))
1147  ALLOCATE (now1(7))
1148  ALLOCATE (before1(7))
1149  ALLOCATE (trig2(2, ctrig_length))
1150  ALLOCATE (after2(7))
1151  ALLOCATE (now2(7))
1152  ALLOCATE (before2(7))
1153  ALLOCATE (trig3(2, ctrig_length))
1154  ALLOCATE (after3(7))
1155  ALLOCATE (now3(7))
1156  ALLOCATE (before3(7))
1157  ALLOCATE (zw(2, ncache/4, 2))
1158  ALLOCATE (zt(2, lzt, n1))
1159  ALLOCATE (zmpi2(2, n1, nd2/nproc, nd3))
1160  ALLOCATE (cosinarr(2, n3/2))
1161  IF (nproc .GT. 1) THEN
1162  ALLOCATE (zmpi1(2, n1, nd2/nproc, nd3/nproc, nproc))
1163  zmpi1 = 0.0_dp
1164  END IF
1165 
1166  zmpi2 = 0.0_dp
1167  !calculating the FFT work arrays (beware on the HalFFT in n3 dimension)
1168  CALL ctrig(n3/2, trig3, after3, before3, now3, 1, ic3)
1169  CALL ctrig(n1, trig1, after1, before1, now1, 1, ic1)
1170  CALL ctrig(n2, trig2, after2, before2, now2, 1, ic2)
1171 
1172  !Calculating array of phases for HalFFT decoding
1173  twopion = 8._dp*atan(1._dp)/real(n3, kind=dp)
1174  DO i3 = 1, n3/2
1175  cosinarr(1, i3) = cos(twopion*(i3 - 1))
1176  cosinarr(2, i3) = -sin(twopion*(i3 - 1))
1177  END DO
1178 
1179  !transform along z axis
1180 
1181  lot = ncache/(2*n3)
1182  cpassert(lot .GE. 1)
1183 
1184  DO j2 = 1, nd2/nproc
1185  !this condition ensures that we manage only the interesting part for the FFT
1186  IF (iproc*(nd2/nproc) + j2 .LE. n2) THEN
1187  DO i1 = 1, n1, lot
1188  ma = i1
1189  mb = min(i1 + (lot - 1), n1)
1190  nfft = mb - ma + 1
1191 
1192  !inserting real data into complex array of half length
1193  !input: I1,I3,J2,(Jp2)
1194 
1195  CALL inserthalf(n1, n3, lot, nfft, i1, zf(1, 1, j2), zw(1, 1, 1))
1196 
1197  !performing FFT
1198  inzee = 1
1199  DO i = 1, ic3
1200  CALL fftstp(lot, nfft, n3/2, lot, n3/2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1201  trig3, after3(i), now3(i), before3(i), 1)
1202  inzee = 3 - inzee
1203  END DO
1204  !output: I1,i3,J2,(Jp2)
1205 
1206  !unpacking FFT in order to restore correct result,
1207  !while exchanging components
1208  !input: I1,i3,J2,(Jp2)
1209  CALL scramble_unpack(i1, j2, lot, nfft, n1, n3, nd2, nproc, nd3, zw(1, 1, inzee), zmpi2, cosinarr)
1210  !output: I1,J2,i3,(Jp2)
1211  END DO
1212  END IF
1213  END DO
1214 
1215  !Interprocessor data transposition
1216  !input: I1,J2,j3,jp3,(Jp2)
1217  IF (nproc .GT. 1) THEN
1218  !communication scheduling
1219  CALL mpi_group%alltoall(zmpi2, &!2*n1*(nd2/nproc)*(nd3/nproc), &
1220  zmpi1, 2*n1*(nd2/nproc)*(nd3/nproc))
1221  ! output: I1,J2,j3,Jp2,(jp3)
1222  END IF
1223 
1224  DO j3 = 1, nd3/nproc
1225  !this condition ensures that we manage only the interesting part for the FFT
1226  IF (iproc*(nd3/nproc) + j3 .LE. n3/2 + 1) THEN
1227  jp2st = 1
1228  j2st = 1
1229 
1230  !transform along x axis
1231  lot = ncache/(4*n1)
1232  cpassert(lot .GE. 1)
1233 
1234  DO j = 1, n2, lot
1235  ma = j
1236  mb = min(j + (lot - 1), n2)
1237  nfft = mb - ma + 1
1238 
1239  !reverse ordering
1240  !input: I1,J2,j3,Jp2,(jp3)
1241  IF (nproc .EQ. 1) THEN
1242  CALL mpiswitch(j3, nfft, jp2st, j2st, lot, n1, nd2, nd3, nproc, zmpi2, zw(1, 1, 1))
1243  ELSE
1244  CALL mpiswitch(j3, nfft, jp2st, j2st, lot, n1, nd2, nd3, nproc, zmpi1, zw(1, 1, 1))
1245  END IF
1246  !output: J2,Jp2,I1,j3,(jp3)
1247 
1248  !performing FFT
1249  !input: I2,I1,j3,(jp3)
1250  inzee = 1
1251  DO i = 1, ic1 - 1
1252  CALL fftstp(lot, nfft, n1, lot, n1, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1253  trig1, after1(i), now1(i), before1(i), 1)
1254  inzee = 3 - inzee
1255  END DO
1256  !storing the last step into zt
1257  i = ic1
1258  CALL fftstp(lot, nfft, n1, lzt, n1, zw(1, 1, inzee), zt(1, j, 1), &
1259  trig1, after1(i), now1(i), before1(i), 1)
1260  !output: I2,i1,j3,(jp3)
1261  END DO
1262 
1263  !transform along y axis, and taking only the first half
1264  lot = ncache/(4*n2)
1265  cpassert(lot .GE. 1)
1266 
1267  DO j = 1, nk1, lot
1268  ma = j
1269  mb = min(j + (lot - 1), nk1)
1270  nfft = mb - ma + 1
1271 
1272  !reverse ordering
1273  !input: I2,i1,j3,(jp3)
1274  CALL switch(nfft, n2, lot, n1, lzt, zt(1, 1, j), zw(1, 1, 1))
1275  !output: i1,I2,j3,(jp3)
1276 
1277  !performing FFT
1278  !input: i1,I2,j3,(jp3)
1279  inzee = 1
1280  DO i = 1, ic2
1281  CALL fftstp(lot, nfft, n2, lot, n2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1282  trig2, after2(i), now2(i), before2(i), 1)
1283  inzee = 3 - inzee
1284  END DO
1285 
1286  CALL realcopy(lot, nfft, n2, nk1, nk2, zw(1, 1, inzee), zr(j, 1, j3))
1287 
1288  END DO
1289  !output: i1,i2,j3,(jp3)
1290  END IF
1291  END DO
1292 
1293  !De-allocations
1294  DEALLOCATE (trig1)
1295  DEALLOCATE (after1)
1296  DEALLOCATE (now1)
1297  DEALLOCATE (before1)
1298  DEALLOCATE (trig2)
1299  DEALLOCATE (after2)
1300  DEALLOCATE (now2)
1301  DEALLOCATE (before2)
1302  DEALLOCATE (trig3)
1303  DEALLOCATE (after3)
1304  DEALLOCATE (now3)
1305  DEALLOCATE (before3)
1306  DEALLOCATE (zmpi2)
1307  DEALLOCATE (zw)
1308  DEALLOCATE (zt)
1309  DEALLOCATE (cosinarr)
1310  IF (nproc .GT. 1) DEALLOCATE (zmpi1)
1311 
1312  END SUBROUTINE kernelfft
1313 
1314 ! **************************************************************************************************
1315 !> \brief ...
1316 !> \param lot ...
1317 !> \param nfft ...
1318 !> \param n2 ...
1319 !> \param nk1 ...
1320 !> \param nk2 ...
1321 !> \param zin ...
1322 !> \param zout ...
1323 ! **************************************************************************************************
1324  SUBROUTINE realcopy(lot, nfft, n2, nk1, nk2, zin, zout)
1325  INTEGER, INTENT(in) :: lot, nfft, n2, nk1, nk2
1326  REAL(kind=dp), DIMENSION(2, lot, n2), INTENT(in) :: zin
1327  REAL(kind=dp), DIMENSION(nk1, nk2), INTENT(inout) :: zout
1328 
1329  INTEGER :: i, j
1330 
1331  DO i = 1, nk2
1332  DO j = 1, nfft
1333  zout(j, i) = zin(1, j, i)
1334  END DO
1335  END DO
1336 
1337  END SUBROUTINE realcopy
1338 
1339 ! **************************************************************************************************
1340 !> \brief ...
1341 !> \param nfft ...
1342 !> \param n2 ...
1343 !> \param lot ...
1344 !> \param n1 ...
1345 !> \param lzt ...
1346 !> \param zt ...
1347 !> \param zw ...
1348 ! **************************************************************************************************
1349  SUBROUTINE switch(nfft, n2, lot, n1, lzt, zt, zw)
1350  INTEGER :: nfft, n2, lot, n1, lzt
1351  REAL(kind=dp) :: zt(2, lzt, n1), zw(2, lot, n2)
1352 
1353  INTEGER :: i, j
1354 
1355  DO 200, j = 1, nfft
1356  DO 100, i = 1, n2
1357  zw(1, j, i) = zt(1, i, j)
1358  zw(2, j, i) = zt(2, i, j)
1359 100 CONTINUE
1360 200 CONTINUE
1361  RETURN
1362  END SUBROUTINE switch
1363 
1364 ! **************************************************************************************************
1365 !> \brief ...
1366 !> \param j3 ...
1367 !> \param nfft ...
1368 !> \param Jp2st ...
1369 !> \param J2st ...
1370 !> \param lot ...
1371 !> \param n1 ...
1372 !> \param nd2 ...
1373 !> \param nd3 ...
1374 !> \param nproc ...
1375 !> \param zmpi1 ...
1376 !> \param zw ...
1377 ! **************************************************************************************************
1378  SUBROUTINE mpiswitch(j3, nfft, Jp2st, J2st, lot, n1, nd2, nd3, nproc, zmpi1, zw)
1379  INTEGER :: j3, nfft, jp2st, j2st, lot, n1, nd2, &
1380  nd3, nproc
1381  REAL(kind=dp) :: zmpi1(2, n1, nd2/nproc, nd3/nproc, nproc), zw(2, lot, n1)
1382 
1383  INTEGER :: i1, j2, jp2, mfft
1384 
1385  mfft = 0
1386  DO 300, jp2 = jp2st, nproc
1387  DO 200, j2 = j2st, nd2/nproc
1388  mfft = mfft + 1
1389  IF (mfft .GT. nfft) THEN
1390  jp2st = jp2
1391  j2st = j2
1392  RETURN
1393  END IF
1394  DO 100, i1 = 1, n1
1395  zw(1, mfft, i1) = zmpi1(1, i1, j2, j3, jp2)
1396  zw(2, mfft, i1) = zmpi1(2, i1, j2, j3, jp2)
1397 100 CONTINUE
1398 200 CONTINUE
1399  j2st = 1
1400 300 CONTINUE
1401  END SUBROUTINE mpiswitch
1402 
1403 ! **************************************************************************************************
1404 !> \brief ...
1405 !> \param p ...
1406 !> \param w ...
1407 !> \param urange ...
1408 !> \param drange ...
1409 !> \param acc ...
1410 ! **************************************************************************************************
1411  SUBROUTINE gequad(p, w, urange, drange, acc)
1412 !
1413  REAL(kind=dp) :: p(*), w(*), urange, drange, acc
1414 
1415 !
1416 !
1417 ! range [10^(-9),1] and accuracy ~10^(-8);
1418 !
1419 !
1420 
1421  p(1) = 4.96142640560223544e19_dp
1422  p(2) = 1.37454269147978052e19_dp
1423  p(3) = 7.58610013441204679e18_dp
1424  p(4) = 4.42040691347806996e18_dp
1425  p(5) = 2.61986077948367892e18_dp
1426  p(6) = 1.56320138155496681e18_dp
1427  p(7) = 9.35645215863028402e17_dp
1428  p(8) = 5.60962910452691703e17_dp
1429  p(9) = 3.3666225119686761e17_dp
1430  p(10) = 2.0218253197947866e17_dp
1431  p(11) = 1.21477756091902017e17_dp
1432  p(12) = 7.3012982513608503e16_dp
1433  p(13) = 4.38951893556421099e16_dp
1434  p(14) = 2.63949482512262325e16_dp
1435  p(15) = 1.58742054072786174e16_dp
1436  p(16) = 9.54806587737665531e15_dp
1437  p(17) = 5.74353712364571709e15_dp
1438  p(18) = 3.455214877389445e15_dp
1439  p(19) = 2.07871658520326804e15_dp
1440  p(20) = 1.25064667315629928e15_dp
1441  p(21) = 7.52469429541933745e14_dp
1442  p(22) = 4.5274603337253175e14_dp
1443  p(23) = 2.72414006900059548e14_dp
1444  p(24) = 1.63912168349216752e14_dp
1445  p(25) = 9.86275802590865738e13_dp
1446  p(26) = 5.93457701624974985e13_dp
1447  p(27) = 3.5709554322296296e13_dp
1448  p(28) = 2.14872890367310454e13_dp
1449  p(29) = 1.29294719957726902e13_dp
1450  p(30) = 7.78003375426361016e12_dp
1451  p(31) = 4.68148199759876704e12_dp
1452  p(32) = 2.8169955024829868e12_dp
1453  p(33) = 1.69507790481958464e12_dp
1454  p(34) = 1.01998486064607581e12_dp
1455  p(35) = 6.13759486539856459e11_dp
1456  p(36) = 3.69320183828682544e11_dp
1457  p(37) = 2.22232783898905102e11_dp
1458  p(38) = 1.33725247623668682e11_dp
1459  p(39) = 8.0467192739036288e10_dp
1460  p(40) = 4.84199582415144143e10_dp
1461  p(41) = 2.91360091170559564e10_dp
1462  p(42) = 1.75321747475309216e10_dp
1463  p(43) = 1.0549735552210995e10_dp
1464  p(44) = 6.34815321079006586e9_dp
1465  p(45) = 3.81991113733594231e9_dp
1466  p(46) = 2.29857747533101109e9_dp
1467  p(47) = 1.38313653595483694e9_dp
1468  p(48) = 8.32282908580025358e8_dp
1469  p(49) = 5.00814519374587467e8_dp
1470  p(50) = 3.01358090773319025e8_dp
1471  p(51) = 1.81337994217503535e8_dp
1472  p(52) = 1.09117589961086823e8_dp
1473  p(53) = 6.56599771718640323e7_dp
1474  p(54) = 3.95099693638497164e7_dp
1475  p(55) = 2.37745694710665991e7_dp
1476  p(56) = 1.43060135285912813e7_dp
1477  p(57) = 8.60844290313506695e6_dp
1478  p(58) = 5.18000974075383424e6_dp
1479  p(59) = 3.116998193057466e6_dp
1480  p(60) = 1.87560993870024029e6_dp
1481  p(61) = 1.12862197183979562e6_dp
1482  p(62) = 679132.441326077231_dp
1483  p(63) = 408658.421279877969_dp
1484  p(64) = 245904.473450669789_dp
1485  p(65) = 147969.568088321005_dp
1486  p(66) = 89038.612357311147_dp
1487  p(67) = 53577.7362552358895_dp
1488  p(68) = 32239.6513926914668_dp
1489  p(69) = 19399.7580852362791_dp
1490  p(70) = 11673.5323603058634_dp
1491  p(71) = 7024.38438577707758_dp
1492  p(72) = 4226.82479307685999_dp
1493  p(73) = 2543.43254175354295_dp
1494  p(74) = 1530.47486269122675_dp
1495  p(75) = 920.941785160749482_dp
1496  p(76) = 554.163803906291646_dp
1497  p(77) = 333.46029740785694_dp
1498  p(78) = 200.6550575335041_dp
1499  p(79) = 120.741366914147284_dp
1500  p(80) = 72.6544243200329916_dp
1501  p(81) = 43.7187810415471025_dp
1502  p(82) = 26.3071631447061043_dp
1503  p(83) = 15.8299486353816329_dp
1504  p(84) = 9.52493152341244004_dp
1505  p(85) = 5.72200417067776041_dp
1506  p(86) = 3.36242234070940928_dp
1507  p(87) = 1.75371394604499472_dp
1508  p(88) = 0.64705932650658966_dp
1509  p(89) = 0.072765905943708247_dp
1510  !
1511  w(1) = 47.67445484528304247e10_dp
1512  w(2) = 11.37485774750442175e9_dp
1513  w(3) = 78.64340976880190239e8_dp
1514  w(4) = 46.27335788759590498e8_dp
1515  w(5) = 24.7380464827152951e8_dp
1516  w(6) = 13.62904116438987719e8_dp
1517  w(7) = 92.79560029045882433e8_dp
1518  w(8) = 52.15931216254660251e8_dp
1519  w(9) = 31.67018011061666244e8_dp
1520  w(10) = 1.29291036801493046e8_dp
1521  w(11) = 1.00139319988015862e8_dp
1522  w(12) = 7.75892350510188341e7_dp
1523  w(13) = 6.01333567950731271e7_dp
1524  w(14) = 4.66141178654796875e7_dp
1525  w(15) = 3.61398903394911448e7_dp
1526  w(16) = 2.80225846672956389e7_dp
1527  w(17) = 2.1730509180930247e7_dp
1528  w(18) = 1.68524482625876965e7_dp
1529  w(19) = 1.30701489345870338e7_dp
1530  w(20) = 1.01371784832269282e7_dp
1531  w(21) = 7.86264116300379329e6_dp
1532  w(22) = 6.09861667912273717e6_dp
1533  w(23) = 4.73045784039455683e6_dp
1534  w(24) = 3.66928949951594161e6_dp
1535  w(25) = 2.8462050836230259e6_dp
1536  w(26) = 2.20777394798527011e6_dp
1537  w(27) = 1.71256191589205524e6_dp
1538  w(28) = 1.32843556197737076e6_dp
1539  w(29) = 1.0304731275955989e6_dp
1540  w(30) = 799345.206572271448_dp
1541  w(31) = 620059.354143595343_dp
1542  w(32) = 480986.704107449333_dp
1543  w(33) = 373107.167700228515_dp
1544  w(34) = 289424.08337412132_dp
1545  w(35) = 224510.248231581788_dp
1546  w(36) = 174155.825690028966_dp
1547  w(37) = 135095.256919654065_dp
1548  w(38) = 104795.442776800312_dp
1549  w(39) = 81291.4458222430418_dp
1550  w(40) = 63059.0493649328682_dp
1551  w(41) = 48915.9040455329689_dp
1552  w(42) = 37944.8484018048756_dp
1553  w(43) = 29434.4290473253969_dp
1554  w(44) = 22832.7622054490044_dp
1555  w(45) = 17711.743950151233_dp
1556  w(46) = 13739.287867104177_dp
1557  w(47) = 10657.7895710752585_dp
1558  w(48) = 8267.42141053961834_dp
1559  w(49) = 6413.17397520136448_dp
1560  w(50) = 4974.80402838654277_dp
1561  w(51) = 3859.03698188553047_dp
1562  w(52) = 2993.51824493299154_dp
1563  w(53) = 2322.1211966811754_dp
1564  w(54) = 1801.30750964719641_dp
1565  w(55) = 1397.30379659817038_dp
1566  w(56) = 1083.91149143250697_dp
1567  w(57) = 840.807939169209188_dp
1568  w(58) = 652.228524366749422_dp
1569  w(59) = 505.944376983506128_dp
1570  w(60) = 392.469362317941064_dp
1571  w(61) = 304.444930257324312_dp
1572  w(62) = 236.162932842453601_dp
1573  w(63) = 183.195466078603525_dp
1574  w(64) = 142.107732186551471_dp
1575  w(65) = 110.23530215723992_dp
1576  w(66) = 85.5113346705382257_dp
1577  w(67) = 66.3325469806696621_dp
1578  w(68) = 51.4552463353841373_dp
1579  w(69) = 39.9146798429449273_dp
1580  w(70) = 30.9624728409162095_dp
1581  w(71) = 24.018098812215013_dp
1582  w(72) = 18.6312338024296588_dp
1583  w(73) = 14.4525541233150501_dp
1584  w(74) = 11.2110836519105938_dp
1585  w(75) = 8.69662175848497178_dp
1586  w(76) = 6.74611236165731961_dp
1587  w(77) = 5.23307018057529994_dp
1588  w(78) = 4.05937850501539556_dp
1589  w(79) = 3.14892659076635714_dp
1590  w(80) = 2.44267408211071604_dp
1591  w(81) = 1.89482240522855261_dp
1592  w(82) = 1.46984505907050079_dp
1593  w(83) = 1.14019261330527007_dp
1594  w(84) = 0.884791217422925293_dp
1595  w(85) = 0.692686387080616483_dp
1596  w(86) = 0.585244576897023282_dp
1597  w(87) = 0.576182522545327589_dp
1598  w(88) = 0.596688817388997178_dp
1599  w(89) = 0.607879901151108771_dp
1600  !
1601  !
1602  urange = 1._dp
1603  drange = 1e-08_dp
1604  acc = 1e-08_dp
1605  !
1606  RETURN
1607  END SUBROUTINE gequad
1608 
1609  END MODULE ps_wavelet_kernel
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
for
for(int lxp=0;lxp<=lp;lxp++)
Definition: grid_cpu_collint.h:66
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
ps_wavelet_base
Creates the wavelet kernel for the wavelet based poisson solver.
Definition: ps_wavelet_base.F:12
ps_wavelet_base::scramble_unpack
subroutine, public scramble_unpack(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zw, zmpi2, cosinarr)
(Based on suitable modifications of S.Goedecker routines) Assign the correct planes to the work array...
Definition: ps_wavelet_base.F:1440
ps_wavelet_fft3d
Definition: ps_wavelet_fft3d.F:9
ps_wavelet_fft3d::ctrig_length
integer, parameter, public ctrig_length
Definition: ps_wavelet_fft3d.F:20
ps_wavelet_fft3d::fftstp
subroutine, public fftstp(mm, nfft, m, nn, n, zin, zout, trig, after, now, before, isign)
...
Definition: ps_wavelet_fft3d.F:280
ps_wavelet_fft3d::ctrig
subroutine, public ctrig(n, trig, after, before, now, isign, ic)
...
Definition: ps_wavelet_fft3d.F:105
ps_wavelet_kernel
Creates the wavelet kernel for the wavelet based poisson solver.
Definition: ps_wavelet_kernel.F:12
ps_wavelet_kernel::createkernel
subroutine, public createkernel(geocode, n01, n02, n03, hx, hy, hz, itype_scf, iproc, nproc, kernel, mpi_group)
Allocate a pointer which corresponds to the zero-padded FFT slice needed for calculating the convolut...
Definition: ps_wavelet_kernel.F:80
ps_wavelet_scaling_function
Creates the wavelet kernel for the wavelet based poisson solver.
Definition: ps_wavelet_scaling_function.F:12
ps_wavelet_scaling_function::scf_recursion
subroutine, public scf_recursion(itype, n_iter, n_range, kernel_scf, kern_1_scf)
Do iterations to go from p0gauss to pgauss order interpolating scaling function.
Definition: ps_wavelet_scaling_function.F:151
ps_wavelet_scaling_function::scaling_function
subroutine, public scaling_function(itype, nd, nrange, a, x)
Calculate the values of a scaling function in real uniform grid.
Definition: ps_wavelet_scaling_function.F:35
ps_wavelet_util
Performs a wavelet based solution of the Poisson equation.
Definition: ps_wavelet_util.F:12
ps_wavelet_util::s_fft_dimensions
subroutine, public s_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
Calculate four sets of dimension needed for the calculation of the convolution for the surface system...
Definition: ps_wavelet_util.F:334
ps_wavelet_util::f_fft_dimensions
subroutine, public f_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
Calculate four sets of dimension needed for the calculation of the zero-padded convolution.
Definition: ps_wavelet_util.F:441