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