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