Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
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 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 :
37 : CONTAINS
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 856 : 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 856 : hgrid = MAX(hx, hy, hz)
95 :
96 856 : IF (geocode == 'P') THEN
97 :
98 : CALL F_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, &
99 326 : md1, md2, md3, nd1, nd2, nd3, nproc)
100 :
101 326 : ALLOCATE (kernel(1))
102 326 : nlimd = n2
103 652 : nlimk = 0
104 :
105 530 : ELSE IF (geocode == 'S') THEN
106 :
107 : CALL S_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, &
108 6 : md1, md2, md3, nd1, nd2, nd3, nproc)
109 :
110 18 : 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 6 : itype_scf, kernel, iproc, nproc, mpi_group)
116 : !last plane calculated for the density and the kernel
117 :
118 6 : nlimd = n2
119 12 : nlimk = n3/2 + 1
120 524 : 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 524 : md1, md2, md3, nd1, nd2, nd3, nproc)
126 1572 : 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 524 : itype_scf, iproc, nproc, kernel, mpi_group)
131 :
132 : !last plane calculated for the density and the kernel
133 524 : nlimd = n2/2
134 1048 : nlimk = n3/2 + 1
135 :
136 : ELSE
137 :
138 0 : 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 856 : 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 6 : SUBROUTINE Surfaces_Kernel(n1, n2, n3, m3, nker1, nker2, nker3, h1, h2, h3, &
222 6 : 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 6 : 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 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: kernel_scf, x_scf, y_scf
242 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: btrig, cossinarr
243 6 : 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 136590 : karray = 0.0_dp
254 486 : cpol(:, :) = 1._dp
255 :
256 6 : cpol(1, 2) = .25_dp
257 :
258 6 : cpol(1, 3) = 1._dp/3._dp
259 :
260 6 : cpol(1, 4) = 7._dp/12._dp
261 6 : cpol(2, 4) = 8._dp/3._dp
262 :
263 6 : cpol(1, 5) = 19._dp/50._dp
264 6 : cpol(2, 5) = 3._dp/2._dp
265 :
266 6 : cpol(1, 6) = 41._dp/272._dp
267 6 : cpol(2, 6) = 27._dp/34._dp
268 6 : cpol(3, 6) = 27._dp/272._dp
269 :
270 6 : cpol(1, 7) = 751._dp/2989._dp
271 6 : cpol(2, 7) = 73._dp/61._dp
272 6 : cpol(3, 7) = 27._dp/61._dp
273 :
274 6 : cpol(1, 8) = -989._dp/4540._dp
275 6 : cpol(2, 8) = -1472._dp/1135._dp
276 6 : cpol(3, 8) = 232._dp/1135._dp
277 6 : cpol(4, 8) = -2624._dp/1135._dp
278 :
279 : !renormalize values
280 6 : cpol(1, 1) = .5_dp*cpol(1, 1)
281 18 : cpol(1:2, 2) = 2._dp/3._dp*cpol(1:2, 2)
282 18 : cpol(1:2, 3) = 3._dp/8._dp*cpol(1:2, 3)
283 24 : cpol(1:3, 4) = 2._dp/15._dp*cpol(1:3, 4)
284 24 : cpol(1:3, 5) = 25._dp/144._dp*cpol(1:3, 5)
285 30 : cpol(1:4, 6) = 34._dp/105._dp*cpol(1:4, 6)
286 30 : cpol(1:4, 7) = 2989._dp/17280._dp*cpol(1:4, 7)
287 36 : cpol(1:5, 8) = -454._dp/2835._dp*cpol(1:5, 8)
288 :
289 : !assign the complete values
290 6 : cpol(2, 1) = cpol(1, 1)
291 :
292 6 : cpol(3, 2) = cpol(1, 2)
293 :
294 6 : cpol(3, 3) = cpol(2, 3)
295 6 : cpol(4, 3) = cpol(1, 3)
296 :
297 6 : cpol(4, 4) = cpol(2, 4)
298 6 : cpol(5, 4) = cpol(1, 4)
299 :
300 6 : cpol(4, 5) = cpol(3, 5)
301 6 : cpol(5, 5) = cpol(2, 5)
302 6 : cpol(6, 5) = cpol(1, 5)
303 :
304 6 : cpol(5, 6) = cpol(3, 6)
305 6 : cpol(6, 6) = cpol(2, 6)
306 6 : cpol(7, 6) = cpol(1, 6)
307 :
308 6 : cpol(5, 7) = cpol(4, 7)
309 6 : cpol(6, 7) = cpol(3, 7)
310 6 : cpol(7, 7) = cpol(2, 7)
311 6 : cpol(8, 7) = cpol(1, 7)
312 :
313 6 : cpol(6, 8) = cpol(4, 8)
314 6 : cpol(7, 8) = cpol(3, 8)
315 6 : cpol(8, 8) = cpol(2, 8)
316 6 : cpol(9, 8) = cpol(1, 8)
317 :
318 : !Number of integration points : 2*itype_scf*n_points
319 6 : n_scf = 2*itype_scf*n_points
320 : !Allocations
321 18 : ALLOCATE (x_scf(0:n_scf))
322 12 : ALLOCATE (y_scf(0:n_scf))
323 :
324 : !Build the scaling function
325 6 : CALL scaling_function(itype_scf, n_scf, n_range, x_scf, y_scf)
326 : !Step grid for the integration
327 6 : 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 6 : n_cell = m3
330 6 : n_range = MAX(n_cell, n_range)
331 :
332 : !Allocations
333 6 : ncache = ncache_optimal
334 : !the HalFFT must be performed only in the third dimension,
335 : !and nker3=n3/2+1, hence
336 6 : IF (ncache <= (nker3 - 1)*4) ncache = nker3 - 1*4
337 :
338 : !enlarge the second dimension of the kernel to be compatible with nproc
339 6 : nact2 = nker2
340 0 : enlarge_ydim: DO
341 6 : IF (nproc*(nact2/nproc) /= nact2) THEN
342 0 : 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 30 : ALLOCATE (kernel(nker1, nact2/nproc, nker3))
350 36 : ALLOCATE (kernel_mpi(nker1, nact2/nproc, nker3/nproc, nproc))
351 18 : ALLOCATE (kernel_scf(n_range))
352 24 : ALLOCATE (halfft_cache(2, ncache/4, 2))
353 18 : ALLOCATE (cossinarr(2, n3/2 - 1))
354 6 : ALLOCATE (btrig(2, ctrig_length))
355 6 : ALLOCATE (after(7))
356 6 : ALLOCATE (now(7))
357 6 : ALLOCATE (before(7))
358 :
359 : !arrays for the halFFT
360 6 : CALL ctrig(n3/2, btrig, after, before, now, 1, ic)
361 :
362 : !build the phases for the HalFFT reconstruction
363 6 : pion = 2._dp*pi/REAL(n3, KIND=dp)
364 324 : DO i3 = 2, n3/2
365 318 : x = REAL(i3 - 1, KIND=dp)*pion
366 318 : cossinarr(1, i3 - 1) = COS(x)
367 324 : 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 136758 : kernel = HUGE(0._dp)
372 136770 : 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 6 : num_of_mus = ncache/(2*n3)
377 :
378 6 : diff = 0._dp
379 : !order of the polynomial to be used for integration (must be a power of two)
380 6 : 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 6 : j2st = iproc*(nact2/nproc)
385 6 : j2nd = MIN((iproc + 1)*(nact2/nproc), n2/2 + 1)
386 :
387 72 : DO ind2 = (n1/2 + 1)*j2st + 1, (n1/2 + 1)*j2nd, num_of_mus
388 66 : istart = ind2
389 66 : iend = MIN(ind2 + (num_of_mus - 1), (n1/2 + 1)*j2nd)
390 66 : nfft = iend - istart + 1
391 66 : shift = 0
392 :
393 : !initialization of the interesting part of the cache array
394 66 : halfft_cache(:, :, :) = 0._dp
395 :
396 66 : IF (istart == 1) THEN
397 : !i2=1
398 3 : shift = 1
399 :
400 : CALL calculates_green_opt_muzero(n_range, n_scf, ipolyord, x_scf, y_scf, &
401 3 : cpol(1, ipolyord), dx, kernel_scf)
402 :
403 : !copy of the first zero value
404 3 : halfft_cache(1, 1, 1) = 0._dp
405 :
406 165 : DO i3 = 1, m3
407 :
408 162 : value = 0.5_dp*h3*kernel_scf(i3)
409 : !index in where to copy the value of the kernel
410 162 : CALL indices(ireim, num_of_mus, n3/2 + i3, 1, ind1)
411 : !index in where to copy the symmetric value
412 162 : CALL indices(jreim, num_of_mus, n3/2 + 2 - i3, 1, jnd1)
413 162 : halfft_cache(ireim, ind1, 1) = value
414 165 : halfft_cache(jreim, jnd1, 1) = value
415 :
416 : END DO
417 :
418 : END IF
419 :
420 2415 : 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 2349 : i1 = MOD(imu, n1/2 + 1)
429 2349 : IF (i1 == 0) i1 = n1/2 + 1
430 2349 : i2 = (imu - i1)/(n1/2 + 1) + 1
431 2349 : ponx = REAL(i1 - 1, KIND=dp)/REAL(n1, KIND=dp)
432 2349 : pony = REAL(i2 - 1, KIND=dp)/REAL(n2, KIND=dp)
433 :
434 2349 : 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 2349 : 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 2349 : halfft_cache(1, imu - istart + 1, 1) = 0._dp
443 129261 : DO i3 = 1, m3
444 126846 : 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 126846 : CALL indices(ireim, num_of_mus, n3/2 + i3, imu - istart + 1, ind1)
448 : !index in where to copy the symmetric value
449 126846 : CALL indices(jreim, num_of_mus, n3/2 + 2 - i3, imu - istart + 1, jnd1)
450 126846 : halfft_cache(ireim, ind1, 1) = value
451 129195 : 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 66 : inzee = 1
458 264 : 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 198 : btrig, after(i), now(i), before(i), 1)
462 264 : inzee = 3 - inzee
463 : END DO
464 : !assign the values of the FFT array
465 : !and compare with the good results
466 2424 : DO imu = istart, iend
467 :
468 : !corresponding value of i1 and i2
469 2352 : i1 = MOD(imu, n1/2 + 1)
470 2352 : IF (i1 == 0) i1 = n1/2 + 1
471 2352 : i2 = (imu - i1)/(n1/2 + 1) + 1
472 :
473 2352 : j2 = i2 - j2st
474 :
475 2352 : a = halfft_cache(1, imu - istart + 1, inzee)
476 2352 : b = halfft_cache(2, imu - istart + 1, inzee)
477 2352 : kernel(i1, j2, 1) = a + b
478 2352 : kernel(i1, j2, n3/2 + 1) = a - b
479 :
480 127074 : DO i3 = 2, n3/2
481 124656 : ind1 = imu - istart + 1 + num_of_mus*(i3 - 1)
482 124656 : jnd1 = imu - istart + 1 + num_of_mus*(n3/2 + 2 - i3 - 1)
483 124656 : cp = cossinarr(1, i3 - 1)
484 124656 : sp = cossinarr(2, i3 - 1)
485 124656 : a = halfft_cache(1, ind1, inzee)
486 124656 : b = halfft_cache(2, ind1, inzee)
487 124656 : c = halfft_cache(1, jnd1, inzee)
488 124656 : d = halfft_cache(2, jnd1, inzee)
489 124656 : feR = .5_dp*(a + c)
490 124656 : feI = .5_dp*(b - d)
491 124656 : foR = .5_dp*(a - c)
492 124656 : foI = .5_dp*(b + d)
493 124656 : fR = feR + cp*foI - sp*foR
494 127008 : 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 6 : IF (nproc > 1) THEN
502 : CALL mpi_group%alltoall(kernel, &!nker1*(nact2/nproc)*(nker3/nproc), &
503 6 : kernel_mpi, nker1*(nact2/nproc)*(nker3/nproc))
504 :
505 18 : DO jp2 = 1, nproc
506 354 : DO i3 = 1, nker3/nproc
507 5052 : DO i2 = 1, nact2/nproc
508 4704 : j2 = i2 + (jp2 - 1)*(nact2/nproc)
509 5040 : IF (j2 <= nker2) THEN
510 136416 : DO i1 = 1, nker1
511 : karray(i1, j2, i3) = &
512 136416 : 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 0 : karray(1:nker1, 1:nker2, 1:nker3) = kernel(1:nker1, 1:nker2, 1:nker3)
521 : END IF
522 :
523 : !De-allocations
524 6 : DEALLOCATE (kernel)
525 6 : DEALLOCATE (kernel_mpi)
526 6 : DEALLOCATE (btrig)
527 6 : DEALLOCATE (after)
528 6 : DEALLOCATE (now)
529 6 : DEALLOCATE (before)
530 6 : DEALLOCATE (halfft_cache)
531 6 : DEALLOCATE (kernel_scf)
532 6 : DEALLOCATE (x_scf)
533 6 : DEALLOCATE (y_scf)
534 :
535 12 : 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 2349 : 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 2349 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: green, green1
564 :
565 129195 : g_mu = 0.0_dp
566 : !We calculate the number of iterations to go from mu0 to mu0_ref
567 2349 : IF (mu <= mu_max) THEN
568 9 : n_iter = 0
569 9 : mu0 = mu
570 : ELSE
571 2340 : n_iter = 1
572 6888 : loop_iter: DO
573 9228 : ratio = REAL(2**n_iter, KIND=dp)
574 9228 : mu0 = mu/ratio
575 9228 : IF (mu0 <= mu_max) THEN
576 : EXIT loop_iter
577 : END IF
578 6888 : 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 2349 : nrec = 2**n_iter*n
584 :
585 7047 : ALLOCATE (green(-nrec:nrec))
586 :
587 : !initialization of the branching value
588 1956717 : ikern = 0
589 1956717 : izero = 0
590 1954368 : initialization: DO
591 1956717 : IF (xval(izero) >= REAL(ikern, KIND=dp) .OR. izero == n_scf) EXIT initialization
592 1954368 : izero = izero + 1
593 : END DO initialization
594 2349 : green = 0._dp
595 : !now perform the interpolation in right direction
596 2349 : ivalue = izero
597 2349 : gright = 0._dp
598 279531 : loop_right: DO
599 281880 : IF (ivalue >= n_scf - intorder - 1) EXIT loop_right
600 2795310 : DO i = 1, intorder + 1
601 2515779 : x = xval(ivalue) - REAL(ikern, KIND=dp)
602 2515779 : f = yval(ivalue)*EXP(-mu0*x)
603 2515779 : filter = intorder*c(i)
604 2515779 : gright = gright + filter*f
605 2795310 : ivalue = ivalue + 1
606 : END DO
607 279531 : ivalue = ivalue - 1
608 : END DO loop_right
609 2349 : iend = n_scf - ivalue
610 21141 : DO i = 1, iend
611 18792 : x = xval(ivalue) - REAL(ikern, KIND=dp)
612 18792 : f = yval(ivalue)*EXP(-mu0*x)
613 18792 : filter = intorder*c(i)
614 18792 : gright = gright + filter*f
615 21141 : ivalue = ivalue + 1
616 : END DO
617 2349 : gright = hres*gright
618 :
619 : !the scaling function is symmetric, so the same for the other direction
620 2349 : gleft = gright
621 :
622 2349 : green(ikern) = gleft + gright
623 :
624 : !now the loop until the last value
625 758880 : DO ikern = 1, nrec
626 1005372 : gltmp = 0._dp
627 1005372 : grtmp = 0._dp
628 1005372 : ivalue = izero
629 1005372 : x0 = xval(izero)
630 246645 : loop_integration: DO
631 1005372 : IF (izero == n_scf) EXIT loop_integration
632 2783565 : DO i = 1, intorder + 1
633 2536920 : x = xval(ivalue)
634 2536920 : fl = yval(ivalue)*EXP(mu0*x)
635 2536920 : fr = yval(ivalue)*EXP(-mu0*x)
636 2536920 : filter = intorder*c(i)
637 2536920 : gltmp = gltmp + filter*fl
638 2536920 : grtmp = grtmp + filter*fr
639 2536920 : ivalue = ivalue + 1
640 2536920 : IF (xval(izero) >= REAL(ikern, KIND=dp) .OR. izero == n_scf) THEN
641 : x1 = xval(izero)
642 : EXIT loop_integration
643 : END IF
644 2748330 : izero = izero + 1
645 : END DO
646 246645 : ivalue = ivalue - 1
647 246645 : izero = izero - 1
648 : END DO loop_integration
649 758727 : gleft = EXP(-mu0)*(gleft + hres*EXP(-mu0*REAL(ikern - 1, KIND=dp))*gltmp)
650 758727 : IF (izero == n_scf) THEN
651 : gright = 0._dp
652 : ELSE
653 32886 : gright = EXP(mu0)*(gright - hres*EXP(mu0*REAL(ikern - 1, KIND=dp))*grtmp)
654 : END IF
655 758727 : green(ikern) = gleft + gright
656 758727 : green(-ikern) = gleft + gright
657 758880 : IF (ABS(green(ikern)) <= 1.e-20_dp) THEN
658 2196 : nrec = ikern
659 2196 : 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 7047 : ALLOCATE (green1(-nrec:nrec))
665 :
666 : !Start the iteration to go from mu0 to mu
667 2349 : CALL scf_recursion(itype_scf, n_iter, nrec, green(-nrec), green1(-nrec))
668 :
669 129195 : DO i = 1, MIN(n, nrec)
670 129195 : g_mu(i) = green(i - 1)
671 : END DO
672 2349 : DO i = MIN(n, nrec) + 1, n
673 2349 : g_mu(i) = 0._dp
674 : END DO
675 :
676 2349 : DEALLOCATE (green, green1)
677 :
678 2349 : 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 3 : 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 3 : ikern = 0
704 3 : izero = 0
705 2496 : initialization: DO
706 2499 : IF (xval(izero) >= REAL(ikern, KIND=dp) .OR. izero == n_scf) EXIT initialization
707 2496 : izero = izero + 1
708 : END DO initialization
709 165 : green = 0._dp
710 : !first case, ikern=0
711 : !now perform the interpolation in right direction
712 : ivalue = izero
713 : gr1 = 0._dp
714 357 : loop_right: DO
715 360 : IF (ivalue >= n_scf - intorder - 1) EXIT loop_right
716 3570 : DO i = 1, intorder + 1
717 3213 : x = xval(ivalue)
718 3213 : y = yval(ivalue)
719 3213 : filter = intorder*c(i)
720 3213 : gr1 = gr1 + filter*x*y
721 3570 : ivalue = ivalue + 1
722 : END DO
723 357 : ivalue = ivalue - 1
724 : END DO loop_right
725 3 : iend = n_scf - ivalue
726 27 : DO i = 1, iend
727 24 : x = xval(ivalue)
728 24 : y = yval(ivalue)
729 24 : filter = intorder*c(i)
730 24 : gr1 = gr1 + filter*x*y
731 27 : ivalue = ivalue + 1
732 : END DO
733 3 : gr1 = hres*gr1
734 : !the scaling function is symmetric
735 3 : gl1 = -gr1
736 3 : gl0 = 0.5_dp
737 3 : gr0 = 0.5_dp
738 :
739 3 : green(1) = 2._dp*gr1
740 :
741 : !now the loop until the last value
742 162 : DO ikern = 1, n - 1
743 : c0 = 0._dp
744 : c1 = 0._dp
745 : ivalue = izero
746 315 : loop_integration: DO
747 474 : IF (izero == n_scf) EXIT loop_integration
748 3555 : DO i = 1, intorder + 1
749 3240 : x = xval(ivalue)
750 3240 : y = yval(ivalue)
751 3240 : filter = intorder*c(i)
752 3240 : c0 = c0 + filter*y
753 3240 : c1 = c1 + filter*y*x
754 3240 : ivalue = ivalue + 1
755 3240 : IF (xval(izero) >= REAL(ikern, KIND=dp) .OR. izero == n_scf) THEN
756 : EXIT loop_integration
757 : END IF
758 3510 : izero = izero + 1
759 : END DO
760 315 : ivalue = ivalue - 1
761 360 : izero = izero - 1
762 : END DO loop_integration
763 159 : c0 = hres*c0
764 159 : c1 = hres*c1
765 :
766 159 : gl0 = gl0 + c0
767 159 : gl1 = gl1 + c1
768 159 : gr0 = gr0 - c0
769 159 : gr1 = gr1 - c1
770 : !general case
771 162 : 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 3 : 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 254016 : 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 254016 : var_realimag = 2 - MOD(intrn, 2)
796 : !actual index over half the length
797 :
798 254016 : i = (intrn + 1)/2
799 : !check
800 254016 : IF (2*(i - 1) + var_realimag /= intrn) THEN
801 0 : PRINT *, 'error, index=', intrn, 'var_realimag=', var_realimag, 'i=', i
802 : END IF
803 : !complete index to be assigned
804 254016 : index = extrn + nelem*(i - 1)
805 :
806 254016 : 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 524 : SUBROUTINE Free_Kernel(n01, n02, n03, nfft1, nfft2, nfft3, n1k, n2k, n3k, &
836 524 : 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 524 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: kern_1_scf, kernel_scf, x_scf, y_scf
857 524 : 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 524 : n_scf = 2*itype_scf*n_points
867 : !Set karray
868 37584312 : 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 524 : n1h = nfft1/2
874 524 : n2h = nfft2/2
875 524 : n3h = nfft3/2
876 :
877 : !Auxiliary dimensions only for building the FFT part
878 524 : nker1 = nfft1
879 524 : nker2 = nfft2
880 524 : nker3 = nfft3/2 + 1
881 :
882 : !adjusting the last two dimensions to be multiples of nproc
883 0 : DO
884 524 : IF (MODULO(nker2, nproc) == 0) EXIT
885 0 : nker2 = nker2 + 1
886 : END DO
887 360 : DO
888 884 : IF (MODULO(nker3, nproc) == 0) EXIT
889 360 : nker3 = nker3 + 1
890 : END DO
891 :
892 : !this will be the array of the kernel in the real space
893 3144 : 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 524 : istart = iproc*nker2/nproc + 1
899 524 : iend = MIN((iproc + 1)*nker2/nproc, n2h + n03)
900 :
901 524 : istart1 = istart
902 524 : IF (iproc == 0) istart1 = n2h - n03 + 2
903 :
904 : !Allocations
905 1572 : ALLOCATE (x_scf(0:n_scf))
906 1048 : ALLOCATE (y_scf(0:n_scf))
907 :
908 : !Build the scaling function
909 524 : CALL scaling_function(itype_scf, n_scf, n_range, x_scf, y_scf)
910 : !Step grid for the integration
911 524 : 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 524 : n_cell = MAX(n01, n02, n03)
914 524 : n_range = MAX(n_cell, n_range)
915 :
916 : !Allocations
917 1572 : ALLOCATE (kernel_scf(-n_range:n_range))
918 1048 : ALLOCATE (kern_1_scf(-n_range:n_range))
919 :
920 : !Lengthes of the box (use FFT dimension)
921 524 : a1 = hgrid*REAL(n01, KIND=dp)
922 524 : a2 = hgrid*REAL(n02, KIND=dp)
923 524 : a3 = hgrid*REAL(n03, KIND=dp)
924 :
925 2670616 : x_scf(:) = hgrid*x_scf(:)
926 2670616 : y_scf(:) = 1._dp/hgrid*y_scf(:)
927 524 : dx = hgrid*dx
928 : !To have a correct integration
929 524 : p0_cell = p0_ref/(hgrid*hgrid)
930 :
931 : !Initialization of the gaussian (Beylkin)
932 524 : 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 524 : a_range = SQRT(a1*a1 + a2*a2 + a3*a3)
937 524 : factor = 1._dp/a_range
938 : !factor2 = factor*factor
939 524 : factor2 = 1._dp/(a1*a1 + a2*a2 + a3*a3)
940 47160 : DO i_gauss = 1, n_gauss
941 47160 : p_gauss(i_gauss) = factor2*p_gauss(i_gauss)
942 : END DO
943 47160 : DO i_gauss = 1, n_gauss
944 47160 : w_gauss(i_gauss) = factor*w_gauss(i_gauss)
945 : END DO
946 :
947 524 : kp(:, :, :) = 0._dp
948 : !Use in this order (better for accuracy).
949 47160 : loop_gauss: DO i_gauss = n_gauss, 1, -1
950 : !Gaussian
951 46636 : pgauss = p_gauss(i_gauss)
952 :
953 : !We calculate the number of iterations to go from pgauss to p0_ref
954 46636 : n_iter = NINT((LOG(pgauss) - LOG(p0_cell))/LOG(4._dp))
955 46636 : IF (n_iter <= 0) THEN
956 9944 : n_iter = 0
957 9944 : p0gauss = pgauss
958 : ELSE
959 36692 : 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 46636 : kernel_scf(:) = 0._dp
965 1471904 : DO i_kern = 0, n_range
966 : kern = 0._dp
967 7492650776 : DO i = 0, n_scf
968 7491182988 : absci = x_scf(i) - REAL(i_kern, KIND=dp)*hgrid
969 7491182988 : absci = absci*absci
970 7492650776 : kern = kern + y_scf(i)*EXP(-p0gauss*absci)*dx
971 : END DO
972 1467788 : kernel_scf(i_kern) = kern
973 1467788 : kernel_scf(-i_kern) = kern
974 1471904 : 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 46636 : 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 2468138 : DO i3 = istart1, iend
986 2420978 : 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 2420978 : 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 0 : "(see ps_wavelet_fft3d.F) exactly.")
998 : END IF
999 114386983 : DO i2 = 1, n02
1000 111919369 : i02 = i2 - 1
1001 6051005158 : DO i1 = 1, n01
1002 5936664811 : i01 = i1 - 1
1003 : kp(i1, i2, i3 - istart + 1) = kp(i1, i2, i3 - istart + 1) + w_gauss(i_gauss)* &
1004 6048584180 : 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 524 : DEALLOCATE (kernel_scf)
1013 524 : DEALLOCATE (kern_1_scf)
1014 524 : DEALLOCATE (x_scf)
1015 524 : 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 524 : kp, karray, mpi_group)
1023 :
1024 : !De-allocations
1025 524 : DEALLOCATE (kp)
1026 :
1027 524 : 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 78784 : 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 483188944 : zw = 0.0_dp
1049 78784 : i3 = 0
1050 4852144 : DO l3 = 1, n3, 2
1051 4773360 : i3 = i3 + 1
1052 4773360 : i03r = ABS(l3 - n3/2 - 1) + 1
1053 4773360 : i03i = ABS(l3 - n3/2) + 1
1054 142616224 : DO l1 = 1, nfft
1055 137764080 : i01 = ABS(l1 - 1 + i1 - n1/2 - 1) + 1
1056 137764080 : zw(1, l1, i3) = zf(i01, i03r)
1057 142537440 : zw(2, l1, i3) = zf(i01, i03i)
1058 : END DO
1059 : END DO
1060 :
1061 78784 : 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 524 : 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 524 : INTEGER, ALLOCATABLE, DIMENSION(:) :: after1, after2, after3, before1, &
1112 524 : before2, before3, now1, now2, now3
1113 : REAL(kind=dp) :: twopion
1114 524 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cosinarr, trig1, trig2, trig3
1115 524 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zt, zw
1116 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: zmpi2
1117 : REAL(KIND=dp), ALLOCATABLE, &
1118 524 : 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 524 : CPASSERT(nd1 >= n1)
1129 524 : CPASSERT(nd2 >= n2)
1130 524 : CPASSERT(nd3 >= n3/2 + 1)
1131 524 : CPASSERT(MOD(nd3, nproc) == 0)
1132 524 : CPASSERT(MOD(nd2, nproc) == 0)
1133 : MARK_USED(nd1)
1134 :
1135 : !defining work arrays dimensions
1136 524 : ncache = ncache_optimal
1137 524 : IF (ncache <= MAX(n1, n2, n3/2)*4) ncache = MAX(n1, n2, n3/2)*4
1138 524 : lzt = n2
1139 524 : IF (MOD(n2, 2) == 0) lzt = lzt + 1
1140 524 : IF (MOD(n2, 4) == 0) lzt = lzt + 1
1141 :
1142 : !Allocations
1143 524 : ALLOCATE (trig1(2, ctrig_length))
1144 524 : ALLOCATE (after1(7))
1145 524 : ALLOCATE (now1(7))
1146 524 : ALLOCATE (before1(7))
1147 524 : ALLOCATE (trig2(2, ctrig_length))
1148 524 : ALLOCATE (after2(7))
1149 524 : ALLOCATE (now2(7))
1150 524 : ALLOCATE (before2(7))
1151 524 : ALLOCATE (trig3(2, ctrig_length))
1152 524 : ALLOCATE (after3(7))
1153 524 : ALLOCATE (now3(7))
1154 524 : ALLOCATE (before3(7))
1155 2096 : ALLOCATE (zw(2, ncache/4, 2))
1156 2096 : ALLOCATE (zt(2, lzt, n1))
1157 2620 : ALLOCATE (zmpi2(2, n1, nd2/nproc, nd3))
1158 2096 : ALLOCATE (cosinarr(2, n3/2))
1159 524 : IF (nproc > 1) THEN
1160 2160 : ALLOCATE (zmpi1(2, n1, nd2/nproc, nd3/nproc, nproc))
1161 360 : zmpi1 = 0.0_dp
1162 : END IF
1163 :
1164 524 : zmpi2 = 0.0_dp
1165 : !calculating the FFT work arrays (beware on the HalFFT in n3 dimension)
1166 524 : CALL ctrig(n3/2, trig3, after3, before3, now3, 1, ic3)
1167 524 : CALL ctrig(n1, trig1, after1, before1, now1, 1, ic1)
1168 524 : CALL ctrig(n2, trig2, after2, before2, now2, 1, ic2)
1169 :
1170 : !Calculating array of phases for HalFFT decoding
1171 524 : twopion = 8._dp*ATAN(1._dp)/REAL(n3, KIND=dp)
1172 23208 : DO i3 = 1, n3/2
1173 22684 : cosinarr(1, i3) = COS(twopion*(i3 - 1))
1174 23208 : cosinarr(2, i3) = -SIN(twopion*(i3 - 1))
1175 : END DO
1176 :
1177 : !transform along z axis
1178 :
1179 524 : lot = ncache/(2*n3)
1180 524 : CPASSERT(lot >= 1)
1181 :
1182 28226 : DO j2 = 1, nd2/nproc
1183 : !this condition ensures that we manage only the interesting part for the FFT
1184 28226 : IF (iproc*(nd2/nproc) + j2 <= n2) THEN
1185 106486 : DO i1 = 1, n1, lot
1186 78784 : ma = i1
1187 78784 : mb = MIN(i1 + (lot - 1), n1)
1188 78784 : nfft = mb - ma + 1
1189 :
1190 : !inserting real data into complex array of half length
1191 : !input: I1,I3,J2,(Jp2)
1192 :
1193 78784 : CALL inserthalf(n1, n3, lot, nfft, i1, zf(1, 1, j2), zw(1, 1, 1))
1194 :
1195 : !performing FFT
1196 78784 : inzee = 1
1197 287648 : DO i = 1, ic3
1198 : CALL fftstp(lot, nfft, n3/2, lot, n3/2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1199 208864 : trig3, after3(i), now3(i), before3(i), 1)
1200 287648 : 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 106486 : 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 524 : IF (nproc > 1) THEN
1216 : !communication scheduling
1217 : CALL mpi_group%alltoall(zmpi2, &!2*n1*(nd2/nproc)*(nd3/nproc), &
1218 360 : zmpi1, 2*n1*(nd2/nproc)*(nd3/nproc))
1219 : ! output: I1,J2,j3,Jp2,(jp3)
1220 : END IF
1221 :
1222 15096 : DO j3 = 1, nd3/nproc
1223 : !this condition ensures that we manage only the interesting part for the FFT
1224 15096 : IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
1225 14392 : Jp2st = 1
1226 14392 : J2st = 1
1227 :
1228 : !transform along x axis
1229 14392 : lot = ncache/(4*n1)
1230 14392 : CPASSERT(lot >= 1)
1231 :
1232 92197 : DO j = 1, n2, lot
1233 77805 : ma = j
1234 77805 : mb = MIN(j + (lot - 1), n2)
1235 77805 : nfft = mb - ma + 1
1236 :
1237 : !reverse ordering
1238 : !input: I1,J2,j3,Jp2,(jp3)
1239 77805 : IF (nproc == 1) THEN
1240 19380 : CALL mpiswitch(j3, nfft, Jp2st, J2st, lot, n1, nd2, nd3, nproc, zmpi2, zw(1, 1, 1))
1241 : ELSE
1242 58425 : 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 77805 : inzee = 1
1249 253447 : DO i = 1, ic1 - 1
1250 : CALL fftstp(lot, nfft, n1, lot, n1, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1251 175642 : trig1, after1(i), now1(i), before1(i), 1)
1252 253447 : inzee = 3 - inzee
1253 : END DO
1254 : !storing the last step into zt
1255 77805 : i = ic1
1256 : CALL fftstp(lot, nfft, n1, lzt, n1, zw(1, 1, inzee), zt(1, j, 1), &
1257 92197 : 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 14392 : lot = ncache/(4*n2)
1263 14392 : CPASSERT(lot >= 1)
1264 :
1265 58965 : DO j = 1, nk1, lot
1266 44573 : ma = j
1267 44573 : mb = MIN(j + (lot - 1), nk1)
1268 44573 : nfft = mb - ma + 1
1269 :
1270 : !reverse ordering
1271 : !input: I2,i1,j3,(jp3)
1272 44573 : 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 44573 : inzee = 1
1278 190890 : DO i = 1, ic2
1279 : CALL fftstp(lot, nfft, n2, lot, n2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1280 146317 : trig2, after2(i), now2(i), before2(i), 1)
1281 190890 : inzee = 3 - inzee
1282 : END DO
1283 :
1284 58965 : 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 524 : DEALLOCATE (trig1)
1293 524 : DEALLOCATE (after1)
1294 524 : DEALLOCATE (now1)
1295 524 : DEALLOCATE (before1)
1296 524 : DEALLOCATE (trig2)
1297 524 : DEALLOCATE (after2)
1298 524 : DEALLOCATE (now2)
1299 524 : DEALLOCATE (before2)
1300 524 : DEALLOCATE (trig3)
1301 524 : DEALLOCATE (after3)
1302 524 : DEALLOCATE (now3)
1303 524 : DEALLOCATE (before3)
1304 524 : DEALLOCATE (zmpi2)
1305 524 : DEALLOCATE (zw)
1306 524 : DEALLOCATE (zt)
1307 524 : DEALLOCATE (cosinarr)
1308 524 : IF (nproc > 1) DEALLOCATE (zmpi1)
1309 :
1310 524 : 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 44573 : 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 2737452 : DO i = 1, nk2
1330 39158011 : DO j = 1, nfft
1331 39113438 : zout(j, i) = zin(1, j, i)
1332 : END DO
1333 : END DO
1334 :
1335 44573 : 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 44573 : 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 721084 : DO 200, j = 1, nfft
1354 72164607 : DO 100, i = 1, n2
1355 71488096 : zw(1, j, i) = zt(1, i, j)
1356 71488096 : zw(2, j, i) = zt(2, i, j)
1357 676511 : 100 CONTINUE
1358 44573 : 200 CONTINUE
1359 44573 : 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 77805 : 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 77805 : mfft = 0
1384 101013 : DO 300, Jp2 = Jp2st, nproc
1385 1410859 : DO 200, J2 = J2st, nd2/nproc
1386 1387651 : mfft = mfft + 1
1387 1387651 : IF (mfft > nfft) THEN
1388 63413 : Jp2st = Jp2
1389 63413 : J2st = J2
1390 63413 : RETURN
1391 : END IF
1392 141651954 : DO 100, I1 = 1, n1
1393 140327716 : zw(1, mfft, I1) = zmpi1(1, I1, J2, j3, Jp2)
1394 140327716 : zw(2, mfft, I1) = zmpi1(2, I1, J2, j3, Jp2)
1395 1324238 : 100 CONTINUE
1396 23208 : 200 CONTINUE
1397 23208 : J2st = 1
1398 14392 : 300 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 524 : 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 524 : p(1) = 4.96142640560223544e19_dp
1420 524 : p(2) = 1.37454269147978052e19_dp
1421 524 : p(3) = 7.58610013441204679e18_dp
1422 524 : p(4) = 4.42040691347806996e18_dp
1423 524 : p(5) = 2.61986077948367892e18_dp
1424 524 : p(6) = 1.56320138155496681e18_dp
1425 524 : p(7) = 9.35645215863028402e17_dp
1426 524 : p(8) = 5.60962910452691703e17_dp
1427 524 : p(9) = 3.3666225119686761e17_dp
1428 524 : p(10) = 2.0218253197947866e17_dp
1429 524 : p(11) = 1.21477756091902017e17_dp
1430 524 : p(12) = 7.3012982513608503e16_dp
1431 524 : p(13) = 4.38951893556421099e16_dp
1432 524 : p(14) = 2.63949482512262325e16_dp
1433 524 : p(15) = 1.58742054072786174e16_dp
1434 524 : p(16) = 9.54806587737665531e15_dp
1435 524 : p(17) = 5.74353712364571709e15_dp
1436 524 : p(18) = 3.455214877389445e15_dp
1437 524 : p(19) = 2.07871658520326804e15_dp
1438 524 : p(20) = 1.25064667315629928e15_dp
1439 524 : p(21) = 7.52469429541933745e14_dp
1440 524 : p(22) = 4.5274603337253175e14_dp
1441 524 : p(23) = 2.72414006900059548e14_dp
1442 524 : p(24) = 1.63912168349216752e14_dp
1443 524 : p(25) = 9.86275802590865738e13_dp
1444 524 : p(26) = 5.93457701624974985e13_dp
1445 524 : p(27) = 3.5709554322296296e13_dp
1446 524 : p(28) = 2.14872890367310454e13_dp
1447 524 : p(29) = 1.29294719957726902e13_dp
1448 524 : p(30) = 7.78003375426361016e12_dp
1449 524 : p(31) = 4.68148199759876704e12_dp
1450 524 : p(32) = 2.8169955024829868e12_dp
1451 524 : p(33) = 1.69507790481958464e12_dp
1452 524 : p(34) = 1.01998486064607581e12_dp
1453 524 : p(35) = 6.13759486539856459e11_dp
1454 524 : p(36) = 3.69320183828682544e11_dp
1455 524 : p(37) = 2.22232783898905102e11_dp
1456 524 : p(38) = 1.33725247623668682e11_dp
1457 524 : p(39) = 8.0467192739036288e10_dp
1458 524 : p(40) = 4.84199582415144143e10_dp
1459 524 : p(41) = 2.91360091170559564e10_dp
1460 524 : p(42) = 1.75321747475309216e10_dp
1461 524 : p(43) = 1.0549735552210995e10_dp
1462 524 : p(44) = 6.34815321079006586e9_dp
1463 524 : p(45) = 3.81991113733594231e9_dp
1464 524 : p(46) = 2.29857747533101109e9_dp
1465 524 : p(47) = 1.38313653595483694e9_dp
1466 524 : p(48) = 8.32282908580025358e8_dp
1467 524 : p(49) = 5.00814519374587467e8_dp
1468 524 : p(50) = 3.01358090773319025e8_dp
1469 524 : p(51) = 1.81337994217503535e8_dp
1470 524 : p(52) = 1.09117589961086823e8_dp
1471 524 : p(53) = 6.56599771718640323e7_dp
1472 524 : p(54) = 3.95099693638497164e7_dp
1473 524 : p(55) = 2.37745694710665991e7_dp
1474 524 : p(56) = 1.43060135285912813e7_dp
1475 524 : p(57) = 8.60844290313506695e6_dp
1476 524 : p(58) = 5.18000974075383424e6_dp
1477 524 : p(59) = 3.116998193057466e6_dp
1478 524 : p(60) = 1.87560993870024029e6_dp
1479 524 : p(61) = 1.12862197183979562e6_dp
1480 524 : p(62) = 679132.441326077231_dp
1481 524 : p(63) = 408658.421279877969_dp
1482 524 : p(64) = 245904.473450669789_dp
1483 524 : p(65) = 147969.568088321005_dp
1484 524 : p(66) = 89038.612357311147_dp
1485 524 : p(67) = 53577.7362552358895_dp
1486 524 : p(68) = 32239.6513926914668_dp
1487 524 : p(69) = 19399.7580852362791_dp
1488 524 : p(70) = 11673.5323603058634_dp
1489 524 : p(71) = 7024.38438577707758_dp
1490 524 : p(72) = 4226.82479307685999_dp
1491 524 : p(73) = 2543.43254175354295_dp
1492 524 : p(74) = 1530.47486269122675_dp
1493 524 : p(75) = 920.941785160749482_dp
1494 524 : p(76) = 554.163803906291646_dp
1495 524 : p(77) = 333.46029740785694_dp
1496 524 : p(78) = 200.6550575335041_dp
1497 524 : p(79) = 120.741366914147284_dp
1498 524 : p(80) = 72.6544243200329916_dp
1499 524 : p(81) = 43.7187810415471025_dp
1500 524 : p(82) = 26.3071631447061043_dp
1501 524 : p(83) = 15.8299486353816329_dp
1502 524 : p(84) = 9.52493152341244004_dp
1503 524 : p(85) = 5.72200417067776041_dp
1504 524 : p(86) = 3.36242234070940928_dp
1505 524 : p(87) = 1.75371394604499472_dp
1506 524 : p(88) = 0.64705932650658966_dp
1507 524 : p(89) = 0.072765905943708247_dp
1508 : !
1509 524 : w(1) = 47.67445484528304247e10_dp
1510 524 : w(2) = 11.37485774750442175e9_dp
1511 524 : w(3) = 78.64340976880190239e8_dp
1512 524 : w(4) = 46.27335788759590498e8_dp
1513 524 : w(5) = 24.7380464827152951e8_dp
1514 524 : w(6) = 13.62904116438987719e8_dp
1515 524 : w(7) = 92.79560029045882433e8_dp
1516 524 : w(8) = 52.15931216254660251e8_dp
1517 524 : w(9) = 31.67018011061666244e8_dp
1518 524 : w(10) = 1.29291036801493046e8_dp
1519 524 : w(11) = 1.00139319988015862e8_dp
1520 524 : w(12) = 7.75892350510188341e7_dp
1521 524 : w(13) = 6.01333567950731271e7_dp
1522 524 : w(14) = 4.66141178654796875e7_dp
1523 524 : w(15) = 3.61398903394911448e7_dp
1524 524 : w(16) = 2.80225846672956389e7_dp
1525 524 : w(17) = 2.1730509180930247e7_dp
1526 524 : w(18) = 1.68524482625876965e7_dp
1527 524 : w(19) = 1.30701489345870338e7_dp
1528 524 : w(20) = 1.01371784832269282e7_dp
1529 524 : w(21) = 7.86264116300379329e6_dp
1530 524 : w(22) = 6.09861667912273717e6_dp
1531 524 : w(23) = 4.73045784039455683e6_dp
1532 524 : w(24) = 3.66928949951594161e6_dp
1533 524 : w(25) = 2.8462050836230259e6_dp
1534 524 : w(26) = 2.20777394798527011e6_dp
1535 524 : w(27) = 1.71256191589205524e6_dp
1536 524 : w(28) = 1.32843556197737076e6_dp
1537 524 : w(29) = 1.0304731275955989e6_dp
1538 524 : w(30) = 799345.206572271448_dp
1539 524 : w(31) = 620059.354143595343_dp
1540 524 : w(32) = 480986.704107449333_dp
1541 524 : w(33) = 373107.167700228515_dp
1542 524 : w(34) = 289424.08337412132_dp
1543 524 : w(35) = 224510.248231581788_dp
1544 524 : w(36) = 174155.825690028966_dp
1545 524 : w(37) = 135095.256919654065_dp
1546 524 : w(38) = 104795.442776800312_dp
1547 524 : w(39) = 81291.4458222430418_dp
1548 524 : w(40) = 63059.0493649328682_dp
1549 524 : w(41) = 48915.9040455329689_dp
1550 524 : w(42) = 37944.8484018048756_dp
1551 524 : w(43) = 29434.4290473253969_dp
1552 524 : w(44) = 22832.7622054490044_dp
1553 524 : w(45) = 17711.743950151233_dp
1554 524 : w(46) = 13739.287867104177_dp
1555 524 : w(47) = 10657.7895710752585_dp
1556 524 : w(48) = 8267.42141053961834_dp
1557 524 : w(49) = 6413.17397520136448_dp
1558 524 : w(50) = 4974.80402838654277_dp
1559 524 : w(51) = 3859.03698188553047_dp
1560 524 : w(52) = 2993.51824493299154_dp
1561 524 : w(53) = 2322.1211966811754_dp
1562 524 : w(54) = 1801.30750964719641_dp
1563 524 : w(55) = 1397.30379659817038_dp
1564 524 : w(56) = 1083.91149143250697_dp
1565 524 : w(57) = 840.807939169209188_dp
1566 524 : w(58) = 652.228524366749422_dp
1567 524 : w(59) = 505.944376983506128_dp
1568 524 : w(60) = 392.469362317941064_dp
1569 524 : w(61) = 304.444930257324312_dp
1570 524 : w(62) = 236.162932842453601_dp
1571 524 : w(63) = 183.195466078603525_dp
1572 524 : w(64) = 142.107732186551471_dp
1573 524 : w(65) = 110.23530215723992_dp
1574 524 : w(66) = 85.5113346705382257_dp
1575 524 : w(67) = 66.3325469806696621_dp
1576 524 : w(68) = 51.4552463353841373_dp
1577 524 : w(69) = 39.9146798429449273_dp
1578 524 : w(70) = 30.9624728409162095_dp
1579 524 : w(71) = 24.018098812215013_dp
1580 524 : w(72) = 18.6312338024296588_dp
1581 524 : w(73) = 14.4525541233150501_dp
1582 524 : w(74) = 11.2110836519105938_dp
1583 524 : w(75) = 8.69662175848497178_dp
1584 524 : w(76) = 6.74611236165731961_dp
1585 524 : w(77) = 5.23307018057529994_dp
1586 524 : w(78) = 4.05937850501539556_dp
1587 524 : w(79) = 3.14892659076635714_dp
1588 524 : w(80) = 2.44267408211071604_dp
1589 524 : w(81) = 1.89482240522855261_dp
1590 524 : w(82) = 1.46984505907050079_dp
1591 524 : w(83) = 1.14019261330527007_dp
1592 524 : w(84) = 0.884791217422925293_dp
1593 524 : w(85) = 0.692686387080616483_dp
1594 524 : w(86) = 0.585244576897023282_dp
1595 524 : w(87) = 0.576182522545327589_dp
1596 524 : w(88) = 0.596688817388997178_dp
1597 524 : w(89) = 0.607879901151108771_dp
1598 : !
1599 : !
1600 524 : urange = 1._dp
1601 524 : drange = 1e-08_dp
1602 524 : acc = 1e-08_dp
1603 : !
1604 524 : RETURN
1605 : END SUBROUTINE gequad
1606 :
1607 : END MODULE ps_wavelet_kernel
|
