(git:71c3ab0)
Loading...
Searching...
No Matches
ps_wavelet_base.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Creates the wavelet kernel for the wavelet based poisson solver.
10!> \author Florian Schiffmann (09.2007,fschiff)
11! **************************************************************************************************
12MODULE ps_wavelet_base
13
14 USE kinds, ONLY: dp
15 USE mathconstants, ONLY: pi
16 USE message_passing, ONLY: mp_comm_type
17 USE ps_wavelet_fft3d, ONLY: ctrig,&
18 ctrig_length,&
19 fftstp
20#include "../base/base_uses.f90"
21
22 IMPLICIT NONE
23
24 PRIVATE
25
26 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ps_wavelet_base'
27
28 PUBLIC :: scramble_unpack, p_poissonsolver, s_poissonsolver, f_poissonsolver
29
30CONTAINS
31
32! **************************************************************************************************
33!> \brief ...
34!> \param n1 ...
35!> \param n2 ...
36!> \param n3 ...
37!> \param nd1 ...
38!> \param nd2 ...
39!> \param nd3 ...
40!> \param md1 ...
41!> \param md2 ...
42!> \param md3 ...
43!> \param nproc ...
44!> \param iproc ...
45!> \param zf ...
46!> \param scal ...
47!> \param hx ...
48!> \param hy ...
49!> \param hz ...
50!> \param mpi_group ...
51! **************************************************************************************************
52 SUBROUTINE p_poissonsolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, zf &
53 , scal, hx, hy, hz, mpi_group)
54 INTEGER, INTENT(in) :: n1, n2, n3, nd1, nd2, nd3, md1, md2, &
55 md3, nproc, iproc
56 REAL(kind=dp), DIMENSION(md1, md3, md2/nproc), &
57 INTENT(inout) :: zf
58 REAL(kind=dp), INTENT(in) :: scal, hx, hy, hz
59
60 CLASS(mp_comm_type), INTENT(in) :: mpi_group
61
62 INTEGER, PARAMETER :: ncache_optimal = 8*1024
63
64 INTEGER :: i, i1, i3, ic1, ic2, ic3, inzee, j, j2, &
65 j2stb, j2stf, j3, jp2stb, jp2stf, lot, &
66 lzt, ma, mb, ncache, nfft
67 INTEGER, ALLOCATABLE, DIMENSION(:) :: after1, after2, after3, before1, &
68 before2, before3, now1, now2, now3
69 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: btrig1, btrig2, btrig3, ftrig1, ftrig2, &
70 ftrig3
71 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zt, zw
72 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: zmpi2
73 REAL(kind=dp), ALLOCATABLE, &
74 DIMENSION(:, :, :, :, :) :: zmpi1
75
76 IF (nd1 < n1/2 + 1) cpabort("Parallel convolution:ERROR:nd1")
77 IF (nd2 < n2/2 + 1) cpabort("Parallel convolution:ERROR:nd2")
78 IF (nd3 < n3/2 + 1) cpabort("Parallel convolution:ERROR:nd3")
79 IF (md1 < n1) cpabort("Parallel convolution:ERROR:md1")
80 IF (md2 < n2) cpabort("Parallel convolution:ERROR:md2")
81 IF (md3 < n3) cpabort("Parallel convolution:ERROR:md3")
82 IF (mod(nd3, nproc) /= 0) cpabort("Parallel convolution:ERROR:nd3")
83 IF (mod(md2, nproc) /= 0) cpabort("Parallel convolution:ERROR:md2")
84
85 !defining work arrays dimensions
86 ncache = ncache_optimal
87 IF (ncache <= max(n1, n2, n3)*4) ncache = max(n1, n2, n3)*4
88
89 lzt = n2
90 IF (mod(n2, 2) == 0) lzt = lzt + 1
91 IF (mod(n2, 4) == 0) lzt = lzt + 1 !maybe this is useless
92
93 !Allocations
94 ALLOCATE (btrig1(2, ctrig_length))
95 ALLOCATE (ftrig1(2, ctrig_length))
96 ALLOCATE (after1(7))
97 ALLOCATE (now1(7))
98 ALLOCATE (before1(7))
99 ALLOCATE (btrig2(2, ctrig_length))
100 ALLOCATE (ftrig2(2, ctrig_length))
101 ALLOCATE (after2(7))
102 ALLOCATE (now2(7))
103 ALLOCATE (before2(7))
104 ALLOCATE (btrig3(2, ctrig_length))
105 ALLOCATE (ftrig3(2, ctrig_length))
106 ALLOCATE (after3(7))
107 ALLOCATE (now3(7))
108 ALLOCATE (before3(7))
109 ALLOCATE (zw(2, ncache/4, 2))
110 ALLOCATE (zt(2, lzt, n1))
111 ALLOCATE (zmpi2(2, n1, md2/nproc, nd3))
112 IF (nproc > 1) ALLOCATE (zmpi1(2, n1, md2/nproc, nd3/nproc, nproc))
113
114 !calculating the FFT work arrays (beware on the HalFFT in n3 dimension)
115 CALL ctrig(n3, btrig3, after3, before3, now3, 1, ic3)
116 CALL ctrig(n1, btrig1, after1, before1, now1, 1, ic1)
117 CALL ctrig(n2, btrig2, after2, before2, now2, 1, ic2)
118 DO j = 1, n1
119 ftrig1(1, j) = btrig1(1, j)
120 ftrig1(2, j) = -btrig1(2, j)
121 END DO
122 DO j = 1, n2
123 ftrig2(1, j) = btrig2(1, j)
124 ftrig2(2, j) = -btrig2(2, j)
125 END DO
126 DO j = 1, n3
127 ftrig3(1, j) = btrig3(1, j)
128 ftrig3(2, j) = -btrig3(2, j)
129 END DO
130
131 ! transform along z axis
132 lot = ncache/(4*n3)
133 IF (lot < 1) THEN
134 CALL cp_abort(__location__, &
135 'convolxc_off:ncache has to be enlarged to be able to hold at '// &
136 'least one 1-d FFT of this size even though this will '// &
137 'reduce the performance for shorter transform lengths')
138 END IF
139 DO j2 = 1, md2/nproc
140 !this condition ensures that we manage only the interesting part for the FFT
141 IF (iproc*(md2/nproc) + j2 <= n2) THEN
142 DO i1 = 1, n1, lot
143 ma = i1
144 mb = min(i1 + (lot - 1), n1)
145 nfft = mb - ma + 1
146 !inserting real data into complex array of half length
147 CALL p_fill_upcorn(md1, md3, lot, nfft, n3, zf(i1, 1, j2), zw(1, 1, 1))
148
149 !performing FFT
150 !input: I1,I3,J2,(Jp2)
151 inzee = 1
152 DO i = 1, ic3
153 CALL fftstp(lot, nfft, n3, lot, n3, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
154 btrig3, after3(i), now3(i), before3(i), 1)
155 inzee = 3 - inzee
156 END DO
157
158 !output: I1,i3,J2,(Jp2)
159 !exchanging components
160 !input: I1,i3,J2,(Jp2)
161 CALL scramble_p(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zw(1, 1, inzee), zmpi2)
162 !output: I1,J2,i3,(Jp2)
163 END DO
164 END IF
165 END DO
166
167 !Interprocessor data transposition
168 !input: I1,J2,j3,jp3,(Jp2)
169 IF (nproc > 1) THEN
170 !communication scheduling
171 CALL mpi_group%alltoall(zmpi2, zmpi1, 2*n1*(md2/nproc)*(nd3/nproc))
172 END IF
173 !output: I1,J2,j3,Jp2,(jp3)
174
175 !now each process perform complete convolution of its planes
176 DO j3 = 1, nd3/nproc
177 !this condition ensures that we manage only the interesting part for the FFT
178 IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
179 jp2stb = 1
180 j2stb = 1
181 jp2stf = 1
182 j2stf = 1
183
184 ! transform along x axis
185 lot = ncache/(4*n1)
186 IF (lot < 1) THEN
187 CALL cp_abort(__location__, &
188 'convolxc_off:ncache has to be enlarged to be able to hold at '// &
189 'least one 1-d FFT of this size even though this will '// &
190 'reduce the performance for shorter transform lengths')
191 END IF
192
193 DO j = 1, n2, lot
194 ma = j
195 mb = min(j + (lot - 1), n2)
196 nfft = mb - ma + 1
197
198 !reverse index ordering, leaving the planes to be transformed at the end
199 !input: I1,J2,j3,Jp2,(jp3)
200 IF (nproc == 1) THEN
201 CALL p_mpiswitch_upcorn(j3, nfft, jp2stb, j2stb, lot, n1, md2, nd3, nproc, zmpi2, zw(1, 1, 1))
202 ELSE
203 CALL p_mpiswitch_upcorn(j3, nfft, jp2stb, j2stb, lot, n1, md2, nd3, nproc, zmpi1, zw(1, 1, 1))
204 END IF
205 !output: J2,Jp2,I1,j3,(jp3)
206
207 !performing FFT
208 !input: I2,I1,j3,(jp3)
209 inzee = 1
210 DO i = 1, ic1 - 1
211 CALL fftstp(lot, nfft, n1, lot, n1, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
212 btrig1, after1(i), now1(i), before1(i), 1)
213 inzee = 3 - inzee
214 END DO
215
216 !storing the last step into zt array
217 i = ic1
218 CALL fftstp(lot, nfft, n1, lzt, n1, zw(1, 1, inzee), zt(1, j, 1), &
219 btrig1, after1(i), now1(i), before1(i), 1)
220 !output: I2,i1,j3,(jp3)
221 END DO
222
223 !transform along y axis
224 lot = ncache/(4*n2)
225 IF (lot < 1) THEN
226 CALL cp_abort(__location__, &
227 'convolxc_off:ncache has to be enlarged to be able to hold at '// &
228 'least one 1-d FFT of this size even though this will '// &
229 'reduce the performance for shorter transform lengths')
230 END IF
231
232 DO j = 1, n1, lot
233 ma = j
234 mb = min(j + (lot - 1), n1)
235 nfft = mb - ma + 1
236
237 !reverse ordering
238 !input: I2,i1,j3,(jp3)
239 CALL p_switch_upcorn(nfft, n2, lot, n1, lzt, zt(1, 1, j), zw(1, 1, 1))
240 !output: i1,I2,j3,(jp3)
241
242 !performing FFT
243 !input: i1,I2,j3,(jp3)
244 inzee = 1
245 DO i = 1, ic2
246 CALL fftstp(lot, nfft, n2, lot, n2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
247 btrig2, after2(i), now2(i), before2(i), 1)
248 inzee = 3 - inzee
249 END DO
250 !output: i1,i2,j3,(jp3)
251
252 !Multiply with kernel in fourier space
253 i3 = iproc*(nd3/nproc) + j3
254 CALL p_multkernel(n1, n2, n3, lot, nfft, j, i3, zw(1, 1, inzee), hx, hy, hz)
255
256 !TRANSFORM BACK IN REAL SPACE
257
258 !transform along y axis
259 !input: i1,i2,j3,(jp3)
260 DO i = 1, ic2
261 CALL fftstp(lot, nfft, n2, lot, n2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
262 ftrig2, after2(i), now2(i), before2(i), -1)
263 inzee = 3 - inzee
264 END DO
265
266 !reverse ordering
267 !input: i1,I2,j3,(jp3)
268 CALL p_unswitch_downcorn(nfft, n2, lot, n1, lzt, zw(1, 1, inzee), zt(1, 1, j))
269 !output: I2,i1,j3,(jp3)
270 END DO
271
272 !transform along x axis
273 !input: I2,i1,j3,(jp3)
274 lot = ncache/(4*n1)
275 DO j = 1, n2, lot
276 ma = j
277 mb = min(j + (lot - 1), n2)
278 nfft = mb - ma + 1
279
280 !performing FFT
281 i = 1
282 CALL fftstp(lzt, nfft, n1, lot, n1, zt(1, j, 1), zw(1, 1, 1), &
283 ftrig1, after1(i), now1(i), before1(i), -1)
284
285 inzee = 1
286 DO i = 2, ic1
287 CALL fftstp(lot, nfft, n1, lot, n1, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
288 ftrig1, after1(i), now1(i), before1(i), -1)
289 inzee = 3 - inzee
290 END DO
291 !output: I2,I1,j3,(jp3)
292
293 !reverse ordering
294 !input: J2,Jp2,I1,j3,(jp3)
295 IF (nproc == 1) THEN
296 CALL p_unmpiswitch_downcorn(j3, nfft, jp2stf, j2stf, lot, n1, md2, nd3, nproc, zw(1, 1, inzee), zmpi2)
297 ELSE
298 CALL p_unmpiswitch_downcorn(j3, nfft, jp2stf, j2stf, lot, n1, md2, nd3, nproc, zw(1, 1, inzee), zmpi1)
299 END IF
300 ! output: I1,J2,j3,Jp2,(jp3)
301 END DO
302 END IF
303 END DO
304
305 !Interprocessor data transposition
306 !input: I1,J2,j3,Jp2,(jp3)
307 IF (nproc > 1) THEN
308 !communication scheduling
309 CALL mpi_group%alltoall(zmpi1, zmpi2, 2*n1*(md2/nproc)*(nd3/nproc))
310 END IF
311 !output: I1,J2,j3,jp3,(Jp2)
312 !transform along z axis
313 !input: I1,J2,i3,(Jp2)
314 lot = ncache/(4*n3)
315 DO j2 = 1, md2/nproc
316 !this condition ensures that we manage only the interesting part for the FFT
317 IF (iproc*(md2/nproc) + j2 <= n2) THEN
318 DO i1 = 1, n1, lot
319 ma = i1
320 mb = min(i1 + (lot - 1), n1)
321 nfft = mb - ma + 1
322
323 !reverse ordering
324 !input: I1,J2,i3,(Jp2)
325 CALL unscramble_p(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zmpi2, zw(1, 1, 1))
326 !output: I1,i3,J2,(Jp2)
327
328 !performing FFT
329 !input: I1,i3,J2,(Jp2)
330 inzee = 1
331 DO i = 1, ic3
332 CALL fftstp(lot, nfft, n3, lot, n3, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
333 ftrig3, after3(i), now3(i), before3(i), -1)
334 inzee = 3 - inzee
335 END DO
336 !output: I1,I3,J2,(Jp2)
337
338 !rebuild the output array
339 CALL p_unfill_downcorn(md1, md3, lot, nfft, n3, zw(1, 1, inzee), zf(i1, 1, j2), scal)
340
341 END DO
342 END IF
343 END DO
344
345 !De-allocations
346 DEALLOCATE (btrig1)
347 DEALLOCATE (ftrig1)
348 DEALLOCATE (after1)
349 DEALLOCATE (now1)
350 DEALLOCATE (before1)
351 DEALLOCATE (btrig2)
352 DEALLOCATE (ftrig2)
353 DEALLOCATE (after2)
354 DEALLOCATE (now2)
355 DEALLOCATE (before2)
356 DEALLOCATE (btrig3)
357 DEALLOCATE (ftrig3)
358 DEALLOCATE (after3)
359 DEALLOCATE (now3)
360 DEALLOCATE (before3)
361 DEALLOCATE (zmpi2)
362 DEALLOCATE (zw)
363 DEALLOCATE (zt)
364 IF (nproc > 1) DEALLOCATE (zmpi1)
365 ! call timing(iproc,'PSolv_comput ','OF')
366 END SUBROUTINE p_poissonsolver
367
368! **************************************************************************************************
369!> \brief ...
370!> \param j3 ...
371!> \param nfft ...
372!> \param Jp2stb ...
373!> \param J2stb ...
374!> \param lot ...
375!> \param n1 ...
376!> \param md2 ...
377!> \param nd3 ...
378!> \param nproc ...
379!> \param zmpi1 ...
380!> \param zw ...
381! **************************************************************************************************
382 SUBROUTINE p_mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot, n1, md2, nd3, nproc, zmpi1, zw)
383 INTEGER, INTENT(in) :: j3, nfft
384 INTEGER, INTENT(inout) :: jp2stb, j2stb
385 INTEGER, INTENT(in) :: lot, n1, md2, nd3, nproc
386 REAL(kind=dp), &
387 DIMENSION(2, n1, md2/nproc, nd3/nproc, nproc), &
388 INTENT(in) :: zmpi1
389 REAL(kind=dp), DIMENSION(2, lot, n1), &
390 INTENT(inout) :: zw
391
392 INTEGER :: i1, j2, jp2, mfft
393
394 mfft = 0
395 DO jp2 = jp2stb, nproc
396 DO j2 = j2stb, md2/nproc
397 mfft = mfft + 1
398 IF (mfft > nfft) THEN
399 jp2stb = jp2
400 j2stb = j2
401 RETURN
402 END IF
403 DO i1 = 1, n1
404 zw(1, mfft, i1) = zmpi1(1, i1, j2, j3, jp2)
405 zw(2, mfft, i1) = zmpi1(2, i1, j2, j3, jp2)
406 END DO
407 END DO
408 j2stb = 1
409 END DO
410 END SUBROUTINE p_mpiswitch_upcorn
411
412! **************************************************************************************************
413!> \brief ...
414!> \param nfft ...
415!> \param n2 ...
416!> \param lot ...
417!> \param n1 ...
418!> \param lzt ...
419!> \param zt ...
420!> \param zw ...
421! **************************************************************************************************
422 SUBROUTINE p_switch_upcorn(nfft, n2, lot, n1, lzt, zt, zw)
423 INTEGER, INTENT(in) :: nfft, n2, lot, n1, lzt
424 REAL(kind=dp), DIMENSION(2, lzt, n1), INTENT(in) :: zt
425 REAL(kind=dp), DIMENSION(2, lot, n2), &
426 INTENT(inout) :: zw
427
428 INTEGER :: i, j
429
430 DO j = 1, nfft
431 DO i = 1, n2
432 zw(1, j, i) = zt(1, i, j)
433 zw(2, j, i) = zt(2, i, j)
434 END DO
435 END DO
436
437 END SUBROUTINE p_switch_upcorn
438
439! **************************************************************************************************
440!> \brief ...
441!> \param nfft ...
442!> \param n2 ...
443!> \param lot ...
444!> \param n1 ...
445!> \param lzt ...
446!> \param zw ...
447!> \param zt ...
448! **************************************************************************************************
449 SUBROUTINE p_unswitch_downcorn(nfft, n2, lot, n1, lzt, zw, zt)
450 INTEGER, INTENT(in) :: nfft, n2, lot, n1, lzt
451 REAL(kind=dp), DIMENSION(2, lot, n2), INTENT(in) :: zw
452 REAL(kind=dp), DIMENSION(2, lzt, n1), &
453 INTENT(inout) :: zt
454
455 INTEGER :: i, j
456
457 DO j = 1, nfft
458 DO i = 1, n2
459 zt(1, i, j) = zw(1, j, i)
460 zt(2, i, j) = zw(2, j, i)
461 END DO
462 END DO
463
464 END SUBROUTINE p_unswitch_downcorn
465
466! **************************************************************************************************
467!> \brief ...
468!> \param j3 ...
469!> \param nfft ...
470!> \param Jp2stf ...
471!> \param J2stf ...
472!> \param lot ...
473!> \param n1 ...
474!> \param md2 ...
475!> \param nd3 ...
476!> \param nproc ...
477!> \param zw ...
478!> \param zmpi1 ...
479! **************************************************************************************************
480 SUBROUTINE p_unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot, n1, md2, nd3, nproc, zw, zmpi1)
481 INTEGER, INTENT(in) :: j3, nfft
482 INTEGER, INTENT(inout) :: jp2stf, j2stf
483 INTEGER, INTENT(in) :: lot, n1, md2, nd3, nproc
484 REAL(kind=dp), DIMENSION(2, lot, n1), INTENT(in) :: zw
485 REAL(kind=dp), &
486 DIMENSION(2, n1, md2/nproc, nd3/nproc, nproc), &
487 INTENT(inout) :: zmpi1
488
489 INTEGER :: i1, j2, jp2, mfft
490
491 mfft = 0
492 DO jp2 = jp2stf, nproc
493 DO j2 = j2stf, md2/nproc
494 mfft = mfft + 1
495 IF (mfft > nfft) THEN
496 jp2stf = jp2
497 j2stf = j2
498 RETURN
499 END IF
500 DO i1 = 1, n1
501 zmpi1(1, i1, j2, j3, jp2) = zw(1, mfft, i1)
502 zmpi1(2, i1, j2, j3, jp2) = zw(2, mfft, i1)
503 END DO
504 END DO
505 j2stf = 1
506 END DO
507 END SUBROUTINE p_unmpiswitch_downcorn
508
509! **************************************************************************************************
510!> \brief (Based on suitable modifications of S.Goedecker routines)
511!> Restore data into output array
512!> \param md1 Dimensions of the undistributed part of the real grid
513!> \param md3 Dimensions of the undistributed part of the real grid
514!> \param lot ...
515!> \param nfft number of planes
516!> \param n3 (twice the) dimension of the last FFTtransform
517!> \param zw FFT work array
518!> \param zf Original distributed density as well as
519!> Distributed solution of the poisson equation (inout)
520!> \param scal Needed to achieve unitarity and correct dimensions
521!> \date February 2006
522!> \author S. Goedecker, L. Genovese
523!> \note Assuming that high frequencies are in the corners
524!> and that n3 is multiple of 4
525!>
526!> RESTRICTIONS on USAGE
527!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
528!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
529!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
530!> This file is distributed under the terms of the
531!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
532! **************************************************************************************************
533 SUBROUTINE p_unfill_downcorn(md1, md3, lot, nfft, n3, zw, zf, scal)
534 INTEGER, INTENT(in) :: md1, md3, lot, nfft, n3
535 REAL(kind=dp), DIMENSION(2, lot, n3), INTENT(in) :: zw
536 REAL(kind=dp), DIMENSION(md1, md3), INTENT(inout) :: zf
537 REAL(kind=dp), INTENT(in) :: scal
538
539 INTEGER :: i1, i3
540 REAL(kind=dp) :: pot1
541
542 DO i3 = 1, n3
543 DO i1 = 1, nfft
544 pot1 = scal*zw(1, i1, i3)
545 zf(i1, i3) = pot1
546 END DO
547 END DO
548
549 END SUBROUTINE p_unfill_downcorn
550
551! **************************************************************************************************
552!> \brief ...
553!> \param md1 ...
554!> \param md3 ...
555!> \param lot ...
556!> \param nfft ...
557!> \param n3 ...
558!> \param zf ...
559!> \param zw ...
560! **************************************************************************************************
561 SUBROUTINE p_fill_upcorn(md1, md3, lot, nfft, n3, zf, zw)
562 INTEGER, INTENT(in) :: md1, md3, lot, nfft, n3
563 REAL(kind=dp), DIMENSION(md1, md3), INTENT(in) :: zf
564 REAL(kind=dp), DIMENSION(2, lot, n3), &
565 INTENT(inout) :: zw
566
567 INTEGER :: i1, i3
568
569 DO i3 = 1, n3
570 DO i1 = 1, nfft
571 zw(1, i1, i3) = zf(i1, i3)
572 zw(2, i1, i3) = 0._dp
573 END DO
574 END DO
575
576 END SUBROUTINE p_fill_upcorn
577
578! **************************************************************************************************
579!> \brief (Based on suitable modifications of S.Goedecker routines)
580!> Assign the correct planes to the work array zmpi2
581!> in order to prepare for interprocessor data transposition.
582!> \param i1 Starting points of the plane and number of remaining lines
583!> \param j2 Starting points of the plane and number of remaining lines
584!> \param lot Starting points of the plane and number of remaining lines
585!> \param nfft Starting points of the plane and number of remaining lines
586!> \param n1 logical dimension of the FFT transform, reference for work arrays
587!> \param n3 logical dimension of the FFT transform, reference for work arrays
588!> \param md2 Dimensions of real grid
589!> \param nproc ...
590!> \param nd3 Dimensions of the kernel
591!> \param zw Work array (input)
592!> \param zmpi2 Work array for multiprocessor manipulation (output)
593!> \date February 2006
594!> \author S. Goedecker, L. Genovese
595!> \note
596!> RESTRICTIONS on USAGE
597!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
598!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
599!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
600!> This file is distributed under the terms of the
601!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
602! **************************************************************************************************
603 SUBROUTINE scramble_p(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zw, zmpi2)
604 INTEGER, INTENT(in) :: i1, j2, lot, nfft, n1, n3, md2, nproc, &
605 nd3
606 REAL(kind=dp), DIMENSION(2, lot, n3), INTENT(in) :: zw
607 REAL(kind=dp), DIMENSION(2, n1, md2/nproc, nd3), &
608 INTENT(inout) :: zmpi2
609
610 INTEGER :: i, i3
611
612 DO i3 = 1, n3/2 + 1
613 DO i = 0, nfft - 1
614 zmpi2(1, i1 + i, j2, i3) = zw(1, i + 1, i3)
615 zmpi2(2, i1 + i, j2, i3) = zw(2, i + 1, i3)
616 END DO
617 END DO
618
619 END SUBROUTINE scramble_p
620
621! **************************************************************************************************
622!> \brief (Based on suitable modifications of S.Goedecker routines)
623!> Insert the correct planes of the work array zmpi2
624!> in order to prepare for backward FFT transform
625!> \param i1 Starting points of the plane and number of remaining lines
626!> \param j2 Starting points of the plane and number of remaining lines
627!> \param lot Starting points of the plane and number of remaining lines
628!> \param nfft Starting points of the plane and number of remaining lines
629!> \param n1 logical dimension of the FFT transform, reference for work arrays
630!> \param n3 logical dimension of the FFT transform, reference for work arrays
631!> \param md2 Dimensions of real grid
632!> \param nproc ...
633!> \param nd3 Dimensions of the kernel
634!> \param zmpi2 Work array for multiprocessor manipulation (output)
635!> \param zw Work array (input)
636!> \date February 2006
637!> \author S. Goedecker, L. Genovese
638!> \note
639!> RESTRICTIONS on USAGE
640!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
641!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
642!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
643!> This file is distributed under the terms of the
644!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
645! **************************************************************************************************
646 SUBROUTINE unscramble_p(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zmpi2, zw)
647 INTEGER, INTENT(in) :: i1, j2, lot, nfft, n1, n3, md2, nproc, &
648 nd3
649 REAL(kind=dp), DIMENSION(2, n1, md2/nproc, nd3), &
650 INTENT(in) :: zmpi2
651 REAL(kind=dp), DIMENSION(2, lot, n3), &
652 INTENT(inout) :: zw
653
654 INTEGER :: i, i3, j3
655
656 i3 = 1
657 DO i = 0, nfft - 1
658 zw(1, i + 1, i3) = zmpi2(1, i1 + i, j2, i3)
659 zw(2, i + 1, i3) = zmpi2(2, i1 + i, j2, i3)
660 END DO
661
662 DO i3 = 2, n3/2 + 1
663 j3 = n3 + 2 - i3
664 DO i = 0, nfft - 1
665 zw(1, i + 1, j3) = zmpi2(1, i1 + i, j2, i3)
666 zw(2, i + 1, j3) = -zmpi2(2, i1 + i, j2, i3)
667 zw(1, i + 1, i3) = zmpi2(1, i1 + i, j2, i3)
668 zw(2, i + 1, i3) = zmpi2(2, i1 + i, j2, i3)
669 END DO
670 END DO
671
672 END SUBROUTINE unscramble_p
673
674! **************************************************************************************************
675!> \brief (Based on suitable modifications of S.Goedecker routines)
676!> Multiply with the kernel taking into account its symmetry
677!> Conceived to be used into convolution loops
678!> \param n1 ...
679!> \param n2 ...
680!> \param n3 ...
681!> \param lot ...
682!> \param nfft ...
683!> \param jS ...
684!> \param i3 ...
685!> \param zw Work array (input/output)
686!> n1,n2: logical dimension of the FFT transform, reference for zw
687!> nd1,nd2: Dimensions of POT
688!> jS,j3,nfft: starting point of the plane and number of remaining lines
689!>
690!> RESTRICTIONS on USAGE
691!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
692!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
693!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
694!> This file is distributed under the terms of the
695!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
696!> \param hx ...
697!> \param hy ...
698!> \param hz ...
699!> \date February 2006
700!> \author S. Goedecker, L. Genovese
701! **************************************************************************************************
702 SUBROUTINE p_multkernel(n1, n2, n3, lot, nfft, jS, i3, zw, hx, hy, hz)
703 INTEGER, INTENT(in) :: n1, n2, n3, lot, nfft, js, i3
704 REAL(kind=dp), DIMENSION(2, lot, n2), &
705 INTENT(inout) :: zw
706 REAL(kind=dp), INTENT(in) :: hx, hy, hz
707
708 INTEGER :: i1, i2, j1, j2, j3
709 REAL(kind=dp) :: fourpi2, ker, mu3, p1, p2
710
711 fourpi2 = 4._dp*pi**2
712 j3 = i3 !n3/2+1-abs(n3/2+2-i3)
713 mu3 = real(j3 - 1, kind=dp)/real(n3, kind=dp)
714 mu3 = (mu3/hy)**2 !beware of the exchanged dimension
715 !Body
716 !generic case
717 DO i2 = 1, n2
718 DO i1 = 1, nfft
719 j1 = i1 + js - 1
720 j1 = j1 - (j1/(n1/2 + 2))*n1 !n1/2+1-abs(n1/2+2-jS-i1)
721 j2 = i2 - (i2/(n2/2 + 2))*n2 !n2/2+1-abs(n2/2+1-i2)
722 p1 = real(j1 - 1, kind=dp)/real(n1, kind=dp)
723 p2 = real(j2 - 1, kind=dp)/real(n2, kind=dp)
724 ker = -fourpi2*((p1/hx)**2 + (p2/hz)**2 + mu3) !beware of the exchanged dimension
725 IF (ker /= 0._dp) ker = 1._dp/ker
726 zw(1, i1, i2) = zw(1, i1, i2)*ker
727 zw(2, i1, i2) = zw(2, i1, i2)*ker
728 END DO
729 END DO
730
731 END SUBROUTINE p_multkernel
732
733! **************************************************************************************************
734!> \brief (Based on suitable modifications of S.Goedecker routines)
735!> Multiply with the kernel taking into account its symmetry
736!> Conceived to be used into convolution loops
737!> \param nd1 ...
738!> \param nd2 ...
739!> \param n1 ...
740!> \param n2 ...
741!> \param lot ...
742!> \param nfft ...
743!> \param jS ...
744!> \param pot Kernel, symmetric and real, half the length
745!> \param zw Work array (input/output)
746!> n1,n2: logical dimension of the FFT transform, reference for zw
747!> nd1,nd2: Dimensions of POT
748!> jS, nfft: starting point of the plane and number of remaining lines
749!>
750!> RESTRICTIONS on USAGE
751!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
752!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
753!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
754!> This file is distributed under the terms of the
755!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
756!> \date February 2006
757!> \author S. Goedecker, L. Genovese
758! **************************************************************************************************
759 SUBROUTINE multkernel(nd1, nd2, n1, n2, lot, nfft, jS, pot, zw)
760 INTEGER, INTENT(in) :: nd1, nd2, n1, n2, lot, nfft, js
761 REAL(kind=dp), DIMENSION(nd1, nd2), INTENT(in) :: pot
762 REAL(kind=dp), DIMENSION(2, lot, n2), &
763 INTENT(inout) :: zw
764
765 INTEGER :: i2, j, j1, j2
766
767 DO j = 1, nfft
768 j1 = j + js - 1
769 j1 = j1 + (j1/(n1/2 + 2))*(n1 + 2 - 2*j1)
770 zw(1, j, 1) = zw(1, j, 1)*pot(j1, 1)
771 zw(2, j, 1) = zw(2, j, 1)*pot(j1, 1)
772 END DO
773
774 !generic case
775 DO i2 = 2, n2/2
776 DO j = 1, nfft
777 j1 = j + js - 1
778 j1 = j1 + (j1/(n1/2 + 2))*(n1 + 2 - 2*j1)
779 j2 = n2 + 2 - i2
780 zw(1, j, i2) = zw(1, j, i2)*pot(j1, i2)
781 zw(2, j, i2) = zw(2, j, i2)*pot(j1, i2)
782 zw(1, j, j2) = zw(1, j, j2)*pot(j1, i2)
783 zw(2, j, j2) = zw(2, j, j2)*pot(j1, i2)
784 END DO
785 END DO
786
787 !case i2=n2/2+1
788 DO j = 1, nfft
789 j1 = j + js - 1
790 j1 = j1 + (j1/(n1/2 + 2))*(n1 + 2 - 2*j1)
791 j2 = n2/2 + 1
792 zw(1, j, j2) = zw(1, j, j2)*pot(j1, j2)
793 zw(2, j, j2) = zw(2, j, j2)*pot(j1, j2)
794 END DO
795
796 END SUBROUTINE multkernel
797
798! **************************************************************************************************
799!> \brief !HERE POT MUST BE THE KERNEL (BEWARE THE HALF DIMENSION)
800!> ****h* BigDFT/S_PoissonSolver
801!> (Based on suitable modifications of S.Goedecker routines)
802!> Applies the local FFT space Kernel to the density in Real space.
803!> Does NOT calculate the LDA exchange-correlation terms
804!> \param n1 logical dimension of the transform.
805!> \param n2 logical dimension of the transform.
806!> \param n3 logical dimension of the transform.
807!> \param nd1 Dimension of POT
808!> \param nd2 Dimension of POT
809!> \param nd3 Dimension of POT
810!> \param md1 Dimension of ZF
811!> \param md2 Dimension of ZF
812!> \param md3 Dimension of ZF
813!> \param nproc number of processors used as returned by MPI_COMM_SIZE
814!> \param iproc [0:nproc-1] number of processor as returned by MPI_COMM_RANK
815!> \param pot Kernel, only the distributed part (REAL)
816!> POT(i1,i2,i3)
817!> i1=1,nd1 , i2=1,nd2 , i3=1,nd3/nproc
818!> \param zf Density (input/output)
819!> ZF(i1,i3,i2)
820!> i1=1,md1 , i2=1,md2/nproc , i3=1,md3
821!> \param scal factor of renormalization of the FFT in order to acheve unitarity
822!> and the correct dimension
823!> \param mpi_group ...
824!> \date October 2006
825!> \author S. Goedecker, L. Genovese
826!> \note As transform lengths most products of the prime factors 2,3,5 are allowed.
827!> The detailed table with allowed transform lengths can
828!> be found in subroutine CTRIG
829!> \note
830!> RESTRICTIONS on USAGE
831!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
832!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
833!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
834!> This file is distributed under the terms of the
835!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
836! **************************************************************************************************
837 SUBROUTINE s_poissonsolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, pot, zf, &
838 scal, mpi_group)
839 INTEGER, INTENT(in) :: n1, n2, n3, nd1, nd2, nd3, md1, md2, &
840 md3, nproc, iproc
841 REAL(kind=dp), DIMENSION(nd1, nd2, nd3/nproc), &
842 INTENT(in) :: pot
843 REAL(kind=dp), DIMENSION(md1, md3, md2/nproc), &
844 INTENT(inout) :: zf
845 REAL(kind=dp), INTENT(in) :: scal
846
847 CLASS(mp_comm_type), INTENT(in) :: mpi_group
848
849 CHARACTER(len=*), PARAMETER :: routinen = 'S_PoissonSolver'
850 INTEGER, PARAMETER :: ncache_optimal = 8*1024
851
852 INTEGER :: handle, i, i1, i3, ic1, ic2, ic3, inzee, &
853 j, j2, j2stb, j2stf, j3, jp2stb, &
854 jp2stf, lot, lzt, ma, mb, ncache, nfft
855 INTEGER, ALLOCATABLE, DIMENSION(:) :: after1, after2, after3, before1, &
856 before2, before3, now1, now2, now3
857 REAL(kind=dp) :: twopion
858 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: btrig1, btrig2, btrig3, cosinarr, &
859 ftrig1, ftrig2, ftrig3
860 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zt, zw
861 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: zmpi2
862 REAL(kind=dp), ALLOCATABLE, &
863 DIMENSION(:, :, :, :, :) :: zmpi1
864
865 CALL timeset(routinen, handle)
866 ! check input
867 IF (mod(n3, 2) /= 0) cpabort("Parallel convolution:ERROR:n3")
868 IF (nd1 < n1/2 + 1) cpabort("Parallel convolution:ERROR:nd1")
869 IF (nd2 < n2/2 + 1) cpabort("Parallel convolution:ERROR:nd2")
870 IF (nd3 < n3/2 + 1) cpabort("Parallel convolution:ERROR:nd3")
871 IF (md1 < n1) cpabort("Parallel convolution:ERROR:md1")
872 IF (md2 < n2) cpabort("Parallel convolution:ERROR:md2")
873 IF (md3 < n3/2) cpabort("Parallel convolution:ERROR:md3")
874 IF (mod(nd3, nproc) /= 0) cpabort("Parallel convolution:ERROR:nd3")
875 IF (mod(md2, nproc) /= 0) cpabort("Parallel convolution:ERROR:md2")
876
877 !defining work arrays dimensions
878 ncache = ncache_optimal
879 IF (ncache <= max(n1, n2, n3/2)*4) ncache = max(n1, n2, n3/2)*4
880
881 ! if (timing_flag == 1 .and. iproc ==0) print *,'parallel ncache=',ncache
882
883 lzt = n2
884 IF (mod(n2, 2) == 0) lzt = lzt + 1
885 IF (mod(n2, 4) == 0) lzt = lzt + 1 !maybe this is useless
886
887 !Allocations
888 ALLOCATE (btrig1(2, ctrig_length))
889 ALLOCATE (ftrig1(2, ctrig_length))
890 ALLOCATE (after1(7))
891 ALLOCATE (now1(7))
892 ALLOCATE (before1(7))
893 ALLOCATE (btrig2(2, ctrig_length))
894 ALLOCATE (ftrig2(2, ctrig_length))
895 ALLOCATE (after2(7))
896 ALLOCATE (now2(7))
897 ALLOCATE (before2(7))
898 ALLOCATE (btrig3(2, ctrig_length))
899 ALLOCATE (ftrig3(2, ctrig_length))
900 ALLOCATE (after3(7))
901 ALLOCATE (now3(7))
902 ALLOCATE (before3(7))
903 ALLOCATE (zw(2, ncache/4, 2))
904 ALLOCATE (zt(2, lzt, n1))
905 ALLOCATE (zmpi2(2, n1, md2/nproc, nd3))
906 ALLOCATE (cosinarr(2, n3/2))
907 IF (nproc > 1) THEN
908 ALLOCATE (zmpi1(2, n1, md2/nproc, nd3/nproc, nproc))
909 zmpi1 = 0.0_dp
910 END IF
911 zmpi2 = 0.0_dp
912 zw = 0.0_dp
913 zt = 0.0_dp
914
915 !calculating the FFT work arrays (beware on the HalFFT in n3 dimension)
916 CALL ctrig(n3/2, btrig3, after3, before3, now3, 1, ic3)
917 CALL ctrig(n1, btrig1, after1, before1, now1, 1, ic1)
918 CALL ctrig(n2, btrig2, after2, before2, now2, 1, ic2)
919 DO j = 1, n1
920 ftrig1(1, j) = btrig1(1, j)
921 ftrig1(2, j) = -btrig1(2, j)
922 END DO
923 DO j = 1, n2
924 ftrig2(1, j) = btrig2(1, j)
925 ftrig2(2, j) = -btrig2(2, j)
926 END DO
927 DO j = 1, n3/2
928 ftrig3(1, j) = btrig3(1, j)
929 ftrig3(2, j) = -btrig3(2, j)
930 END DO
931
932 !Calculating array of phases for HalFFT decoding
933 twopion = 8._dp*atan(1._dp)/real(n3, kind=dp)
934 DO i3 = 1, n3/2
935 cosinarr(1, i3) = cos(twopion*(i3 - 1))
936 cosinarr(2, i3) = -sin(twopion*(i3 - 1))
937 END DO
938
939 !initializing integral
940 !ehartree=0._dp
941
942 ! transform along z axis
943 lot = ncache/(2*n3)
944 IF (lot < 1) THEN
945 CALL cp_abort(__location__, &
946 'convolxc_off:ncache has to be enlarged to be able to hold at '// &
947 'least one 1-d FFT of this size even though this will '// &
948 'reduce the performance for shorter transform lengths')
949 END IF
950
951 DO j2 = 1, md2/nproc
952 !this condition ensures that we manage only the interesting part for the FFT
953 IF (iproc*(md2/nproc) + j2 <= n2) THEN
954 DO i1 = 1, n1, lot
955 ma = i1
956 mb = min(i1 + (lot - 1), n1)
957 nfft = mb - ma + 1
958
959 !inserting real data into complex array of half length
960 CALL halfill_upcorn(md1, md3, lot, nfft, n3, zf(i1, 1, j2), zw(1, 1, 1))
961
962 !performing FFT
963 !input: I1,I3,J2,(Jp2)
964 inzee = 1
965 DO i = 1, ic3
966 CALL fftstp(lot, nfft, n3/2, lot, n3/2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
967 btrig3, after3(i), now3(i), before3(i), 1)
968 inzee = 3 - inzee
969 END DO
970 !output: I1,i3,J2,(Jp2)
971 !unpacking FFT in order to restore correct result,
972 !while exchanging components
973 !input: I1,i3,J2,(Jp2)
974 CALL scramble_unpack(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zw(1, 1, inzee), zmpi2, cosinarr)
975 !output: I1,J2,i3,(Jp2)
976 END DO
977 END IF
978 END DO
979 !Interprocessor data transposition
980 !input: I1,J2,j3,jp3,(Jp2)
981 IF (nproc > 1) THEN
982 CALL mpi_group%alltoall(zmpi2, zmpi1, 2*n1*(md2/nproc)*(nd3/nproc))
983 END IF
984 !output: I1,J2,j3,Jp2,(jp3)
985
986 !now each process perform complete convolution of its planes
987 DO j3 = 1, nd3/nproc
988 !this condition ensures that we manage only the interesting part for the FFT
989 IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
990 jp2stb = 1
991 j2stb = 1
992 jp2stf = 1
993 j2stf = 1
994
995 ! transform along x axis
996 lot = ncache/(4*n1)
997 IF (lot < 1) THEN
998 CALL cp_abort(__location__, &
999 'convolxc_off:ncache has to be enlarged to be able to hold at '// &
1000 'least one 1-d FFT of this size even though this will '// &
1001 'reduce the performance for shorter transform lengths')
1002 END IF
1003
1004 DO j = 1, n2, lot
1005 ma = j
1006 mb = min(j + (lot - 1), n2)
1007 nfft = mb - ma + 1
1008
1009 !reverse index ordering, leaving the planes to be transformed at the end
1010 !input: I1,J2,j3,Jp2,(jp3)
1011 IF (nproc == 1) THEN
1012 CALL s_mpiswitch_upcorn(j3, nfft, jp2stb, j2stb, lot, n1, md2, nd3, nproc, zmpi2, zw(1, 1, 1))
1013 ELSE
1014 CALL s_mpiswitch_upcorn(j3, nfft, jp2stb, j2stb, lot, n1, md2, nd3, nproc, zmpi1, zw(1, 1, 1))
1015 END IF
1016 !output: J2,Jp2,I1,j3,(jp3)
1017
1018 !performing FFT
1019 !input: I2,I1,j3,(jp3)
1020 inzee = 1
1021 DO i = 1, ic1 - 1
1022 CALL fftstp(lot, nfft, n1, lot, n1, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1023 btrig1, after1(i), now1(i), before1(i), 1)
1024 inzee = 3 - inzee
1025 END DO
1026
1027 !storing the last step into zt array
1028 i = ic1
1029 CALL fftstp(lot, nfft, n1, lzt, n1, zw(1, 1, inzee), zt(1, j, 1), &
1030 btrig1, after1(i), now1(i), before1(i), 1)
1031 !output: I2,i1,j3,(jp3)
1032 END DO
1033
1034 !transform along y axis
1035 lot = ncache/(4*n2)
1036 IF (lot < 1) THEN
1037 CALL cp_abort(__location__, &
1038 'convolxc_off:ncache has to be enlarged to be able to hold at '// &
1039 'least one 1-d FFT of this size even though this will '// &
1040 'reduce the performance for shorter transform lengths')
1041 END IF
1042
1043 DO j = 1, n1, lot
1044 ma = j
1045 mb = min(j + (lot - 1), n1)
1046 nfft = mb - ma + 1
1047
1048 !reverse ordering
1049 !input: I2,i1,j3,(jp3)
1050 CALL s_switch_upcorn(nfft, n2, lot, n1, lzt, zt(1, 1, j), zw(1, 1, 1))
1051 !output: i1,I2,j3,(jp3)
1052
1053 !performing FFT
1054 !input: i1,I2,j3,(jp3)
1055 inzee = 1
1056 DO i = 1, ic2
1057 CALL fftstp(lot, nfft, n2, lot, n2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1058 btrig2, after2(i), now2(i), before2(i), 1)
1059 inzee = 3 - inzee
1060 END DO
1061 !output: i1,i2,j3,(jp3)
1062
1063 !Multiply with kernel in fourier space
1064 CALL multkernel(nd1, nd2, n1, n2, lot, nfft, j, pot(1, 1, j3), zw(1, 1, inzee))
1065
1066 !TRANSFORM BACK IN REAL SPACE
1067
1068 !transform along y axis
1069 !input: i1,i2,j3,(jp3)
1070 DO i = 1, ic2
1071 CALL fftstp(lot, nfft, n2, lot, n2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1072 ftrig2, after2(i), now2(i), before2(i), -1)
1073 inzee = 3 - inzee
1074 END DO
1075
1076 !reverse ordering
1077 !input: i1,I2,j3,(jp3)
1078 CALL s_unswitch_downcorn(nfft, n2, lot, n1, lzt, zw(1, 1, inzee), zt(1, 1, j))
1079 !output: I2,i1,j3,(jp3)
1080 END DO
1081
1082 !transform along x axis
1083 !input: I2,i1,j3,(jp3)
1084 lot = ncache/(4*n1)
1085 DO j = 1, n2, lot
1086 ma = j
1087 mb = min(j + (lot - 1), n2)
1088 nfft = mb - ma + 1
1089
1090 !performing FFT
1091 i = 1
1092 CALL fftstp(lzt, nfft, n1, lot, n1, zt(1, j, 1), zw(1, 1, 1), &
1093 ftrig1, after1(i), now1(i), before1(i), -1)
1094
1095 inzee = 1
1096 DO i = 2, ic1
1097 CALL fftstp(lot, nfft, n1, lot, n1, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1098 ftrig1, after1(i), now1(i), before1(i), -1)
1099 inzee = 3 - inzee
1100 END DO
1101 !output: I2,I1,j3,(jp3)
1102
1103 !reverse ordering
1104 !input: J2,Jp2,I1,j3,(jp3)
1105 IF (nproc == 1) THEN
1106 CALL s_unmpiswitch_downcorn(j3, nfft, jp2stf, j2stf, lot, n1, md2, nd3, nproc, zw(1, 1, inzee), zmpi2)
1107 ELSE
1108 CALL s_unmpiswitch_downcorn(j3, nfft, jp2stf, j2stf, lot, n1, md2, nd3, nproc, zw(1, 1, inzee), zmpi1)
1109 END IF
1110 ! output: I1,J2,j3,Jp2,(jp3)
1111 END DO
1112 END IF
1113 END DO
1114
1115 !Interprocessor data transposition
1116 !input: I1,J2,j3,Jp2,(jp3)
1117 IF (nproc > 1) THEN
1118 !communication scheduling
1119 CALL mpi_group%alltoall(zmpi1, zmpi2, 2*n1*(md2/nproc)*(nd3/nproc))
1120 END IF
1121
1122 !output: I1,J2,j3,jp3,(Jp2)
1123
1124 !transform along z axis
1125 !input: I1,J2,i3,(Jp2)
1126 lot = ncache/(2*n3)
1127 DO j2 = 1, md2/nproc
1128 !this condition ensures that we manage only the interesting part for the FFT
1129 IF (iproc*(md2/nproc) + j2 <= n2) THEN
1130 DO i1 = 1, n1, lot
1131 ma = i1
1132 mb = min(i1 + (lot - 1), n1)
1133 nfft = mb - ma + 1
1134
1135 !reverse ordering and repack the FFT data in order to be backward HalFFT transformed
1136 !input: I1,J2,i3,(Jp2)
1137 CALL unscramble_pack(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zmpi2, zw(1, 1, 1), cosinarr)
1138 !output: I1,i3,J2,(Jp2)
1139
1140 !performing FFT
1141 !input: I1,i3,J2,(Jp2)
1142 inzee = 1
1143 DO i = 1, ic3
1144 CALL fftstp(lot, nfft, n3/2, lot, n3/2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1145 ftrig3, after3(i), now3(i), before3(i), -1)
1146 inzee = 3 - inzee
1147 END DO
1148 !output: I1,I3,J2,(Jp2)
1149
1150 !rebuild the output array
1151 CALL unfill_downcorn(md1, md3, lot, nfft, n3, zw(1, 1, inzee), zf(i1, 1, j2) &
1152 , scal)
1153
1154 !integrate local pieces together
1155 !ehartree=ehartree+0.5_dp*ehartreetmp*hx*hy*hz
1156 END DO
1157 END IF
1158 END DO
1159
1160 !De-allocations
1161 DEALLOCATE (btrig1)
1162 DEALLOCATE (ftrig1)
1163 DEALLOCATE (after1)
1164 DEALLOCATE (now1)
1165 DEALLOCATE (before1)
1166 DEALLOCATE (btrig2)
1167 DEALLOCATE (ftrig2)
1168 DEALLOCATE (after2)
1169 DEALLOCATE (now2)
1170 DEALLOCATE (before2)
1171 DEALLOCATE (btrig3)
1172 DEALLOCATE (ftrig3)
1173 DEALLOCATE (after3)
1174 DEALLOCATE (now3)
1175 DEALLOCATE (before3)
1176 DEALLOCATE (zmpi2)
1177 DEALLOCATE (zw)
1178 DEALLOCATE (zt)
1179 DEALLOCATE (cosinarr)
1180 IF (nproc > 1) DEALLOCATE (zmpi1)
1181
1182 ! call timing(iproc,'PSolv_comput ','OF')
1183 CALL timestop(handle)
1184 END SUBROUTINE s_poissonsolver
1185
1186! **************************************************************************************************
1187!> \brief ...
1188!> \param j3 ...
1189!> \param nfft ...
1190!> \param Jp2stb ...
1191!> \param J2stb ...
1192!> \param lot ...
1193!> \param n1 ...
1194!> \param md2 ...
1195!> \param nd3 ...
1196!> \param nproc ...
1197!> \param zmpi1 ...
1198!> \param zw ...
1199! **************************************************************************************************
1200 SUBROUTINE s_mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot, n1, md2, nd3, nproc, zmpi1, zw)
1201 INTEGER, INTENT(in) :: j3, nfft
1202 INTEGER, INTENT(inout) :: jp2stb, j2stb
1203 INTEGER, INTENT(in) :: lot, n1, md2, nd3, nproc
1204 REAL(kind=dp), &
1205 DIMENSION(2, n1, md2/nproc, nd3/nproc, nproc), &
1206 INTENT(in) :: zmpi1
1207 REAL(kind=dp), DIMENSION(2, lot, n1), &
1208 INTENT(inout) :: zw
1209
1210 INTEGER :: i1, j2, jp2, mfft
1211
1212 mfft = 0
1213 DO jp2 = jp2stb, nproc
1214 DO j2 = j2stb, md2/nproc
1215 mfft = mfft + 1
1216 IF (mfft > nfft) THEN
1217 jp2stb = jp2
1218 j2stb = j2
1219 RETURN
1220 END IF
1221 DO i1 = 1, n1
1222 zw(1, mfft, i1) = zmpi1(1, i1, j2, j3, jp2)
1223 zw(2, mfft, i1) = zmpi1(2, i1, j2, j3, jp2)
1224 END DO
1225 END DO
1226 j2stb = 1
1227 END DO
1228 END SUBROUTINE s_mpiswitch_upcorn
1229
1230! **************************************************************************************************
1231!> \brief ...
1232!> \param nfft ...
1233!> \param n2 ...
1234!> \param lot ...
1235!> \param n1 ...
1236!> \param lzt ...
1237!> \param zt ...
1238!> \param zw ...
1239! **************************************************************************************************
1240 SUBROUTINE s_switch_upcorn(nfft, n2, lot, n1, lzt, zt, zw)
1241 INTEGER, INTENT(in) :: nfft, n2, lot, n1, lzt
1242 REAL(kind=dp), DIMENSION(2, lzt, n1), INTENT(in) :: zt
1243 REAL(kind=dp), DIMENSION(2, lot, n2), &
1244 INTENT(inout) :: zw
1245
1246 INTEGER :: i, j
1247
1248 DO j = 1, nfft
1249 DO i = 1, n2
1250 zw(1, j, i) = zt(1, i, j)
1251 zw(2, j, i) = zt(2, i, j)
1252 END DO
1253 END DO
1254 END SUBROUTINE s_switch_upcorn
1255
1256! **************************************************************************************************
1257!> \brief ...
1258!> \param nfft ...
1259!> \param n2 ...
1260!> \param lot ...
1261!> \param n1 ...
1262!> \param lzt ...
1263!> \param zw ...
1264!> \param zt ...
1265! **************************************************************************************************
1266 SUBROUTINE s_unswitch_downcorn(nfft, n2, lot, n1, lzt, zw, zt)
1267 INTEGER, INTENT(in) :: nfft, n2, lot, n1, lzt
1268 REAL(kind=dp), DIMENSION(2, lot, n2), INTENT(in) :: zw
1269 REAL(kind=dp), DIMENSION(2, lzt, n1), &
1270 INTENT(inout) :: zt
1271
1272 INTEGER :: i, j
1273
1274 DO j = 1, nfft
1275 DO i = 1, n2
1276 zt(1, i, j) = zw(1, j, i)
1277 zt(2, i, j) = zw(2, j, i)
1278 END DO
1279 END DO
1280 END SUBROUTINE s_unswitch_downcorn
1281
1282! **************************************************************************************************
1283!> \brief ...
1284!> \param j3 ...
1285!> \param nfft ...
1286!> \param Jp2stf ...
1287!> \param J2stf ...
1288!> \param lot ...
1289!> \param n1 ...
1290!> \param md2 ...
1291!> \param nd3 ...
1292!> \param nproc ...
1293!> \param zw ...
1294!> \param zmpi1 ...
1295! **************************************************************************************************
1296 SUBROUTINE s_unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot, n1, md2, nd3, nproc, zw, zmpi1)
1297 INTEGER, INTENT(in) :: j3, nfft
1298 INTEGER, INTENT(inout) :: jp2stf, j2stf
1299 INTEGER, INTENT(in) :: lot, n1, md2, nd3, nproc
1300 REAL(kind=dp), DIMENSION(2, lot, n1), INTENT(in) :: zw
1301 REAL(kind=dp), &
1302 DIMENSION(2, n1, md2/nproc, nd3/nproc, nproc), &
1303 INTENT(inout) :: zmpi1
1304
1305 INTEGER :: i1, j2, jp2, mfft
1306
1307 mfft = 0
1308 DO jp2 = jp2stf, nproc
1309 DO j2 = j2stf, md2/nproc
1310 mfft = mfft + 1
1311 IF (mfft > nfft) THEN
1312 jp2stf = jp2
1313 j2stf = j2
1314 RETURN
1315 END IF
1316 DO i1 = 1, n1
1317 zmpi1(1, i1, j2, j3, jp2) = zw(1, mfft, i1)
1318 zmpi1(2, i1, j2, j3, jp2) = zw(2, mfft, i1)
1319 END DO
1320 END DO
1321 j2stf = 1
1322 END DO
1323 END SUBROUTINE s_unmpiswitch_downcorn
1324
1325! **************************************************************************************************
1326!> \brief (Based on suitable modifications of S.Goedecker routines)
1327!> Restore data into output array, calculating in the meanwhile
1328!> Hartree energy of the potential
1329!> \param md1 Dimensions of the undistributed part of the real grid
1330!> \param md3 Dimensions of the undistributed part of the real grid
1331!> \param lot ...
1332!> \param nfft number of planes
1333!> \param n3 (twice the) dimension of the last FFTtransform.
1334!> \param zw FFT work array
1335!> \param zf Original distributed density as well as
1336!> Distributed solution of the poisson equation (inout)
1337!> \param scal Needed to achieve unitarity and correct dimensions
1338!> \date February 2006
1339!> \author S. Goedecker, L. Genovese
1340!> \note Assuming that high frequencies are in the corners
1341!> and that n3 is multiple of 4
1342!>
1343!> RESTRICTIONS on USAGE
1344!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1345!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1346!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
1347!> This file is distributed under the terms of the
1348!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
1349! **************************************************************************************************
1350 SUBROUTINE unfill_downcorn(md1, md3, lot, nfft, n3, zw, zf, scal)
1351 INTEGER, INTENT(in) :: md1, md3, lot, nfft, n3
1352 REAL(kind=dp), DIMENSION(2, lot, n3/2), INTENT(in) :: zw
1353 REAL(kind=dp), DIMENSION(md1, md3), INTENT(inout) :: zf
1354 REAL(kind=dp), INTENT(in) :: scal
1355
1356 INTEGER :: i1, i3
1357 REAL(kind=dp) :: pot1
1358
1359 DO i3 = 1, n3/4
1360 DO i1 = 1, nfft
1361 pot1 = scal*zw(1, i1, i3)
1362 !ehartreetmp =ehartreetmp + pot1* zf(i1,2*i3-1)
1363 zf(i1, 2*i3 - 1) = pot1
1364 pot1 = scal*zw(2, i1, i3)
1365 !ehartreetmp =ehartreetmp + pot1* zf(i1,2*i3)
1366 zf(i1, 2*i3) = pot1
1367 END DO
1368 END DO
1369 END SUBROUTINE unfill_downcorn
1370
1371! **************************************************************************************************
1372!> \brief ...
1373!> \param md1 ...
1374!> \param md3 ...
1375!> \param lot ...
1376!> \param nfft ...
1377!> \param n3 ...
1378!> \param zf ...
1379!> \param zw ...
1380! **************************************************************************************************
1381 SUBROUTINE halfill_upcorn(md1, md3, lot, nfft, n3, zf, zw)
1382 INTEGER :: md1, md3, lot, nfft, n3
1383 REAL(kind=dp) :: zf(md1, md3), zw(2, lot, n3/2)
1384
1385 INTEGER :: i1, i3
1386
1387 DO i3 = 1, n3/4
1388 ! WARNING: Assuming that high frequencies are in the corners
1389 ! and that n3 is multiple of 4
1390 !in principle we can relax this condition
1391 DO i1 = 1, nfft
1392 zw(1, i1, i3) = 0._dp
1393 zw(2, i1, i3) = 0._dp
1394 END DO
1395 END DO
1396 DO i3 = n3/4 + 1, n3/2
1397 DO i1 = 1, nfft
1398 zw(1, i1, i3) = zf(i1, 2*i3 - 1 - n3/2)
1399 zw(2, i1, i3) = zf(i1, 2*i3 - n3/2)
1400 END DO
1401 END DO
1402
1403 END SUBROUTINE halfill_upcorn
1404
1405! **************************************************************************************************
1406!> \brief (Based on suitable modifications of S.Goedecker routines)
1407!> Assign the correct planes to the work array zmpi2
1408!> in order to prepare for interprocessor data transposition.
1409!> In the meanwhile, it unpacks the data of the HalFFT in order to prepare for
1410!> multiplication with the kernel
1411!> \param i1 Starting points of the plane and number of remaining lines
1412!> \param j2 Starting points of the plane and number of remaining lines
1413!> \param lot Starting points of the plane and number of remaining lines
1414!> \param nfft Starting points of the plane and number of remaining lines
1415!> \param n1 logical dimension of the FFT transform, reference for work arrays
1416!> \param n3 logical dimension of the FFT transform, reference for work arrays
1417!> \param md2 Dimensions of real grid
1418!> \param nproc ...
1419!> \param nd3 Dimensions of the kernel
1420!> \param zw Work array (input)
1421!> \param zmpi2 Work array for multiprocessor manipulation (output)
1422!> \param cosinarr Array of the phases needed for unpacking
1423!> \date February 2006
1424!> \author S. Goedecker, L. Genovese
1425!> \note
1426!> RESTRICTIONS on USAGE
1427!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1428!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1429!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
1430!> This file is distributed under the terms of the
1431!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
1432! **************************************************************************************************
1433 SUBROUTINE scramble_unpack(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zw, zmpi2, cosinarr)
1434 INTEGER, INTENT(in) :: i1, j2, lot, nfft, n1, n3, md2, nproc, &
1435 nd3
1436 REAL(kind=dp), DIMENSION(2, lot, n3/2), INTENT(in) :: zw
1437 REAL(kind=dp), DIMENSION(2, n1, md2/nproc, nd3), &
1438 INTENT(inout) :: zmpi2
1439 REAL(kind=dp), DIMENSION(2, n3/2), INTENT(in) :: cosinarr
1440
1441 INTEGER :: i, i3, ind1, ind2
1442 REAL(kind=dp) :: a, b, c, cp, d, fei, fer, fi, foi, for, &
1443 fr, sp
1444
1445!case i3=1 and i3=n3/2+1
1446
1447 DO i = 0, nfft - 1
1448 a = zw(1, i + 1, 1)
1449 b = zw(2, i + 1, 1)
1450 zmpi2(1, i1 + i, j2, 1) = a + b
1451 zmpi2(2, i1 + i, j2, 1) = 0._dp
1452 zmpi2(1, i1 + i, j2, n3/2 + 1) = a - b
1453 zmpi2(2, i1 + i, j2, n3/2 + 1) = 0._dp
1454 END DO
1455 !case 2<=i3<=n3/2
1456 DO i3 = 2, n3/2
1457 ind1 = i3
1458 ind2 = n3/2 - i3 + 2
1459 cp = cosinarr(1, i3)
1460 sp = cosinarr(2, i3)
1461 DO i = 0, nfft - 1
1462 a = zw(1, i + 1, ind1)
1463 b = zw(2, i + 1, ind1)
1464 c = zw(1, i + 1, ind2)
1465 d = zw(2, i + 1, ind2)
1466 fer = .5_dp*(a + c)
1467 fei = .5_dp*(b - d)
1468 for = .5_dp*(a - c)
1469 foi = .5_dp*(b + d)
1470 fr = fer + cp*foi - sp*for
1471 fi = fei - cp*for - sp*foi
1472 zmpi2(1, i1 + i, j2, ind1) = fr
1473 zmpi2(2, i1 + i, j2, ind1) = fi
1474 END DO
1475 END DO
1476
1477 END SUBROUTINE scramble_unpack
1478
1479! **************************************************************************************************
1480!> \brief (Based on suitable modifications of S.Goedecker routines)
1481!> Insert the correct planes of the work array zmpi2
1482!> in order to prepare for backward FFT transform
1483!> In the meanwhile, it packs the data in order to be transformed with the HalFFT
1484!> procedure
1485!> \param i1 Starting points of the plane and number of remaining lines
1486!> \param j2 Starting points of the plane and number of remaining lines
1487!> \param lot Starting points of the plane and number of remaining lines
1488!> \param nfft Starting points of the plane and number of remaining lines
1489!> \param n1 logical dimension of the FFT transform, reference for work arrays
1490!> \param n3 logical dimension of the FFT transform, reference for work arrays
1491!> \param md2 Dimensions of real grid
1492!> \param nproc ...
1493!> \param nd3 Dimensions of the kernel
1494!> \param zmpi2 Work array for multiprocessor manipulation (output)
1495!> \param zw Work array (inout)
1496!> \param cosinarr Array of the phases needed for packing
1497!> \date February 2006
1498!> \author S. Goedecker, L. Genovese
1499!> \note
1500!> RESTRICTIONS on USAGE
1501!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1502!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1503!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
1504!> This file is distributed under the terms of the
1505!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
1506! **************************************************************************************************
1507 SUBROUTINE unscramble_pack(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zmpi2, zw, cosinarr)
1508 INTEGER, INTENT(in) :: i1, j2, lot, nfft, n1, n3, md2, nproc, &
1509 nd3
1510 REAL(kind=dp), DIMENSION(2, n1, md2/nproc, nd3), &
1511 INTENT(in) :: zmpi2
1512 REAL(kind=dp), DIMENSION(2, lot, n3/2), &
1513 INTENT(inout) :: zw
1514 REAL(kind=dp), DIMENSION(2, n3/2), INTENT(in) :: cosinarr
1515
1516 INTEGER :: i, i3, inda, indb
1517 REAL(kind=dp) :: a, b, c, cp, d, ie, ih, io, re, rh, ro, &
1518 sp
1519
1520 DO i3 = 1, n3/2
1521 inda = i3
1522 indb = n3/2 + 2 - i3
1523 cp = cosinarr(1, i3)
1524 sp = cosinarr(2, i3)
1525 DO i = 0, nfft - 1
1526 a = zmpi2(1, i1 + i, j2, inda)
1527 b = zmpi2(2, i1 + i, j2, inda)
1528 c = zmpi2(1, i1 + i, j2, indb)
1529 d = -zmpi2(2, i1 + i, j2, indb)
1530 re = (a + c)
1531 ie = (b + d)
1532 ro = (a - c)*cp - (b - d)*sp
1533 io = (a - c)*sp + (b - d)*cp
1534 rh = re - io
1535 ih = ie + ro
1536 zw(1, i + 1, inda) = rh
1537 zw(2, i + 1, inda) = ih
1538 END DO
1539 END DO
1540
1541 END SUBROUTINE unscramble_pack
1542
1543! **************************************************************************************************
1544!> \brief (Based on suitable modifications of S.Goedecker routines)
1545!> Applies the local FFT space Kernel to the density in Real space.
1546!> Calculates also the LDA exchange-correlation terms
1547!> \param n1 logical dimension of the transform.
1548!> \param n2 logical dimension of the transform.
1549!> \param n3 logical dimension of the transform.
1550!> \param nd1 Dimension of POT
1551!> \param nd2 Dimension of POT
1552!> \param nd3 Dimension of POT
1553!> \param md1 Dimension of ZF
1554!> \param md2 Dimension of ZF
1555!> \param md3 Dimension of ZF
1556!> \param nproc number of processors used as returned by MPI_COMM_SIZE
1557!> \param iproc [0:nproc-1] number of processor as returned by MPI_COMM_RANK
1558!> \param pot Kernel, only the distributed part (REAL)
1559!> POT(i1,i2,i3)
1560!> i1=1,nd1 , i2=1,nd2 , i3=1,nd3/nproc
1561!> \param zf Density (input/output)
1562!> ZF(i1,i3,i2)
1563!> i1=1,md1 , i2=1,md2/nproc , i3=1,md3
1564!> \param scal factor of renormalization of the FFT in order to acheve unitarity
1565!> and the correct dimension
1566!> \param mpi_group ...
1567!> \date February 2006
1568!> \author S. Goedecker, L. Genovese
1569!> \note As transform lengths
1570!> most products of the prime factors 2,3,5 are allowed.
1571!> The detailed table with allowed transform lengths can
1572!> be found in subroutine CTRIG
1573!> \note
1574!> RESTRICTIONS on USAGE
1575!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1576!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1577!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
1578!> This file is distributed under the terms of the
1579!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
1580! **************************************************************************************************
1581 SUBROUTINE f_poissonsolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, pot, zf, &
1582 scal, mpi_group)
1583 INTEGER, INTENT(in) :: n1, n2, n3, nd1, nd2, nd3, md1, md2, &
1584 md3, nproc, iproc
1585 REAL(kind=dp), DIMENSION(nd1, nd2, nd3/nproc), &
1586 INTENT(in) :: pot
1587 REAL(kind=dp), DIMENSION(md1, md3, md2/nproc), &
1588 INTENT(inout) :: zf
1589 REAL(kind=dp), INTENT(in) :: scal
1590
1591 CLASS(mp_comm_type), INTENT(in) :: mpi_group
1592
1593 INTEGER, PARAMETER :: ncache_optimal = 8*1024
1594
1595 INTEGER :: i, i1, i3, ic1, ic2, ic3, inzee, j, j2, &
1596 j2stb, j2stf, j3, jp2stb, jp2stf, lot, &
1597 lzt, ma, mb, ncache, nfft
1598 INTEGER, ALLOCATABLE, DIMENSION(:) :: after1, after2, after3, before1, &
1599 before2, before3, now1, now2, now3
1600 REAL(kind=dp) :: twopion
1601 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: btrig1, btrig2, btrig3, cosinarr, &
1602 ftrig1, ftrig2, ftrig3
1603 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zt, zw
1604 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: zmpi2
1605 REAL(kind=dp), ALLOCATABLE, &
1606 DIMENSION(:, :, :, :, :) :: zmpi1
1607
1608 IF (mod(n1, 2) /= 0) cpabort("Parallel convolution:ERROR:n1")
1609 IF (mod(n2, 2) /= 0) cpabort("Parallel convolution:ERROR:n2")
1610 IF (mod(n3, 2) /= 0) cpabort("Parallel convolution:ERROR:n3")
1611 IF (nd1 < n1/2 + 1) cpabort("Parallel convolution:ERROR:nd1")
1612 IF (nd2 < n2/2 + 1) cpabort("Parallel convolution:ERROR:nd2")
1613 IF (nd3 < n3/2 + 1) cpabort("Parallel convolution:ERROR:nd3")
1614 IF (md1 < n1/2) cpabort("Parallel convolution:ERROR:md1")
1615 IF (md2 < n2/2) cpabort("Parallel convolution:ERROR:md2")
1616 IF (md3 < n3/2) cpabort("Parallel convolution:ERROR:md3")
1617 IF (mod(nd3, nproc) /= 0) cpabort("Parallel convolution:ERROR:nd3")
1618 IF (mod(md2, nproc) /= 0) cpabort("Parallel convolution:ERROR:md2")
1619
1620 !defining work arrays dimensions
1621
1622 ncache = ncache_optimal
1623 IF (ncache <= max(n1, n2, n3/2)*4) ncache = max(n1, n2, n3/2)*4
1624 lzt = n2/2
1625 IF (mod(n2/2, 2) == 0) lzt = lzt + 1
1626 IF (mod(n2/2, 4) == 0) lzt = lzt + 1
1627
1628 !Allocations
1629 ALLOCATE (btrig1(2, ctrig_length))
1630 ALLOCATE (ftrig1(2, ctrig_length))
1631 ALLOCATE (after1(7))
1632 ALLOCATE (now1(7))
1633 ALLOCATE (before1(7))
1634 ALLOCATE (btrig2(2, ctrig_length))
1635 ALLOCATE (ftrig2(2, ctrig_length))
1636 ALLOCATE (after2(7))
1637 ALLOCATE (now2(7))
1638 ALLOCATE (before2(7))
1639 ALLOCATE (btrig3(2, ctrig_length))
1640 ALLOCATE (ftrig3(2, ctrig_length))
1641 ALLOCATE (after3(7))
1642 ALLOCATE (now3(7))
1643 ALLOCATE (before3(7))
1644 ALLOCATE (zw(2, ncache/4, 2))
1645 ALLOCATE (zt(2, lzt, n1))
1646 ALLOCATE (zmpi2(2, n1, md2/nproc, nd3))
1647 zmpi2 = 1.0e99_dp ! init mpi'ed buffers to silence warnings under valgrind
1648 ALLOCATE (cosinarr(2, n3/2))
1649 IF (nproc > 1) ALLOCATE (zmpi1(2, n1, md2/nproc, nd3/nproc, nproc))
1650
1651 !calculating the FFT work arrays (beware on the HalFFT in n3 dimension)
1652 CALL ctrig(n3/2, btrig3, after3, before3, now3, 1, ic3)
1653 CALL ctrig(n1, btrig1, after1, before1, now1, 1, ic1)
1654 CALL ctrig(n2, btrig2, after2, before2, now2, 1, ic2)
1655 DO j = 1, n1
1656 ftrig1(1, j) = btrig1(1, j)
1657 ftrig1(2, j) = -btrig1(2, j)
1658 END DO
1659 DO j = 1, n2
1660 ftrig2(1, j) = btrig2(1, j)
1661 ftrig2(2, j) = -btrig2(2, j)
1662 END DO
1663 DO j = 1, n3
1664 ftrig3(1, j) = btrig3(1, j)
1665 ftrig3(2, j) = -btrig3(2, j)
1666 END DO
1667
1668 !Calculating array of phases for HalFFT decoding
1669 twopion = 8._dp*atan(1._dp)/real(n3, kind=dp)
1670 DO i3 = 1, n3/2
1671 cosinarr(1, i3) = cos(twopion*(i3 - 1))
1672 cosinarr(2, i3) = -sin(twopion*(i3 - 1))
1673 END DO
1674
1675 ! transform along z axis
1676 lot = ncache/(2*n3)
1677 IF (lot < 1) THEN
1678 CALL cp_abort(__location__, &
1679 'convolxc_on:ncache has to be enlarged to be able to hold at '// &
1680 'least one 1-d FFT of this size even though this will '// &
1681 'reduce the performance for shorter transform lengths')
1682 END IF
1683
1684 DO j2 = 1, md2/nproc
1685 !this condition ensures that we manage only the interesting part for the FFT
1686 IF (iproc*(md2/nproc) + j2 <= n2/2) THEN
1687 DO i1 = 1, (n1/2), lot
1688 ma = i1
1689 mb = min(i1 + (lot - 1), (n1/2))
1690 nfft = mb - ma + 1
1691
1692 !inserting real data into complex array of half length
1693 CALL halfill_upcorn(md1, md3, lot, nfft, n3, zf(i1, 1, j2), zw(1, 1, 1))
1694
1695 !performing FFT
1696 !input: I1,I3,J2,(Jp2)
1697 inzee = 1
1698 DO i = 1, ic3
1699 CALL fftstp(lot, nfft, n3/2, lot, n3/2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1700 btrig3, after3(i), now3(i), before3(i), 1)
1701 inzee = 3 - inzee
1702 END DO
1703 !output: I1,i3,J2,(Jp2)
1704
1705 !unpacking FFT in order to restore correct result,
1706 !while exchanging components
1707 !input: I1,i3,J2,(Jp2)
1708 CALL scramble_unpack(i1, j2, lot, nfft, n1/2, n3, md2, nproc, nd3, zw(1, 1, inzee), zmpi2, cosinarr)
1709 !output: I1,J2,i3,(Jp2)
1710 END DO
1711 END IF
1712 END DO
1713
1714 !Interprocessor data transposition
1715 !input: I1,J2,j3,jp3,(Jp2)
1716 IF (nproc > 1) THEN
1717 !communication scheduling
1718 CALL mpi_group%alltoall(zmpi2, zmpi1, n1*(md2/nproc)*(nd3/nproc))
1719 END IF
1720 !output: I1,J2,j3,Jp2,(jp3)
1721
1722 !now each process perform complete convolution of its planes
1723 DO j3 = 1, nd3/nproc
1724 !this condition ensures that we manage only the interesting part for the FFT
1725 IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
1726 jp2stb = 1
1727 j2stb = 1
1728 jp2stf = 1
1729 j2stf = 1
1730
1731 ! transform along x axis
1732 lot = ncache/(4*n1)
1733 IF (lot < 1) THEN
1734 CALL cp_abort(__location__, &
1735 'convolxc_on:ncache has to be enlarged to be able to hold at '// &
1736 'least one 1-d FFT of this size even though this will '// &
1737 'reduce the performance for shorter transform lengths')
1738 END IF
1739
1740 DO j = 1, n2/2, lot
1741 ma = j
1742 mb = min(j + (lot - 1), n2/2)
1743 nfft = mb - ma + 1
1744
1745 !reverse index ordering, leaving the planes to be transformed at the end
1746 !input: I1,J2,j3,Jp2,(jp3)
1747 IF (nproc == 1) THEN
1748 CALL mpiswitch_upcorn(j3, nfft, jp2stb, j2stb, lot, n1, md2, nd3, nproc, zmpi2, zw(1, 1, 1))
1749 ELSE
1750 CALL mpiswitch_upcorn(j3, nfft, jp2stb, j2stb, lot, n1, md2, nd3, nproc, zmpi1, zw(1, 1, 1))
1751 END IF
1752 !output: J2,Jp2,I1,j3,(jp3)
1753
1754 !performing FFT
1755 !input: I2,I1,j3,(jp3)
1756 inzee = 1
1757 DO i = 1, ic1 - 1
1758 CALL fftstp(lot, nfft, n1, lot, n1, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1759 btrig1, after1(i), now1(i), before1(i), 1)
1760 inzee = 3 - inzee
1761 END DO
1762
1763 !storing the last step into zt array
1764 i = ic1
1765 CALL fftstp(lot, nfft, n1, lzt, n1, zw(1, 1, inzee), zt(1, j, 1), &
1766 btrig1, after1(i), now1(i), before1(i), 1)
1767 !output: I2,i1,j3,(jp3)
1768 END DO
1769
1770 !transform along y axis
1771 lot = ncache/(4*n2)
1772 IF (lot < 1) THEN
1773 CALL cp_abort(__location__, &
1774 'convolxc_on:ncache has to be enlarged to be able to hold at '// &
1775 'least one 1-d FFT of this size even though this will '// &
1776 'reduce the performance for shorter transform lengths')
1777 END IF
1778
1779 DO j = 1, n1, lot
1780 ma = j
1781 mb = min(j + (lot - 1), n1)
1782 nfft = mb - ma + 1
1783
1784 !reverse ordering
1785 !input: I2,i1,j3,(jp3)
1786 CALL switch_upcorn(nfft, n2, lot, n1, lzt, zt(1, 1, j), zw(1, 1, 1))
1787 !output: i1,I2,j3,(jp3)
1788
1789 !performing FFT
1790 !input: i1,I2,j3,(jp3)
1791 inzee = 1
1792 DO i = 1, ic2
1793 CALL fftstp(lot, nfft, n2, lot, n2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1794 btrig2, after2(i), now2(i), before2(i), 1)
1795 inzee = 3 - inzee
1796 END DO
1797 !output: i1,i2,j3,(jp3)
1798
1799 !Multiply with kernel in fourier space
1800 CALL multkernel(nd1, nd2, n1, n2, lot, nfft, j, pot(1, 1, j3), zw(1, 1, inzee))
1801
1802 !TRANSFORM BACK IN REAL SPACE
1803
1804 !transform along y axis
1805 !input: i1,i2,j3,(jp3)
1806 DO i = 1, ic2
1807 CALL fftstp(lot, nfft, n2, lot, n2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1808 ftrig2, after2(i), now2(i), before2(i), -1)
1809 inzee = 3 - inzee
1810 END DO
1811
1812 !reverse ordering
1813 !input: i1,I2,j3,(jp3)
1814 CALL unswitch_downcorn(nfft, n2, lot, n1, lzt, zw(1, 1, inzee), zt(1, 1, j))
1815 !output: I2,i1,j3,(jp3)
1816 END DO
1817
1818 !transform along x axis
1819 !input: I2,i1,j3,(jp3)
1820 lot = ncache/(4*n1)
1821 DO j = 1, n2/2, lot
1822 ma = j
1823 mb = min(j + (lot - 1), n2/2)
1824 nfft = mb - ma + 1
1825
1826 !performing FFT
1827 i = 1
1828 CALL fftstp(lzt, nfft, n1, lot, n1, zt(1, j, 1), zw(1, 1, 1), &
1829 ftrig1, after1(i), now1(i), before1(i), -1)
1830
1831 inzee = 1
1832 DO i = 2, ic1
1833 CALL fftstp(lot, nfft, n1, lot, n1, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1834 ftrig1, after1(i), now1(i), before1(i), -1)
1835 inzee = 3 - inzee
1836 END DO
1837 !output: I2,I1,j3,(jp3)
1838
1839 !reverse ordering
1840 !input: J2,Jp2,I1,j3,(jp3)
1841 IF (nproc == 1) THEN
1842 CALL unmpiswitch_downcorn(j3, nfft, jp2stf, j2stf, lot, n1, md2, nd3, nproc, zw(1, 1, inzee), zmpi2)
1843 ELSE
1844 CALL unmpiswitch_downcorn(j3, nfft, jp2stf, j2stf, lot, n1, md2, nd3, nproc, zw(1, 1, inzee), zmpi1)
1845 END IF
1846 ! output: I1,J2,j3,Jp2,(jp3)
1847 END DO
1848 END IF
1849 END DO
1850
1851 !Interprocessor data transposition
1852 !input: I1,J2,j3,Jp2,(jp3)
1853 IF (nproc > 1) THEN
1854 !communication scheduling
1855 CALL mpi_group%alltoall(zmpi1, zmpi2, n1*(md2/nproc)*(nd3/nproc))
1856 !output: I1,J2,j3,jp3,(Jp2)
1857 END IF
1858
1859 !transform along z axis
1860 !input: I1,J2,i3,(Jp2)
1861 lot = ncache/(2*n3)
1862 DO j2 = 1, md2/nproc
1863 !this condition ensures that we manage only the interesting part for the FFT
1864 IF (iproc*(md2/nproc) + j2 <= n2/2) THEN
1865 DO i1 = 1, (n1/2), lot
1866 ma = i1
1867 mb = min(i1 + (lot - 1), (n1/2))
1868 nfft = mb - ma + 1
1869
1870 !reverse ordering and repack the FFT data in order to be backward HalFFT transformed
1871 !input: I1,J2,i3,(Jp2)
1872 CALL unscramble_pack(i1, j2, lot, nfft, n1/2, n3, md2, nproc, nd3, zmpi2, zw(1, 1, 1), cosinarr)
1873 !output: I1,i3,J2,(Jp2)
1874
1875 !performing FFT
1876 !input: I1,i3,J2,(Jp2)
1877 inzee = 1
1878 DO i = 1, ic3
1879 CALL fftstp(lot, nfft, n3/2, lot, n3/2, zw(1, 1, inzee), zw(1, 1, 3 - inzee), &
1880 ftrig3, after3(i), now3(i), before3(i), -1)
1881 inzee = 3 - inzee
1882 END DO
1883 !output: I1,I3,J2,(Jp2)
1884
1885 !calculates the exchange correlation terms locally and rebuild the output array
1886 CALL unfill_downcorn(md1, md3, lot, nfft, n3, zw(1, 1, inzee), zf(i1, 1, j2) &
1887 , scal) !,ehartreetmp)
1888
1889 !integrate local pieces together
1890 !ehartree=ehartree+0.5_dp*ehartreetmp*hgrid**3
1891 END DO
1892 END IF
1893 END DO
1894
1895 !De-allocations
1896 DEALLOCATE (btrig1)
1897 DEALLOCATE (ftrig1)
1898 DEALLOCATE (after1)
1899 DEALLOCATE (now1)
1900 DEALLOCATE (before1)
1901 DEALLOCATE (btrig2)
1902 DEALLOCATE (ftrig2)
1903 DEALLOCATE (after2)
1904 DEALLOCATE (now2)
1905 DEALLOCATE (before2)
1906 DEALLOCATE (btrig3)
1907 DEALLOCATE (ftrig3)
1908 DEALLOCATE (after3)
1909 DEALLOCATE (now3)
1910 DEALLOCATE (before3)
1911 DEALLOCATE (zmpi2)
1912 DEALLOCATE (zw)
1913 DEALLOCATE (zt)
1914 DEALLOCATE (cosinarr)
1915 IF (nproc > 1) DEALLOCATE (zmpi1)
1916
1917 END SUBROUTINE f_poissonsolver
1918
1919! **************************************************************************************************
1920!> \brief ...
1921!> \param nfft ...
1922!> \param n2 ...
1923!> \param lot ...
1924!> \param n1 ...
1925!> \param lzt ...
1926!> \param zt ...
1927!> \param zw ...
1928! **************************************************************************************************
1929 SUBROUTINE switch_upcorn(nfft, n2, lot, n1, lzt, zt, zw)
1930 INTEGER :: nfft, n2, lot, n1, lzt
1931 REAL(kind=dp) :: zt(2, lzt, n1), zw(2, lot, n2)
1932
1933 INTEGER :: i, j
1934
1935! WARNING: Assuming that high frequencies are in the corners
1936! and that n2 is multiple of 2
1937! Low frequencies
1938
1939 DO j = 1, nfft
1940 DO i = n2/2 + 1, n2
1941 zw(1, j, i) = zt(1, i - n2/2, j)
1942 zw(2, j, i) = zt(2, i - n2/2, j)
1943 END DO
1944 END DO
1945 ! High frequencies
1946 DO i = 1, n2/2
1947 DO j = 1, nfft
1948 zw(1, j, i) = 0._dp
1949 zw(2, j, i) = 0._dp
1950 END DO
1951 END DO
1952 END SUBROUTINE switch_upcorn
1953
1954! **************************************************************************************************
1955!> \brief ...
1956!> \param j3 ...
1957!> \param nfft ...
1958!> \param Jp2stb ...
1959!> \param J2stb ...
1960!> \param lot ...
1961!> \param n1 ...
1962!> \param md2 ...
1963!> \param nd3 ...
1964!> \param nproc ...
1965!> \param zmpi1 ...
1966!> \param zw ...
1967! **************************************************************************************************
1968 SUBROUTINE mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot, n1, md2, nd3, nproc, zmpi1, zw)
1969 INTEGER :: j3, nfft, jp2stb, j2stb, lot, n1, md2, &
1970 nd3, nproc
1971 REAL(kind=dp) :: zmpi1(2, n1/2, md2/nproc, nd3/nproc, nproc), zw(2, lot, n1)
1972
1973 INTEGER :: i1, j2, jp2, mfft
1974
1975! WARNING: Assuming that high frequencies are in the corners
1976! and that n1 is multiple of 2
1977
1978 mfft = 0
1979 main: DO jp2 = jp2stb, nproc
1980 DO j2 = j2stb, md2/nproc
1981 mfft = mfft + 1
1982 IF (mfft > nfft) THEN
1983 jp2stb = jp2
1984 j2stb = j2
1985 EXIT main
1986 END IF
1987 DO i1 = 1, n1/2
1988 zw(1, mfft, i1) = 0._dp
1989 zw(2, mfft, i1) = 0._dp
1990 END DO
1991 DO i1 = n1/2 + 1, n1
1992 zw(1, mfft, i1) = zmpi1(1, i1 - n1/2, j2, j3, jp2)
1993 zw(2, mfft, i1) = zmpi1(2, i1 - n1/2, j2, j3, jp2)
1994 END DO
1995 END DO
1996 j2stb = 1
1997 END DO main
1998 END SUBROUTINE mpiswitch_upcorn
1999
2000! **************************************************************************************************
2001!> \brief ...
2002!> \param nfft ...
2003!> \param n2 ...
2004!> \param lot ...
2005!> \param n1 ...
2006!> \param lzt ...
2007!> \param zw ...
2008!> \param zt ...
2009! **************************************************************************************************
2010 SUBROUTINE unswitch_downcorn(nfft, n2, lot, n1, lzt, zw, zt)
2011 INTEGER :: nfft, n2, lot, n1, lzt
2012 REAL(kind=dp) :: zw(2, lot, n2), zt(2, lzt, n1)
2013
2014 INTEGER :: i, j
2015
2016! WARNING: Assuming that high frequencies are in the corners
2017! and that n2 is multiple of 2
2018! Low frequencies
2019
2020 DO j = 1, nfft
2021 DO i = 1, n2/2
2022 zt(1, i, j) = zw(1, j, i)
2023 zt(2, i, j) = zw(2, j, i)
2024 END DO
2025 END DO
2026 RETURN
2027 END SUBROUTINE unswitch_downcorn
2028
2029! **************************************************************************************************
2030!> \brief ...
2031!> \param j3 ...
2032!> \param nfft ...
2033!> \param Jp2stf ...
2034!> \param J2stf ...
2035!> \param lot ...
2036!> \param n1 ...
2037!> \param md2 ...
2038!> \param nd3 ...
2039!> \param nproc ...
2040!> \param zw ...
2041!> \param zmpi1 ...
2042! **************************************************************************************************
2043 SUBROUTINE unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot, n1, md2, nd3, nproc, zw, zmpi1)
2044 INTEGER :: j3, nfft, jp2stf, j2stf, lot, n1, md2, &
2045 nd3, nproc
2046 REAL(kind=dp) :: zw(2, lot, n1), zmpi1(2, n1/2, md2/nproc, nd3/nproc, nproc)
2047
2048 INTEGER :: i1, j2, jp2, mfft
2049
2050! WARNING: Assuming that high frequencies are in the corners
2051! and that n1 is multiple of 2
2052
2053 mfft = 0
2054 main: DO jp2 = jp2stf, nproc
2055 DO j2 = j2stf, md2/nproc
2056 mfft = mfft + 1
2057 IF (mfft > nfft) THEN
2058 jp2stf = jp2
2059 j2stf = j2
2060 EXIT main
2061 END IF
2062 DO i1 = 1, n1/2
2063 zmpi1(1, i1, j2, j3, jp2) = zw(1, mfft, i1)
2064 zmpi1(2, i1, j2, j3, jp2) = zw(2, mfft, i1)
2065 END DO
2066 END DO
2067 j2stf = 1
2068 END DO main
2069 END SUBROUTINE unmpiswitch_downcorn
2070
2071! **************************************************************************************************
2072!> \brief (Based on suitable modifications of S.Goedecker routines)
2073!> Restore data into output array, calculating in the meanwhile
2074!> Hartree energy of the potential
2075!> \param md1 Dimensions of the undistributed part of the real grid
2076!> \param md3 Dimensions of the undistributed part of the real grid
2077!> \param lot ...
2078!> \param nfft number of planes
2079!> \param n3 (twice the) dimension of the last FFTtransform.
2080!> \param zw FFT work array
2081!> \param zf Original distributed density as well as
2082!> Distributed solution of the poisson equation (inout)
2083!> \param scal Needed to achieve unitarity and correct dimensions
2084!> \param ehartreetmp Hartree energy
2085!> \date February 2006
2086!> \author S. Goedecker, L. Genovese
2087!> \note Assuming that high frequencies are in the corners
2088!> and that n3 is multiple of 4
2089!>
2090!> RESTRICTIONS on USAGE
2091!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
2092!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
2093!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
2094!> This file is distributed under the terms of the
2095!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
2096! **************************************************************************************************
2097 SUBROUTINE f_unfill_downcorn(md1, md3, lot, nfft, n3, zw, zf, scal, ehartreetmp)
2098 INTEGER, INTENT(in) :: md1, md3, lot, nfft, n3
2099 REAL(kind=dp), DIMENSION(2, lot, n3/2), INTENT(in) :: zw
2100 REAL(kind=dp), DIMENSION(md1, md3), INTENT(inout) :: zf
2101 REAL(kind=dp), INTENT(in) :: scal
2102 REAL(kind=dp), INTENT(out) :: ehartreetmp
2103
2104 INTEGER :: i1, i3
2105 REAL(kind=dp) :: pot1
2106
2107 ehartreetmp = 0._dp
2108 DO i3 = 1, n3/4
2109 DO i1 = 1, nfft
2110 pot1 = scal*zw(1, i1, i3)
2111 ehartreetmp = ehartreetmp + pot1*zf(i1, 2*i3 - 1)
2112 zf(i1, 2*i3 - 1) = pot1
2113 pot1 = scal*zw(2, i1, i3)
2114 ehartreetmp = ehartreetmp + pot1*zf(i1, 2*i3)
2115 zf(i1, 2*i3) = pot1
2116 END DO
2117 END DO
2118 END SUBROUTINE f_unfill_downcorn
2119
2120END MODULE ps_wavelet_base
main
int main(int argc, char *argv[])
Stand-alone miniapp for smoke-testing and benchmarking dbm_multiply.
Definition dbm_miniapp.c:248
for
for(int lxp=0;lxp<=lp;lxp++)
Definition grid_cpu_collint.h:66
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::sp
integer, parameter, public sp
Definition kinds.F:33
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
ps_wavelet_base
Creates the wavelet kernel for the wavelet based poisson solver.
Definition ps_wavelet_base.F:12
ps_wavelet_base::s_poissonsolver
subroutine, public s_poissonsolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, pot, zf, scal, mpi_group)
!HERE POT MUST BE THE KERNEL (BEWARE THE HALF DIMENSION) ****h* BigDFT/S_PoissonSolver (Based on suit...
Definition ps_wavelet_base.F:839
ps_wavelet_base::scramble_unpack
subroutine, public scramble_unpack(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zw, zmpi2, cosinarr)
(Based on suitable modifications of S.Goedecker routines) Assign the correct planes to the work array...
Definition ps_wavelet_base.F:1434
ps_wavelet_base::f_poissonsolver
subroutine, public f_poissonsolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, pot, zf, scal, mpi_group)
(Based on suitable modifications of S.Goedecker routines) Applies the local FFT space Kernel to the d...
Definition ps_wavelet_base.F:1583
ps_wavelet_base::p_poissonsolver
subroutine, public p_poissonsolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, zf, scal, hx, hy, hz, mpi_group)
...
Definition ps_wavelet_base.F:54
ps_wavelet_fft3d
Definition ps_wavelet_fft3d.F:9
ps_wavelet_fft3d::ctrig_length
integer, parameter, public ctrig_length
Definition ps_wavelet_fft3d.F:20
ps_wavelet_fft3d::fftstp
subroutine, public fftstp(mm, nfft, m, nn, n, zin, zout, trig, after, now, before, isign)
...
Definition ps_wavelet_fft3d.F:281
ps_wavelet_fft3d::ctrig
subroutine, public ctrig(n, trig, after, before, now, isign, ic)
...
Definition ps_wavelet_fft3d.F:106
message_passing::mp_comm_type
Definition message_passing.F:158