(git:b77b4be)
Loading...
Searching...
No Matches
pw_spline_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief different utils that are useful to manipulate splines on the regular
10!> grid of a pw
11!> \par History
12!> 05.2003 created [fawzi]
13!> 08.2004 removed spline evaluation method using more than 2 read streams
14!> (pw_compose_stripe_rs3), added linear solver based spline
15!> inversion [fawzi]
16!> \author Fawzi Mohamed
17! **************************************************************************************************
18MODULE pw_spline_utils
19
20 USE cp_log_handling, ONLY: cp_get_default_logger,&
21 cp_logger_type,&
22 cp_to_string
23 USE kinds, ONLY: dp
24 USE mathconstants, ONLY: twopi
25 USE message_passing, ONLY: mp_comm_congruent
26 USE pw_grid_types, ONLY: fullspace,&
27 pw_mode_local
28 USE pw_methods, ONLY: pw_axpy,&
29 pw_copy,&
30 pw_integral_ab,&
31 pw_zero
32 USE pw_pool_types, ONLY: pw_pool_release,&
33 pw_pool_type
34 USE pw_types, ONLY: pw_c1d_gs_type,&
35 pw_r3d_rs_type
36#include "../base/base_uses.f90"
37
38 IMPLICIT NONE
39 PRIVATE
40
41 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
42 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_spline_utils'
43
44 INTEGER, PARAMETER, PUBLIC :: no_precond = 0, &
45 precond_spl3_aint = 1, &
46 precond_spl3_1 = 2, &
47 precond_spl3_aint2 = 3, &
48 precond_spl3_2 = 4, &
49 precond_spl3_3 = 5
50
51 REAL(kind=dp), PUBLIC, PARAMETER, DIMENSION(4) :: nn10_coeffs = &
52 (/125._dp/216._dp, 25._dp/432._dp, 5._dp/864._dp, 1._dp/1728._dp/), &
53 spline3_coeffs = &
54 (/8._dp/(27._dp), 2._dp/(27._dp), 1._dp/(27._dp*2._dp), &
55 1._dp/(27._dp*8._dp)/), &
56 spline2_coeffs = &
57 (/27._dp/(64._dp), 9._dp/(64._dp*2_dp), 3._dp/(64._dp*4._dp), &
58 1._dp/(64._dp*8._dp)/), &
59 nn50_coeffs = &
60 (/15625._dp/17576._dp, 625._dp/35152._dp, 25._dp/70304._dp, &
61 1._dp/140608._dp/), &
62 spl3_aint_coeff = &
63 (/46._dp/27._dp, -2._dp/(27._dp), -1._dp/(27._dp*2._dp), &
64 -1._dp/(27._dp*8._dp)/), &
65 spl3_precond1_coeff = &
66 (/64._dp/3._dp, -8._dp/3._dp, -1._dp/3._dp, -1._dp/24._dp/), &
67 spl3_1d_transf_coeffs = &
68 (/2._dp/3._dp, 23._dp/48._dp, 1._dp/6._dp, 1._dp/48._dp/)
69
70 REAL(kind=dp), PUBLIC, PARAMETER, DIMENSION(3) :: spline3_deriv_coeffs = &
71 (/2.0_dp/9.0_dp, 1.0_dp/18.0_dp, 1.0_dp/72.0_dp/), &
72 spline2_deriv_coeffs = &
73 (/9.0_dp/32.0_dp, 3.0_dp/64.0_dp, 1.0_dp/128.0_dp/), &
74 nn10_deriv_coeffs = &
75 (/25._dp/72._dp, 5._dp/144, 1._dp/288._dp/), &
76 nn50_deriv_coeffs = &
77 (/625._dp/1352._dp, 25._dp/2704._dp, 1._dp/5408._dp/), &
78 spl3_1d_coeffs0 = &
79 (/1._dp/6_dp, 2._dp/3._dp, 1._dp/6._dp/), &
80 spl3_1d_transf_border1 = &
81 (/0.517977704_dp, 0.464044595_dp, 0.17977701e-1_dp/)
82
83 PUBLIC :: pw_spline3_interpolate_values_g, &
84 pw_spline3_deriv_g
85 PUBLIC :: pw_spline_scale_deriv
86 PUBLIC :: pw_spline2_interpolate_values_g, &
87 pw_spline2_deriv_g
88 PUBLIC :: pw_nn_smear_r, pw_nn_deriv_r, &
89 spl3_nopbc, spl3_pbc, spl3_nopbct
90 PUBLIC :: add_fine2coarse, add_coarse2fine
91 PUBLIC :: pw_spline_precond_create, &
92 pw_spline_do_precond, &
93 pw_spline_precond_set_kind, &
94 find_coeffs, &
95 pw_spline_precond_release, &
96 pw_spline_precond_type, &
97 eval_interp_spl3_pbc, &
98 eval_d_interp_spl3_pbc
99
100!***
101
102! **************************************************************************************************
103!> \brief stores information for the preconditioner used to calculate the
104!> coeffs of splines
105!> \author fawzi
106! **************************************************************************************************
107 TYPE pw_spline_precond_type
108 INTEGER :: kind = no_precond
109 REAL(kind=dp), DIMENSION(4) :: coeffs = 0.0_dp
110 REAL(kind=dp), DIMENSION(3) :: coeffs_1d = 0.0_dp
111 LOGICAL :: sharpen = .false., normalize = .false., pbc = .false., transpose = .false.
112 TYPE(pw_pool_type), POINTER :: pool => null()
113 END TYPE pw_spline_precond_type
114
115CONTAINS
116
117! **************************************************************************************************
118!> \brief calculates the FFT of the coefficients of the quadratic spline that
119!> interpolates the given values
120!> \param spline_g on entry the FFT of the values to interpolate as cc,
121!> will contain the FFT of the coefficients of the spline
122!> \par History
123!> 06.2003 created [fawzi]
124!> \author Fawzi Mohamed
125!> \note
126!> does not work with spherical cutoff
127! **************************************************************************************************
128 SUBROUTINE pw_spline2_interpolate_values_g(spline_g)
129 TYPE(pw_c1d_gs_type), INTENT(IN) :: spline_g
130
131 CHARACTER(len=*), PARAMETER :: routinen = 'pw_spline2_interpolate_values_g'
132
133 INTEGER :: handle, i, ii, j, k
134 INTEGER, DIMENSION(2, 3) :: gbo
135 INTEGER, DIMENSION(3) :: n_tot
136 REAL(kind=dp) :: c23, coeff
137 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cosivals, cosjvals, coskvals
138
139 CALL timeset(routinen, handle)
140
141 n_tot(1:3) = spline_g%pw_grid%npts(1:3)
142 gbo = spline_g%pw_grid%bounds
143
144 cpassert(.NOT. spline_g%pw_grid%spherical)
145 cpassert(spline_g%pw_grid%grid_span == fullspace)
146
147 ALLOCATE (cosivals(gbo(1, 1):gbo(2, 1)), cosjvals(gbo(1, 2):gbo(2, 2)), &
148 coskvals(gbo(1, 3):gbo(2, 3)))
149
150 coeff = twopi/n_tot(1)
151!$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(cosIVals,coeff,gbo)
152 DO i = gbo(1, 1), gbo(2, 1)
153 cosivals(i) = cos(coeff*real(i, dp))
154 END DO
155 coeff = twopi/n_tot(2)
156!$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(cosJVals,coeff,gbo)
157 DO j = gbo(1, 2), gbo(2, 2)
158 cosjvals(j) = cos(coeff*real(j, dp))
159 END DO
160 coeff = twopi/n_tot(3)
161!$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(cosKVals,coeff,gbo)
162 DO k = gbo(1, 3), gbo(2, 3)
163 coskvals(k) = cos(coeff*real(k, dp))
164 END DO
165
166!$OMP PARALLEL DO PRIVATE(i,j,k,ii,coeff,c23) DEFAULT(NONE) SHARED(spline_g,cosIVals,cosJVals,cosKVals)
167 DO ii = 1, SIZE(spline_g%array)
168 i = spline_g%pw_grid%g_hat(1, ii)
169 j = spline_g%pw_grid%g_hat(2, ii)
170 k = spline_g%pw_grid%g_hat(3, ii)
171
172 c23 = cosjvals(j)*coskvals(k)
173 coeff = 64.0_dp/(cosivals(i)*c23 + &
174 (cosivals(i)*cosjvals(j) + cosivals(i)*coskvals(k) + c23)*3.0_dp + &
175 (cosivals(i) + cosjvals(j) + coskvals(k))*9.0_dp + &
176 27.0_dp)
177
178 spline_g%array(ii) = spline_g%array(ii)*coeff
179
180 END DO
181 DEALLOCATE (cosivals, cosjvals, coskvals)
182
183 CALL timestop(handle)
184 END SUBROUTINE pw_spline2_interpolate_values_g
185
186! **************************************************************************************************
187!> \brief calculates the FFT of the coefficients of the2 cubic spline that
188!> interpolates the given values
189!> \param spline_g on entry the FFT of the values to interpolate as cc,
190!> will contain the FFT of the coefficients of the spline
191!> \par History
192!> 06.2003 created [fawzi]
193!> \author Fawzi Mohamed
194!> \note
195!> does not work with spherical cutoff
196!> stupid distribution for cos calculation, it should calculate only the
197!> needed cos, and avoid the mpi_allreduce
198! **************************************************************************************************
199 SUBROUTINE pw_spline3_interpolate_values_g(spline_g)
200 TYPE(pw_c1d_gs_type), INTENT(IN) :: spline_g
201
202 CHARACTER(len=*), PARAMETER :: routinen = 'pw_spline3_interpolate_values_g'
203
204 INTEGER :: handle, i, ii, j, k
205 INTEGER, DIMENSION(2, 3) :: gbo
206 INTEGER, DIMENSION(3) :: n_tot
207 REAL(kind=dp) :: c23, coeff
208 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cosivals, cosjvals, coskvals
209
210 CALL timeset(routinen, handle)
211
212 n_tot(1:3) = spline_g%pw_grid%npts(1:3)
213 gbo = spline_g%pw_grid%bounds
214
215 cpassert(.NOT. spline_g%pw_grid%spherical)
216 cpassert(spline_g%pw_grid%grid_span == fullspace)
217
218 ALLOCATE (cosivals(gbo(1, 1):gbo(2, 1)), &
219 cosjvals(gbo(1, 2):gbo(2, 2)), &
220 coskvals(gbo(1, 3):gbo(2, 3)))
221
222 coeff = twopi/n_tot(1)
223!$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(cosIVals,coeff,gbo)
224 DO i = gbo(1, 1), gbo(2, 1)
225 cosivals(i) = cos(coeff*real(i, dp))
226 END DO
227 coeff = twopi/n_tot(2)
228!$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(cosJVals,coeff,gbo)
229 DO j = gbo(1, 2), gbo(2, 2)
230 cosjvals(j) = cos(coeff*real(j, dp))
231 END DO
232 coeff = twopi/n_tot(3)
233!$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(cosKVals,coeff,gbo)
234 DO k = gbo(1, 3), gbo(2, 3)
235 coskvals(k) = cos(coeff*real(k, dp))
236 END DO
237
238!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k,ii,coeff,c23) SHARED(spline_g,cosIVals,cosJVals,cosKVals)
239 DO ii = 1, SIZE(spline_g%array)
240 i = spline_g%pw_grid%g_hat(1, ii)
241 j = spline_g%pw_grid%g_hat(2, ii)
242 k = spline_g%pw_grid%g_hat(3, ii)
243 ! no opt
244!FM coeff=1.0/((cosVal(1)*cosVal(2)*cosVal(3))/27.0_dp+&
245!FM (cosVal(1)*cosVal(2)+cosVal(1)*cosVal(3)+&
246!FM cosVal(2)*cosVal(3))*2.0_dp/27.0_dp+&
247!FM (cosVal(1)+cosVal(2)+cosVal(3))*4.0_dp/27.0_dp+&
248!FM 8.0_dp/27.0_dp)
249 ! opt
250 c23 = cosjvals(j)*coskvals(k)
251 coeff = 27.0_dp/(cosivals(i)*c23 + &
252 (cosivals(i)*cosjvals(j) + cosivals(i)*coskvals(k) + c23)*2.0_dp + &
253 (cosivals(i) + cosjvals(j) + coskvals(k))*4.0_dp + &
254 8.0_dp)
255
256 spline_g%array(ii) = spline_g%array(ii)*coeff
257
258 END DO
259 DEALLOCATE (cosivals, cosjvals, coskvals)
260
261 CALL timestop(handle)
262 END SUBROUTINE pw_spline3_interpolate_values_g
263
264! **************************************************************************************************
265!> \brief rescales the derivatives from gridspacing=1 to the real derivatives
266!> \param deriv_vals_r an array of x,y,z derivatives
267!> \param transpose if true applies the transpose of the map (defaults to
268!> false)
269!> \param scale a scaling factor (defaults to 1.0)
270!> \par History
271!> 06.2003 created [fawzi]
272!> \author Fawzi Mohamed
273! **************************************************************************************************
274 SUBROUTINE pw_spline_scale_deriv(deriv_vals_r, transpose, scale)
275 TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(IN) :: deriv_vals_r
276 LOGICAL, INTENT(in), OPTIONAL :: transpose
277 REAL(kind=dp), INTENT(in), OPTIONAL :: scale
278
279 CHARACTER(len=*), PARAMETER :: routinen = 'pw_spline_scale_deriv'
280
281 INTEGER :: handle, i, idir, j, k
282 INTEGER, DIMENSION(2, 3) :: bo
283 INTEGER, DIMENSION(3) :: n_tot
284 LOGICAL :: diag, my_transpose
285 REAL(kind=dp) :: dval1, dval2, dval3, my_scale, scalef
286 REAL(kind=dp), DIMENSION(3, 3) :: dh_inv, h_grid
287 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ddata, ddata2, ddata3
288
289 CALL timeset(routinen, handle)
290
291 my_transpose = .false.
292 IF (PRESENT(transpose)) my_transpose = transpose
293 my_scale = 1.0_dp
294 IF (PRESENT(scale)) my_scale = scale
295 n_tot(1:3) = deriv_vals_r(1)%pw_grid%npts(1:3)
296 bo = deriv_vals_r(1)%pw_grid%bounds_local
297 dh_inv = deriv_vals_r(1)%pw_grid%dh_inv
298
299 ! map grid to real derivative
300 diag = .true.
301 IF (my_transpose) THEN
302 DO j = 1, 3
303 DO i = 1, 3
304 h_grid(j, i) = my_scale*dh_inv(i, j) ! REAL(n_tot(i),dp)*cell_h_inv(i,j)
305 IF (i /= j .AND. h_grid(j, i) /= 0.0_dp) diag = .false.
306 END DO
307 END DO
308 ELSE
309 DO j = 1, 3
310 DO i = 1, 3
311 h_grid(i, j) = my_scale*dh_inv(i, j) ! REAL(n_tot(i),dp)*cell_h_inv(i,j)
312 IF (i /= j .AND. h_grid(i, j) /= 0.0_dp) diag = .false.
313 END DO
314 END DO
315 END IF
316
317 IF (diag) THEN
318 DO idir = 1, 3
319 ddata => deriv_vals_r(idir)%array
320 scalef = h_grid(idir, idir)
321 CALL dscal((bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1), &
322 scalef, ddata, 1)
323 END DO
324 ELSE
325 ddata => deriv_vals_r(1)%array
326 ddata2 => deriv_vals_r(2)%array
327 ddata3 => deriv_vals_r(3)%array
328!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(k,j,i,dVal1,dVal2,dVal3) &
329!$OMP SHARED(ddata,ddata2,ddata3,h_grid,bo)
330 DO k = bo(1, 3), bo(2, 3)
331 DO j = bo(1, 2), bo(2, 2)
332 DO i = bo(1, 1), bo(2, 1)
333
334 dval1 = ddata(i, j, k)
335 dval2 = ddata2(i, j, k)
336 dval3 = ddata3(i, j, k)
337
338 ddata(i, j, k) = h_grid(1, 1)*dval1 + &
339 h_grid(2, 1)*dval2 + h_grid(3, 1)*dval3
340 ddata2(i, j, k) = h_grid(1, 2)*dval1 + &
341 h_grid(2, 2)*dval2 + h_grid(3, 2)*dval3
342 ddata3(i, j, k) = h_grid(1, 3)*dval1 + &
343 h_grid(2, 3)*dval2 + h_grid(3, 3)*dval3
344
345 END DO
346 END DO
347 END DO
348 END IF
349
350 CALL timestop(handle)
351 END SUBROUTINE pw_spline_scale_deriv
352
353! **************************************************************************************************
354!> \brief calculates the FFT of the values of the x,y,z (idir=1,2,3)
355!> derivative of the cubic spline
356!> \param spline_g on entry the FFT of the coefficients of the spline
357!> will contain the FFT of the derivative
358!> \param idir direction of the derivative
359!> \par History
360!> 06.2003 created [fawzi]
361!> \author Fawzi Mohamed
362!> \note
363!> the distance between gridpoints is assumed to be 1
364! **************************************************************************************************
365 SUBROUTINE pw_spline3_deriv_g(spline_g, idir)
366 TYPE(pw_c1d_gs_type), INTENT(IN) :: spline_g
367 INTEGER, INTENT(in) :: idir
368
369 CHARACTER(len=*), PARAMETER :: routinen = 'pw_spline3_deriv_g'
370 REAL(kind=dp), PARAMETER :: inv9 = 1.0_dp/9.0_dp
371
372 INTEGER :: handle, i, ii, j, k
373 INTEGER, DIMENSION(2, 3) :: bo, gbo
374 INTEGER, DIMENSION(3) :: n, n_tot
375 REAL(kind=dp) :: coeff, tmp
376 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: csivals, csjvals, cskvals
377
378 CALL timeset(routinen, handle)
379
380 n(1:3) = spline_g%pw_grid%npts_local(1:3)
381 n_tot(1:3) = spline_g%pw_grid%npts(1:3)
382 bo = spline_g%pw_grid%bounds_local
383 gbo = spline_g%pw_grid%bounds
384
385 cpassert(.NOT. spline_g%pw_grid%spherical)
386 cpassert(spline_g%pw_grid%grid_span == fullspace)
387
388 ALLOCATE (csivals(gbo(1, 1):gbo(2, 1)), &
389 csjvals(gbo(1, 2):gbo(2, 2)), &
390 cskvals(gbo(1, 3):gbo(2, 3)))
391
392 coeff = twopi/n_tot(1)
393 IF (idir == 1) THEN
394!$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(gbo,csIVals,coeff)
395 DO i = gbo(1, 1), gbo(2, 1)
396 csivals(i) = sin(coeff*real(i, dp))
397 END DO
398 ELSE
399!$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(gbo,csIVals,coeff)
400 DO i = gbo(1, 1), gbo(2, 1)
401 csivals(i) = cos(coeff*real(i, dp))
402 END DO
403 END IF
404 coeff = twopi/n_tot(2)
405 IF (idir == 2) THEN
406!$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(gbo,csJVals,coeff)
407 DO j = gbo(1, 2), gbo(2, 2)
408 csjvals(j) = sin(coeff*real(j, dp))
409 END DO
410 ELSE
411!$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(gbo,csJVals,coeff)
412 DO j = gbo(1, 2), gbo(2, 2)
413 csjvals(j) = cos(coeff*real(j, dp))
414 END DO
415 END IF
416 coeff = twopi/n_tot(3)
417 IF (idir == 3) THEN
418!$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(gbo,csKVals,coeff)
419 DO k = gbo(1, 3), gbo(2, 3)
420 cskvals(k) = sin(coeff*real(k, dp))
421 END DO
422 ELSE
423!$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(gbo,csKVals,coeff)
424 DO k = gbo(1, 3), gbo(2, 3)
425 cskvals(k) = cos(coeff*real(k, dp))
426 END DO
427 END IF
428
429 SELECT CASE (idir)
430 CASE (1)
431 ! x deriv
432!$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) DEFAULT(NONE) SHARED(spline_g,csIVals,csJVals,csKVals)
433 DO ii = 1, SIZE(spline_g%array)
434 i = spline_g%pw_grid%g_hat(1, ii)
435 j = spline_g%pw_grid%g_hat(2, ii)
436 k = spline_g%pw_grid%g_hat(3, ii)
437!FM ! formula
438!FM coeff=(sinVal(1)*cosVal(2)*cosVal(3))/9.0_dp+&
439!FM (sinVal(1)*cosVal(2)+sinVal(1)*cosVal(3))*2.0_dp/9.0_dp+&
440!FM sinVal(1)*4.0_dp/9.0_dp
441 tmp = csivals(i)*csjvals(j)
442 coeff = (tmp*cskvals(k) + &
443 (tmp + csivals(i)*cskvals(k))*2.0_dp + &
444 csivals(i)*4.0_dp)*inv9
445
446 spline_g%array(ii) = spline_g%array(ii)* &
447 cmplx(0.0_dp, coeff, dp)
448 END DO
449 CASE (2)
450 ! y deriv
451!$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) DEFAULT(NONE) SHARED(spline_g,csIVals,csJVals,csKVals)
452 DO ii = 1, SIZE(spline_g%array)
453 i = spline_g%pw_grid%g_hat(1, ii)
454 j = spline_g%pw_grid%g_hat(2, ii)
455 k = spline_g%pw_grid%g_hat(3, ii)
456
457 tmp = csivals(i)*csjvals(j)
458 coeff = (tmp*cskvals(k) + &
459 (tmp + csjvals(j)*cskvals(k))*2.0_dp + &
460 csjvals(j)*4.0_dp)*inv9
461
462 spline_g%array(ii) = spline_g%array(ii)* &
463 cmplx(0.0_dp, coeff, dp)
464 END DO
465 CASE (3)
466 ! z deriv
467!$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) DEFAULT(NONE) SHARED(spline_g,csIVals,csJVals,csKVals)
468 DO ii = 1, SIZE(spline_g%array)
469 i = spline_g%pw_grid%g_hat(1, ii)
470 j = spline_g%pw_grid%g_hat(2, ii)
471 k = spline_g%pw_grid%g_hat(3, ii)
472
473 tmp = csivals(i)*cskvals(k)
474 coeff = (tmp*csjvals(j) + &
475 (tmp + csjvals(j)*cskvals(k))*2.0_dp + &
476 cskvals(k)*4.0_dp)*inv9
477
478 spline_g%array(ii) = spline_g%array(ii)* &
479 cmplx(0.0_dp, coeff, dp)
480 END DO
481 END SELECT
482
483 DEALLOCATE (csivals, csjvals, cskvals)
484
485 CALL timestop(handle)
486 END SUBROUTINE pw_spline3_deriv_g
487
488! **************************************************************************************************
489!> \brief calculates the FFT of the values of the x,y,z (idir=1,2,3)
490!> derivative of the quadratic spline
491!> \param spline_g on entry the FFT of the coefficients of the spline
492!> will contain the FFT of the derivative
493!> \param idir direction of the derivative
494!> \par History
495!> 06.2003 created [fawzi]
496!> \author Fawzi Mohamed
497!> \note
498!> the distance between gridpoints is assumed to be 1
499! **************************************************************************************************
500 SUBROUTINE pw_spline2_deriv_g(spline_g, idir)
501 TYPE(pw_c1d_gs_type), INTENT(IN) :: spline_g
502 INTEGER, INTENT(in) :: idir
503
504 CHARACTER(len=*), PARAMETER :: routinen = 'pw_spline2_deriv_g'
505 REAL(kind=dp), PARAMETER :: inv16 = 1.0_dp/16.0_dp
506
507 INTEGER :: handle, i, ii, j, k
508 INTEGER, DIMENSION(2, 3) :: bo
509 INTEGER, DIMENSION(3) :: n, n_tot
510 REAL(kind=dp) :: coeff, tmp
511 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: csivals, csjvals, cskvals
512
513 CALL timeset(routinen, handle)
514
515 n(1:3) = spline_g%pw_grid%npts_local(1:3)
516 n_tot(1:3) = spline_g%pw_grid%npts(1:3)
517 bo = spline_g%pw_grid%bounds
518
519 cpassert(.NOT. spline_g%pw_grid%spherical)
520 cpassert(spline_g%pw_grid%grid_span == fullspace)
521
522 ALLOCATE (csivals(bo(1, 1):bo(2, 1)), csjvals(bo(1, 2):bo(2, 2)), &
523 cskvals(bo(1, 3):bo(2, 3)))
524
525 coeff = twopi/n_tot(1)
526 IF (idir == 1) THEN
527!$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(bo,coeff,csIVals)
528 DO i = bo(1, 1), bo(2, 1)
529 csivals(i) = sin(coeff*real(i, dp))
530 END DO
531 ELSE
532!$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(bo,coeff,csIVals)
533 DO i = bo(1, 1), bo(2, 1)
534 csivals(i) = cos(coeff*real(i, dp))
535 END DO
536 END IF
537 coeff = twopi/n_tot(2)
538 IF (idir == 2) THEN
539!$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(bo,coeff,csJVals)
540 DO j = bo(1, 2), bo(2, 2)
541 csjvals(j) = sin(coeff*real(j, dp))
542 END DO
543 ELSE
544!$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(bo,coeff,csJVals)
545 DO j = bo(1, 2), bo(2, 2)
546 csjvals(j) = cos(coeff*real(j, dp))
547 END DO
548 END IF
549 coeff = twopi/n_tot(3)
550 IF (idir == 3) THEN
551!$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(bo,coeff,csKVals)
552 DO k = bo(1, 3), bo(2, 3)
553 cskvals(k) = sin(coeff*real(k, dp))
554 END DO
555 ELSE
556!$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(bo,coeff,csKVals)
557 DO k = bo(1, 3), bo(2, 3)
558 cskvals(k) = cos(coeff*real(k, dp))
559 END DO
560 END IF
561
562 SELECT CASE (idir)
563 CASE (1)
564 ! x deriv
565!$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) SHARED(spline_g,csIVals,csJVals,csKVals) DEFAULT(NONE)
566 DO ii = 1, SIZE(spline_g%array)
567 i = spline_g%pw_grid%g_hat(1, ii)
568 j = spline_g%pw_grid%g_hat(2, ii)
569 k = spline_g%pw_grid%g_hat(3, ii)
570!FM ! formula
571!FM coeff=(sinVal(1)*cosVal(2)*cosVal(3))/16.0_dp+&
572!FM (sinVal(1)*cosVal(2)+sinVal(1)*cosVal(3))*3.0_dp/16.0_dp+&
573!FM sinVal(1)*9.0_dp/16.0_dp
574 tmp = csivals(i)*csjvals(j)
575 coeff = (tmp*cskvals(k) + &
576 (tmp + csivals(i)*cskvals(k))*3.0_dp + &
577 csivals(i)*9.0_dp)*inv16
578
579 spline_g%array(ii) = spline_g%array(ii)* &
580 cmplx(0.0_dp, coeff, dp)
581 END DO
582 CASE (2)
583 ! y deriv
584!$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) DEFAULT(NONE) SHARED(spline_g,csIVals,csJVals,csKVals)
585 DO ii = 1, SIZE(spline_g%array)
586 i = spline_g%pw_grid%g_hat(1, ii)
587 j = spline_g%pw_grid%g_hat(2, ii)
588 k = spline_g%pw_grid%g_hat(3, ii)
589
590 tmp = csivals(i)*csjvals(j)
591 coeff = (tmp*cskvals(k) + &
592 (tmp + csjvals(j)*cskvals(k))*3.0_dp + &
593 csjvals(j)*9.0_dp)*inv16
594
595 spline_g%array(ii) = spline_g%array(ii)* &
596 cmplx(0.0_dp, coeff, dp)
597 END DO
598 CASE (3)
599 ! z deriv
600!$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) DEFAULT(NONE) SHARED(spline_g,csIVals,csJVals,csKVals)
601 DO ii = 1, SIZE(spline_g%array)
602 i = spline_g%pw_grid%g_hat(1, ii)
603 j = spline_g%pw_grid%g_hat(2, ii)
604 k = spline_g%pw_grid%g_hat(3, ii)
605
606 tmp = csivals(i)*cskvals(k)
607 coeff = (tmp*csjvals(j) + &
608 (tmp + csjvals(j)*cskvals(k))*3.0_dp + &
609 cskvals(k)*9.0_dp)*inv16
610
611 spline_g%array(ii) = spline_g%array(ii)* &
612 cmplx(0.0_dp, coeff, dp)
613 END DO
614 END SELECT
615
616 DEALLOCATE (csivals, csjvals, cskvals)
617
618 CALL timestop(handle)
619 END SUBROUTINE pw_spline2_deriv_g
620
621! **************************************************************************************************
622!> \brief applies a nearest neighbor linear operator to a stripe in x direction:
623!> out_val(i)=sum(weight(j)*in_val(i+j-1),j=0..2)
624!> \param weights the weights of the linear operator
625!> \param in_val the argument to the operator
626!> \param in_val_first the first argument (needed to calculate out_val(1))
627!> \param in_val_last the last argument (needed to calculate out_val(n_el))
628!> \param out_val the place where the result is accumulated
629!> \param n_el the number of elements in in_v and out_v
630!> \par History
631!> 04.2004 created [fawzi]
632!> \author fawzi
633!> \note
634!> uses 2 read streams and 1 write stream
635! **************************************************************************************************
636 SUBROUTINE pw_compose_stripe(weights, in_val, in_val_first, in_val_last, &
637 out_val, n_el)
638 REAL(kind=dp), DIMENSION(0:2), INTENT(in) :: weights
639 REAL(kind=dp), DIMENSION(*), INTENT(in) :: in_val
640 REAL(kind=dp), INTENT(in) :: in_val_first, in_val_last
641 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: out_val
642 INTEGER :: n_el
643
644 INTEGER :: i
645 REAL(kind=dp) :: v0, v1, v2
646
647!1:n_el), &
648!1:n_el), &
649
650 IF (n_el < 1) RETURN
651 v0 = in_val_first
652 v1 = in_val(1)
653 IF (weights(1) == 0.0_dp) THEN
654 ! optimized version for x deriv
655 DO i = 1, n_el - 3, 3
656 v2 = in_val(i + 1)
657 out_val(i) = out_val(i) + &
658 weights(0)*v0 + &
659 weights(2)*v2
660 v0 = in_val(i + 2)
661 out_val(i + 1) = out_val(i + 1) + &
662 weights(0)*v1 + &
663 weights(2)*v0
664 v1 = in_val(i + 3)
665 out_val(i + 2) = out_val(i + 2) + &
666 weights(0)*v2 + &
667 weights(2)*v1
668 END DO
669 ELSE
670 ! generic version
671 DO i = 1, n_el - 3, 3
672 v2 = in_val(i + 1)
673 out_val(i) = out_val(i) + &
674 weights(0)*v0 + &
675 weights(1)*v1 + &
676 weights(2)*v2
677 v0 = in_val(i + 2)
678 out_val(i + 1) = out_val(i + 1) + &
679 weights(0)*v1 + &
680 weights(1)*v2 + &
681 weights(2)*v0
682 v1 = in_val(i + 3)
683 out_val(i + 2) = out_val(i + 2) + &
684 weights(0)*v2 + &
685 weights(1)*v0 + &
686 weights(2)*v1
687 END DO
688 END IF
689 SELECT CASE (modulo(n_el - 1, 3))
690 CASE (0)
691 v2 = in_val_last
692 out_val(n_el) = out_val(n_el) + &
693 weights(0)*v0 + &
694 weights(1)*v1 + &
695 weights(2)*v2
696 CASE (1)
697 v2 = in_val(n_el)
698 out_val(n_el - 1) = out_val(n_el - 1) + &
699 weights(0)*v0 + &
700 weights(1)*v1 + &
701 weights(2)*v2
702 v0 = in_val_last
703 out_val(n_el) = out_val(n_el) + &
704 weights(0)*v1 + &
705 weights(1)*v2 + &
706 weights(2)*v0
707 CASE (2)
708 v2 = in_val(n_el - 1)
709 out_val(n_el - 2) = out_val(n_el - 2) + &
710 weights(0)*v0 + &
711 weights(1)*v1 + &
712 weights(2)*v2
713 v0 = in_val(n_el)
714 out_val(n_el - 1) = out_val(n_el - 1) + &
715 weights(0)*v1 + &
716 weights(1)*v2 + &
717 weights(2)*v0
718 v1 = in_val_last
719 out_val(n_el) = out_val(n_el) + &
720 weights(0)*v2 + &
721 weights(1)*v0 + &
722 weights(2)*v1
723 END SELECT
724
725 END SUBROUTINE pw_compose_stripe
726
727! **************************************************************************************************
728!> \brief private routine that computes pw_nn_compose_r (it seems that without
729!> passing arrays in this way either some compiler do a copyin/out (xlf)
730!> or by inlining suboptimal code is produced (nag))
731!> \param weights a 3x3x3 array with the linear operator
732!> \param in_val the argument for the linear operator
733!> \param out_val place where the value of the linear oprator should be added
734!> \param pw_in pw to be able to get the needed meta data about in_val and
735!> out_val
736!> \param bo boundaries of in_val and out_val
737!> \author fawzi
738! **************************************************************************************************
739 SUBROUTINE pw_nn_compose_r_work(weights, in_val, out_val, pw_in, bo)
740 REAL(kind=dp), DIMENSION(0:2, 0:2, 0:2) :: weights
741 INTEGER, DIMENSION(2, 3) :: bo
742 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
743 REAL(kind=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, & 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(inout) :: out_val
744 REAL(kind=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, & 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(in) :: in_val
745
746 INTEGER :: i, j, jw, k, kw, myj, myk
747 INTEGER, DIMENSION(2, 3) :: gbo
748 INTEGER, DIMENSION(3) :: s
749 LOGICAL :: has_boundary, yderiv, zderiv
750 REAL(kind=dp) :: in_val_f, in_val_l
751 REAL(kind=dp), DIMENSION(:, :), POINTER :: l_boundary, tmp, u_boundary
752
753 zderiv = all(weights(:, :, 1) == 0.0_dp)
754 yderiv = all(weights(:, 1, :) == 0.0_dp)
755 bo = pw_in%pw_grid%bounds_local
756 gbo = pw_in%pw_grid%bounds
757 DO i = 1, 3
758 s(i) = bo(2, i) - bo(1, i) + 1
759 END DO
760 IF (any(s < 1)) RETURN
761 has_boundary = any(pw_in%pw_grid%bounds_local(:, 1) /= &
762 pw_in%pw_grid%bounds(:, 1))
763 IF (has_boundary) THEN
764 ALLOCATE (l_boundary(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)), &
765 u_boundary(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)), &
766 tmp(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
767 tmp(:, :) = pw_in%array(bo(2, 1), :, :)
768 CALL pw_in%pw_grid%para%group%sendrecv(tmp, pw_in%pw_grid%para%pos_of_x( &
769 gbo(1, 1) + modulo(bo(2, 1) + 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)), &
770 l_boundary, pw_in%pw_grid%para%pos_of_x( &
771 gbo(1, 1) + modulo(bo(1, 1) - 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)))
772 tmp(:, :) = pw_in%array(bo(1, 1), :, :)
773 CALL pw_in%pw_grid%para%group%sendrecv(tmp, pw_in%pw_grid%para%pos_of_x( &
774 gbo(1, 1) + modulo(bo(1, 1) - 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)), &
775 u_boundary, pw_in%pw_grid%para%pos_of_x( &
776 gbo(1, 1) + modulo(bo(2, 1) + 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)))
777 DEALLOCATE (tmp)
778 END IF
779
780!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(k,kw,myk,j,jw,myj,in_val_f,&
781!$OMP in_val_l) SHARED(zderiv,yderiv,bo,in_val,out_val,s,l_boundary,&
782!$OMP u_boundary,weights,has_boundary)
783 DO k = 0, s(3) - 1
784 DO kw = 0, 2
785 myk = bo(1, 3) + modulo(k + kw - 1, s(3))
786 IF (zderiv .AND. kw == 1) cycle
787 DO j = 0, s(2) - 1
788 DO jw = 0, 2
789 myj = bo(1, 2) + modulo(j + jw - 1, s(2))
790 IF (yderiv .AND. jw == 1) cycle
791 IF (has_boundary) THEN
792 in_val_f = l_boundary(myj, myk)
793 in_val_l = u_boundary(myj, myk)
794 ELSE
795 in_val_f = in_val(bo(2, 1), myj, myk)
796 in_val_l = in_val(bo(1, 1), myj, myk)
797 END IF
798 CALL pw_compose_stripe(weights=weights(:, jw, kw), &
799 in_val=in_val(:, myj, myk), &
800 in_val_first=in_val_f, in_val_last=in_val_l, &
801 out_val=out_val(:, bo(1, 2) + j, bo(1, 3) + k), n_el=s(1))
802 END DO
803 END DO
804 END DO
805 END DO
806 IF (has_boundary) THEN
807 DEALLOCATE (l_boundary, u_boundary)
808 END IF
809 END SUBROUTINE pw_nn_compose_r_work
810
811! **************************************************************************************************
812!> \brief applies a nearest neighbor linear operator to a pw in real space
813!> \param weights a 3x3x3 array with the linear operator
814!> \param pw_in the argument for the linear operator
815!> \param pw_out place where the value of the linear oprator should be added
816!> \author fawzi
817!> \note
818!> has specialized versions for derivative operator (with central values==0)
819! **************************************************************************************************
820 SUBROUTINE pw_nn_compose_r(weights, pw_in, pw_out)
821 REAL(kind=dp), DIMENSION(0:2, 0:2, 0:2) :: weights
822 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in, pw_out
823
824 CHARACTER(len=*), PARAMETER :: routinen = 'pw_nn_compose_r'
825
826 INTEGER :: handle
827
828 CALL timeset(routinen, handle)
829 IF (.NOT. all(pw_in%pw_grid%bounds_local(:, 2:3) == pw_in%pw_grid%bounds(:, 2:3))) THEN
830 cpabort("wrong pw distribution")
831 END IF
832 CALL pw_nn_compose_r_work(weights=weights, in_val=pw_in%array, &
833 out_val=pw_out%array, pw_in=pw_in, bo=pw_in%pw_grid%bounds_local)
834 CALL timestop(handle)
835 END SUBROUTINE pw_nn_compose_r
836
837! **************************************************************************************************
838!> \brief calculates the values of a nearest neighbor smearing
839!> \param pw_in the argument for the linear operator
840!> \param pw_out place where the smeared values should be added
841!> \param coeffs array with the coefficent of the smearing, ordered with
842!> the distance from the center: coeffs(1) the coeff of the central
843!> element, coeffs(2) the coeff of the 6 element with distance 1,
844!> coeff(3) the coeff of the 12 elements at distance sqrt(2),
845!> coeff(4) the coeff of the 8 elements at distance sqrt(3).
846!> \author Fawzi Mohamed
847!> \note
848!> does not normalize the smear to 1.
849!> with coeff=(/ 8._dp/27._dp, 2._dp/27._dp, 1._dp/54._dp, 1._dp/216._dp /)
850!> is equivalent to pw_spline3_evaluate_values_g, with
851!> coeff=(/ 27._dp/64._dp, 9._dp/128._dp, 3._dp/256._dp, 1._dp/512._dp /)
852! **************************************************************************************************
853 SUBROUTINE pw_nn_smear_r(pw_in, pw_out, coeffs)
854 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in, pw_out
855 REAL(kind=dp), DIMENSION(4), INTENT(in) :: coeffs
856
857 INTEGER :: i, j, k
858 REAL(kind=dp), DIMENSION(-1:1, -1:1, -1:1) :: weights
859
860 DO k = -1, 1
861 DO j = -1, 1
862 DO i = -1, 1
863 weights(i, j, k) = coeffs(abs(i) + abs(j) + abs(k) + 1)
864 END DO
865 END DO
866 END DO
867
868 CALL pw_nn_compose_r(weights=weights, pw_in=pw_in, pw_out=pw_out)
869 END SUBROUTINE pw_nn_smear_r
870
871! **************************************************************************************************
872!> \brief calculates a nearest neighbor central derivative.
873!> for the x dir:
874!> pw_out%array(i,j,k)=( pw_in(i+1,j,k)-pw_in(i-1,j,k) )*coeff(1)+
875!> ( pw_in(i+1,j(+-)1,k)-pw_in(i-1,j(+-)1,k)+
876!> pw_in(i+1,j,k(+-)1)-pw_in(i-1,j,k(+-)1) )*coeff(2)+
877!> ( pw_in(i+1,j(+-)1,k(+-)1)-pw_in(i-1,j(+-)1,k(+-)1)+
878!> pw_in(i+1,j(+-)1,k(-+)1)-pw_in(i-1,j(+-)1,k(-+)1) )*coeff(3)
879!> periodic boundary conditions are applied
880!> \param pw_in the argument for the linear operator
881!> \param pw_out place where the smeared values should be added
882!> \param coeffs array with the coefficent of the front (positive) plane
883!> of the central derivative, ordered with
884!> the distance from the center: coeffs(1) the coeff of the central
885!> element, coeffs(2) the coeff of the 4 element with distance 1,
886!> coeff(3) the coeff of the 4 elements at distance sqrt(2)
887!> \param idir ...
888!> \author Fawzi Mohamed
889!> \note
890!> with coeff=(/ 2.0_dp/9.0_dp, 1.0_dp/18.0_dp, 1.0_dp/72.0_dp /)
891!> is equivalent to pw_spline3_deriv_r, with
892!> coeff=(/ 9.0_dp/32.0_dp, 3.0_dp/64.0_dp, 1.0_dp/128.0_dp /)
893!> to pw_spline2_deriv_r
894!> coeff=(/ 25._dp/72._dp, 5._dp/144, 1._dp/288._dp /)
895! **************************************************************************************************
896 SUBROUTINE pw_nn_deriv_r(pw_in, pw_out, coeffs, idir)
897 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in, pw_out
898 REAL(kind=dp), DIMENSION(3), INTENT(in) :: coeffs
899 INTEGER :: idir
900
901 INTEGER :: i, idirval, j, k
902 REAL(kind=dp), DIMENSION(-1:1, -1:1, -1:1) :: weights
903
904 DO k = -1, 1
905 DO j = -1, 1
906 DO i = -1, 1
907 SELECT CASE (idir)
908 CASE (1)
909 idirval = i
910 CASE (2)
911 idirval = j
912 CASE (3)
913 idirval = k
914 CASE default
915 cpabort("invalid idir ("//trim(cp_to_string(idir))//")")
916 END SELECT
917 IF (idirval == 0) THEN
918 weights(i, j, k) = 0.0_dp
919 ELSE
920 weights(i, j, k) = real(idirval, dp)*coeffs(abs(i) + abs(j) + abs(k))
921 END IF
922 END DO
923 END DO
924 END DO
925
926 CALL pw_nn_compose_r(weights=weights, pw_in=pw_in, pw_out=pw_out)
927 END SUBROUTINE pw_nn_deriv_r
928
929! **************************************************************************************************
930!> \brief low level function that adds a coarse grid
931!> to a fine grid.
932!> If pbc is true periodic boundary conditions are applied
933!>
934!> It will add to
935!>
936!> fine_values(2*coarse_bounds(1,1):2*coarse_bounds(2,1),
937!> 2*coarse_bounds(1,2):2*coarse_bounds(2,2),
938!> 2*coarse_bounds(1,3):2*coarse_bounds(2,3))
939!>
940!> using
941!>
942!> coarse_coeffs(coarse_bounds(1,1):coarse_bounds(2,1),
943!> coarse_bounds(1,2):coarse_bounds(2,2),
944!> coarse_bounds(1,3):coarse_bounds(2,3))
945!>
946!> composed with the weights obtained by the direct product of the
947!> 1d coefficients weights:
948!>
949!> for i,j,k in -3..3
950!> w(i,j,k)=weights_1d(abs(i)+1)*weights_1d(abs(j)+1)*
951!> weights_1d(abs(k)+1)
952!> \param coarse_coeffs_pw the values of the coefficients
953!> \param fine_values_pw where to add the values due to the
954!> coarse coeffs
955!> \param weights_1d the weights of the 1d smearing
956!> \param w_border0 the 1d weight at the border (when pbc is false)
957!> \param w_border1 the 1d weights for a point one off the border
958!> (w_border1(1) is the weight of the coefficent at the border)
959!> (used if pbc is false)
960!> \param pbc if periodic boundary conditions should be applied
961!> \param safe_computation ...
962!> \author fawzi
963!> \note
964!> coarse looping is continuos, I did not check if keeping the fine looping
965!> contiguous is better.
966!> And I ask myself yet again why, why we use x-slice distribution,
967!> z-slice distribution would be much better performancewise
968!> (and would semplify this code enormously).
969!> fine2coarse has much more understandable parallel part (build up of
970!> send/rcv sizes,... but worse if you have really a lot of processors,
971!> probabily irrelevant because it is not critical) [fawzi].
972! **************************************************************************************************
973 SUBROUTINE add_coarse2fine(coarse_coeffs_pw, fine_values_pw, &
974 weights_1d, w_border0, w_border1, pbc, safe_computation)
975 TYPE(pw_r3d_rs_type), INTENT(IN) :: coarse_coeffs_pw, fine_values_pw
976 REAL(kind=dp), DIMENSION(4), INTENT(in) :: weights_1d
977 REAL(kind=dp), INTENT(in) :: w_border0
978 REAL(kind=dp), DIMENSION(3), INTENT(in) :: w_border1
979 LOGICAL, INTENT(in) :: pbc
980 LOGICAL, INTENT(in), OPTIONAL :: safe_computation
981
982 CHARACTER(len=*), PARAMETER :: routinen = 'add_coarse2fine'
983
984 INTEGER :: coarse_slice_size, f_shift(3), fi, fi_lb, fi_ub, fj, fk, handle, handle2, i, ii, &
985 ij, ik, ip, j, k, my_lb, my_ub, n_procs, p, p_lb, p_old, p_ub, rcv_tot_size, rest_b, &
986 s(3), send_tot_size, sf, shift, ss, x, x_att, xx
987 INTEGER, ALLOCATABLE, DIMENSION(:) :: rcv_offset, rcv_size, real_rcv_size, &
988 send_offset, send_size, sent_size
989 INTEGER, DIMENSION(2, 3) :: coarse_bo, coarse_gbo, fine_bo, &
990 fine_gbo, my_coarse_bo
991 INTEGER, DIMENSION(:), POINTER :: pos_of_x
992 LOGICAL :: has_i_lbound, has_i_ubound, is_split, &
993 safe_calc
994 REAL(kind=dp) :: v0, v1, v2, v3, wi, wj, wk
995 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rcv_buf, send_buf
996 REAL(kind=dp), DIMENSION(3) :: w_0, ww0
997 REAL(kind=dp), DIMENSION(4) :: w_1, ww1
998 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: coarse_coeffs, fine_values
999
1000 CALL timeset(routinen, handle)
1001! CALL timeset(routineN//"_pre",handle2)
1002 safe_calc = .false.
1003 IF (PRESENT(safe_computation)) safe_calc = safe_computation
1004 ii = coarse_coeffs_pw%pw_grid%para%group%compare(fine_values_pw%pw_grid%para%group)
1005 IF (ii > mp_comm_congruent) THEN
1006 cpabort("")
1007 END IF
1008 my_coarse_bo = coarse_coeffs_pw%pw_grid%bounds_local
1009 coarse_gbo = coarse_coeffs_pw%pw_grid%bounds
1010 fine_bo = fine_values_pw%pw_grid%bounds_local
1011 fine_gbo = fine_values_pw%pw_grid%bounds
1012 f_shift = fine_gbo(1, :) - 2*coarse_gbo(1, :)
1013 DO j = 2, 3
1014 DO i = 1, 2
1015 coarse_bo(i, j) = floor((fine_bo(i, j) - f_shift(j))/2.)
1016 END DO
1017 END DO
1018 IF (fine_bo(1, 1) <= fine_bo(2, 1)) THEN
1019 coarse_bo(1, 1) = floor((fine_bo(1, 1) - 2 - f_shift(1))/2.)
1020 coarse_bo(2, 1) = floor((fine_bo(2, 1) + 3 - f_shift(1))/2.)
1021 ELSE
1022 coarse_bo(1, 1) = coarse_gbo(2, 1)
1023 coarse_bo(2, 1) = coarse_gbo(2, 1) - 1
1024 END IF
1025 is_split = any(coarse_gbo(:, 1) /= my_coarse_bo(:, 1))
1026 IF (.NOT. is_split .OR. .NOT. pbc) THEN
1027 coarse_bo(1, 1) = max(coarse_gbo(1, 1), coarse_bo(1, 1))
1028 coarse_bo(2, 1) = min(coarse_gbo(2, 1), coarse_bo(2, 1))
1029 END IF
1030 has_i_ubound = (fine_gbo(2, 1) /= fine_bo(2, 1)) .OR. pbc .AND. is_split
1031 has_i_lbound = (fine_gbo(1, 1) /= fine_bo(1, 1)) .OR. pbc .AND. is_split
1032
1033 IF (pbc) THEN
1034 cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1035 cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + 1 + f_shift))
1036 ELSE
1037 cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + f_shift))
1038 cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1039 END IF
1040
1041 coarse_coeffs => coarse_coeffs_pw%array
1042 DO i = 1, 3
1043 s(i) = coarse_gbo(2, i) - coarse_gbo(1, i) + 1
1044 END DO
1045! CALL timestop(handle2)
1046 ! *** parallel case
1047 IF (is_split) THEN
1048 CALL timeset(routinen//"_comm", handle2)
1049 coarse_slice_size = (coarse_bo(2, 2) - coarse_bo(1, 2) + 1)* &
1050 (coarse_bo(2, 3) - coarse_bo(1, 3) + 1)
1051 n_procs = coarse_coeffs_pw%pw_grid%para%group%num_pe
1052 ALLOCATE (send_size(0:n_procs - 1), send_offset(0:n_procs - 1), &
1053 sent_size(0:n_procs - 1), rcv_size(0:n_procs - 1), &
1054 rcv_offset(0:n_procs - 1), real_rcv_size(0:n_procs - 1))
1055
1056 ! ** rcv size count
1057
1058 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
1059 p_old = pos_of_x(coarse_gbo(1, 1) &
1060 + modulo(coarse_bo(1, 1) - coarse_gbo(1, 1), s(1)))
1061 rcv_size = 0
1062 DO x = coarse_bo(1, 1), coarse_bo(2, 1)
1063 p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1064 rcv_size(p) = rcv_size(p) + coarse_slice_size
1065 END DO
1066
1067 ! ** send size count
1068
1069 pos_of_x => fine_values_pw%pw_grid%para%pos_of_x
1070 sf = fine_gbo(2, 1) - fine_gbo(1, 1) + 1
1071 fi_lb = 2*my_coarse_bo(1, 1) - 3 + f_shift(1)
1072 fi_ub = 2*my_coarse_bo(2, 1) + 3 + f_shift(1)
1073 IF (.NOT. pbc) THEN
1074 fi_lb = max(fi_lb, fine_gbo(1, 1))
1075 fi_ub = min(fi_ub, fine_gbo(2, 1))
1076 ELSE
1077 fi_ub = min(fi_ub, fi_lb + sf - 1)
1078 END IF
1079 p_old = pos_of_x(fine_gbo(1, 1) + modulo(fi_lb - fine_gbo(1, 1), sf))
1080 p_lb = floor((fi_lb - 2 - f_shift(1))/2.)
1081 send_size = 0
1082 DO x = fi_lb, fi_ub
1083 p = pos_of_x(fine_gbo(1, 1) + modulo(x - fine_gbo(1, 1), sf))
1084 IF (p /= p_old) THEN
1085 p_ub = floor((x - 1 + 3 - f_shift(1))/2.)
1086
1087 send_size(p_old) = send_size(p_old) + (min(p_ub, my_coarse_bo(2, 1)) &
1088 - max(p_lb, my_coarse_bo(1, 1)) + 1)*coarse_slice_size
1089
1090 IF (pbc) THEN
1091 DO xx = p_lb, coarse_gbo(1, 1) - 1
1092 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1093 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1094 send_size(p_old) = send_size(p_old) + coarse_slice_size
1095 END IF
1096 END DO
1097 DO xx = coarse_gbo(2, 1) + 1, p_ub
1098 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1099 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1100 send_size(p_old) = send_size(p_old) + coarse_slice_size
1101 END IF
1102 END DO
1103 END IF
1104
1105 p_old = p
1106 p_lb = floor((x - 2 - f_shift(1))/2.)
1107 END IF
1108 END DO
1109 p_ub = floor((fi_ub + 3 - f_shift(1))/2.)
1110
1111 send_size(p_old) = send_size(p_old) + (min(p_ub, my_coarse_bo(2, 1)) &
1112 - max(p_lb, my_coarse_bo(1, 1)) + 1)*coarse_slice_size
1113
1114 IF (pbc) THEN
1115 DO xx = p_lb, coarse_gbo(1, 1) - 1
1116 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1117 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1118 send_size(p_old) = send_size(p_old) + coarse_slice_size
1119 END IF
1120 END DO
1121 DO xx = coarse_gbo(2, 1) + 1, p_ub
1122 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1123 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1124 send_size(p_old) = send_size(p_old) + coarse_slice_size
1125 END IF
1126 END DO
1127 END IF
1128 ! ** offsets & alloc send-rcv
1129
1130 send_tot_size = 0
1131 DO ip = 0, n_procs - 1
1132 send_offset(ip) = send_tot_size
1133 send_tot_size = send_tot_size + send_size(ip)
1134 END DO
1135 ALLOCATE (send_buf(0:send_tot_size - 1))
1136
1137 rcv_tot_size = 0
1138 DO ip = 0, n_procs - 1
1139 rcv_offset(ip) = rcv_tot_size
1140 rcv_tot_size = rcv_tot_size + rcv_size(ip)
1141 END DO
1142 IF (.NOT. rcv_tot_size == (coarse_bo(2, 1) - coarse_bo(1, 1) + 1)*coarse_slice_size) THEN
1143 cpabort("Error calculating rcv_tot_size ")
1144 END IF
1145 ALLOCATE (rcv_buf(0:rcv_tot_size - 1))
1146
1147 ! ** fill send buffer
1148
1149 p_old = pos_of_x(fine_gbo(1, 1) + modulo(fi_lb - fine_gbo(1, 1), sf))
1150 p_lb = floor((fi_lb - 2 - f_shift(1))/2.)
1151 sent_size(:) = send_offset
1152 ss = my_coarse_bo(2, 1) - my_coarse_bo(1, 1) + 1
1153 DO x = fi_lb, fi_ub
1154 p = pos_of_x(fine_gbo(1, 1) + modulo(x - fine_gbo(1, 1), sf))
1155 IF (p /= p_old) THEN
1156 shift = floor((fine_gbo(1, 1) + modulo(x - 1 - fine_gbo(1, 1), sf) - f_shift(1))/2._dp) - &
1157 floor((x - 1 - f_shift(1))/2._dp)
1158 p_ub = floor((x - 1 + 3 - f_shift(1))/2._dp)
1159
1160 IF (pbc) THEN
1161 DO xx = p_lb + shift, coarse_gbo(1, 1) - 1
1162 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), sf)
1163 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1164 CALL dcopy(coarse_slice_size, &
1165 coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1166 my_coarse_bo(1, 3)), ss, send_buf(sent_size(p_old)), 1)
1167 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1168 END IF
1169 END DO
1170 END IF
1171
1172 ii = sent_size(p_old)
1173 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1174 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1175 DO i = max(p_lb + shift, my_coarse_bo(1, 1)), min(p_ub + shift, my_coarse_bo(2, 1))
1176 send_buf(ii) = coarse_coeffs(i, j, k)
1177 ii = ii + 1
1178 END DO
1179 END DO
1180 END DO
1181 sent_size(p_old) = ii
1182
1183 IF (pbc) THEN
1184 DO xx = coarse_gbo(2, 1) + 1, p_ub + shift
1185 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1186 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1187 CALL dcopy(coarse_slice_size, &
1188 coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1189 my_coarse_bo(1, 3)), ss, &
1190 send_buf(sent_size(p_old)), 1)
1191 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1192 END IF
1193 END DO
1194 END IF
1195
1196 p_old = p
1197 p_lb = floor((x - 2 - f_shift(1))/2.)
1198 END IF
1199 END DO
1200 shift = floor((fine_gbo(1, 1) + modulo(x - 1 - fine_gbo(1, 1), sf) - f_shift(1))/2._dp) - &
1201 floor((x - 1 - f_shift(1))/2._dp)
1202 p_ub = floor((fi_ub + 3 - f_shift(1))/2.)
1203
1204 IF (pbc) THEN
1205 DO xx = p_lb + shift, coarse_gbo(1, 1) - 1
1206 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1207 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1208 CALL dcopy(coarse_slice_size, &
1209 coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1210 my_coarse_bo(1, 3)), ss, send_buf(sent_size(p_old)), 1)
1211 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1212 END IF
1213 END DO
1214 END IF
1215
1216 ii = sent_size(p_old)
1217 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1218 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1219 DO i = max(p_lb + shift, my_coarse_bo(1, 1)), min(p_ub + shift, my_coarse_bo(2, 1))
1220 send_buf(ii) = coarse_coeffs(i, j, k)
1221 ii = ii + 1
1222 END DO
1223 END DO
1224 END DO
1225 sent_size(p_old) = ii
1226
1227 IF (pbc) THEN
1228 DO xx = coarse_gbo(2, 1) + 1, p_ub + shift
1229 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1230 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1231 CALL dcopy(coarse_slice_size, &
1232 coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1233 my_coarse_bo(1, 3)), ss, send_buf(sent_size(p_old)), 1)
1234 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1235 END IF
1236 END DO
1237 END IF
1238
1239 cpassert(all(sent_size(:n_procs - 2) == send_offset(1:)))
1240 cpassert(sent_size(n_procs - 1) == send_tot_size)
1241 ! test send/rcv sizes
1242 CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(send_size, real_rcv_size, 1)
1243 cpassert(all(real_rcv_size == rcv_size))
1244 ! all2all
1245 CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(sb=send_buf, scount=send_size, sdispl=send_offset, &
1246 rb=rcv_buf, rcount=rcv_size, rdispl=rcv_offset)
1247
1248 ! ** reorder rcv buffer
1249 ! (actually reordering should be needed only with pbc)
1250
1251 ALLOCATE (coarse_coeffs(coarse_bo(1, 1):coarse_bo(2, 1), &
1252 coarse_bo(1, 2):coarse_bo(2, 2), &
1253 coarse_bo(1, 3):coarse_bo(2, 3)))
1254
1255 my_lb = max(coarse_gbo(1, 1), coarse_bo(1, 1))
1256 my_ub = min(coarse_gbo(2, 1), coarse_bo(2, 1))
1257 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
1258 sent_size(:) = rcv_offset
1259 ss = coarse_bo(2, 1) - coarse_bo(1, 1) + 1
1260 DO x = my_ub + 1, coarse_bo(2, 1)
1261 p_old = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1262 CALL dcopy(coarse_slice_size, &
1263 rcv_buf(sent_size(p_old)), 1, &
1264 coarse_coeffs(x, coarse_bo(1, 2), &
1265 coarse_bo(1, 3)), ss)
1266 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1267 END DO
1268 p_old = pos_of_x(coarse_gbo(1, 1) &
1269 + modulo(my_lb - coarse_gbo(1, 1), s(1)))
1270 p_lb = my_lb
1271 DO x = my_lb, my_ub
1272 p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1273 IF (p /= p_old) THEN
1274 p_ub = x - 1
1275
1276 ii = sent_size(p_old)
1277 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1278 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1279 DO i = p_lb, p_ub
1280 coarse_coeffs(i, j, k) = rcv_buf(ii)
1281 ii = ii + 1
1282 END DO
1283 END DO
1284 END DO
1285 sent_size(p_old) = ii
1286
1287 p_lb = x
1288 p_old = p
1289 END IF
1290 rcv_size(p) = rcv_size(p) + coarse_slice_size
1291 END DO
1292 p_ub = my_ub
1293 ii = sent_size(p_old)
1294 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1295 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1296 DO i = p_lb, p_ub
1297 coarse_coeffs(i, j, k) = rcv_buf(ii)
1298 ii = ii + 1
1299 END DO
1300 END DO
1301 END DO
1302 sent_size(p_old) = ii
1303 DO x = coarse_bo(1, 1), my_lb - 1
1304 p_old = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1305 CALL dcopy(coarse_slice_size, &
1306 rcv_buf(sent_size(p_old)), 1, &
1307 coarse_coeffs(x, coarse_bo(1, 2), &
1308 coarse_bo(1, 3)), ss)
1309 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1310 END DO
1311
1312 cpassert(all(sent_size(0:n_procs - 2) == rcv_offset(1:)))
1313 cpassert(sent_size(n_procs - 1) == rcv_tot_size)
1314
1315 ! dealloc
1316 DEALLOCATE (send_size, send_offset, rcv_size, rcv_offset)
1317 DEALLOCATE (send_buf, rcv_buf, real_rcv_size)
1318 CALL timestop(handle2)
1319
1320 END IF
1321 fine_values => fine_values_pw%array
1322 w_0 = (/weights_1d(3), weights_1d(1), weights_1d(3)/)
1323 w_1 = (/weights_1d(4), weights_1d(2), weights_1d(2), weights_1d(4)/)
1324
1325 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1326 DO ik = -3, 3
1327 IF (pbc) THEN
1328 wk = weights_1d(abs(ik) + 1)
1329 fk = fine_gbo(1, 3) + modulo(2*k + ik - fine_gbo(1, 3) + f_shift(3), 2*s(3))
1330 ELSE
1331 fk = 2*k + ik + f_shift(3)
1332 IF (fk <= fine_bo(1, 3) + 1 .OR. fk >= fine_bo(2, 3) - 1) THEN
1333 IF (fk < fine_bo(1, 3) .OR. fk > fine_bo(2, 3)) cycle
1334 IF (fk == fine_bo(1, 3) .OR. fk == fine_bo(2, 3)) THEN
1335 IF (ik /= 0) cycle
1336 wk = w_border0
1337 ELSE IF (fk == 2*coarse_bo(1, 3) + 1 + f_shift(3)) THEN
1338 SELECT CASE (ik)
1339 CASE (1)
1340 wk = w_border1(1)
1341 CASE (-1)
1342 wk = w_border1(2)
1343 CASE (-3)
1344 wk = w_border1(3)
1345 CASE default
1346 cpabort("")
1347 cycle
1348 END SELECT
1349 ELSE
1350 SELECT CASE (ik)
1351 CASE (3)
1352 wk = w_border1(3)
1353 CASE (1)
1354 wk = w_border1(2)
1355 CASE (-1)
1356 wk = w_border1(1)
1357 CASE default
1358 cpabort("")
1359 cycle
1360 END SELECT
1361 END IF
1362 ELSE
1363 wk = weights_1d(abs(ik) + 1)
1364 END IF
1365 END IF
1366 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1367 DO ij = -3, 3
1368 IF (pbc) THEN
1369 wj = weights_1d(abs(ij) + 1)*wk
1370 fj = fine_gbo(1, 2) + modulo(2*j + ij - fine_gbo(1, 2) + f_shift(2), 2*s(2))
1371 ELSE
1372 fj = 2*j + ij + f_shift(2)
1373 IF (fj <= fine_bo(1, 2) + 1 .OR. fj >= fine_bo(2, 2) - 1) THEN
1374 IF (fj < fine_bo(1, 2) .OR. fj > fine_bo(2, 2)) cycle
1375 IF (fj == fine_bo(1, 2) .OR. fj == fine_bo(2, 2)) THEN
1376 IF (ij /= 0) cycle
1377 wj = w_border0*wk
1378 ELSE IF (fj == 2*coarse_bo(1, 2) + 1 + f_shift(2)) THEN
1379 SELECT CASE (ij)
1380 CASE (1)
1381 wj = w_border1(1)*wk
1382 CASE (-1)
1383 wj = w_border1(2)*wk
1384 CASE (-3)
1385 wj = w_border1(3)*wk
1386 CASE default
1387 cycle
1388 END SELECT
1389 ELSE
1390 SELECT CASE (ij)
1391 CASE (-1)
1392 wj = w_border1(1)*wk
1393 CASE (1)
1394 wj = w_border1(2)*wk
1395 CASE (3)
1396 wj = w_border1(3)*wk
1397 CASE default
1398 cycle
1399 END SELECT
1400 END IF
1401 ELSE
1402 wj = weights_1d(abs(ij) + 1)*wk
1403 END IF
1404 END IF
1405
1406 IF (fine_bo(2, 1) - fine_bo(1, 1) < 7 .OR. safe_calc) THEN
1407! CALL timeset(routineN//"_safe",handle2)
1408 DO i = coarse_bo(1, 1), coarse_bo(2, 1)
1409 DO ii = -3, 3
1410 IF (pbc .AND. .NOT. is_split) THEN
1411 wi = weights_1d(abs(ii) + 1)*wj
1412 fi = fine_gbo(1, 1) + modulo(2*i + ii - fine_gbo(1, 1) + f_shift(1), 2*s(1))
1413 ELSE
1414 fi = 2*i + ii + f_shift(1)
1415 IF (fi < fine_bo(1, 1) .OR. fi > fine_bo(2, 1)) cycle
1416 IF (.NOT. pbc .AND. (fi <= fine_gbo(1, 1) + 1 .OR. &
1417 fi >= fine_gbo(2, 1) - 1)) THEN
1418 IF (fi == fine_gbo(1, 1) .OR. fi == fine_gbo(2, 1)) THEN
1419 IF (ii /= 0) cycle
1420 wi = w_border0*wj
1421 ELSE IF (fi == fine_gbo(1, 1) + 1) THEN
1422 SELECT CASE (ii)
1423 CASE (1)
1424 wi = w_border1(1)*wj
1425 CASE (-1)
1426 wi = w_border1(2)*wj
1427 CASE (-3)
1428 wi = w_border1(3)*wj
1429 CASE default
1430 cycle
1431 END SELECT
1432 ELSE
1433 SELECT CASE (ii)
1434 CASE (-1)
1435 wi = w_border1(1)*wj
1436 CASE (1)
1437 wi = w_border1(2)*wj
1438 CASE (3)
1439 wi = w_border1(3)*wj
1440 CASE default
1441 cycle
1442 END SELECT
1443 END IF
1444 ELSE
1445 wi = weights_1d(abs(ii) + 1)*wj
1446 END IF
1447 END IF
1448 fine_values(fi, fj, fk) = &
1449 fine_values(fi, fj, fk) + &
1450 wi*coarse_coeffs(i, j, k)
1451 END DO
1452 END DO
1453! CALL timestop(handle2)
1454 ELSE
1455! CALL timeset(routineN//"_core1",handle2)
1456 ww0 = wj*w_0
1457 ww1 = wj*w_1
1458 IF (pbc .AND. .NOT. is_split) THEN
1459 v3 = coarse_coeffs(coarse_bo(2, 1), j, k)
1460 i = coarse_bo(1, 1)
1461 fi = 2*i + f_shift(1)
1462 v0 = coarse_coeffs(i, j, k)
1463 v1 = coarse_coeffs(i + 1, j, k)
1464 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1465 ww0(1)*v3 + ww0(2)*v0 + ww0(3)*v1
1466 v2 = coarse_coeffs(i + 2, j, k)
1467 fi = fi + 1
1468 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1469 ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1470 ELSE IF (.NOT. has_i_lbound) THEN
1471 i = coarse_bo(1, 1)
1472 fi = 2*i + f_shift(1)
1473 v0 = coarse_coeffs(i, j, k)
1474 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1475 w_border0*wj*v0
1476 v1 = coarse_coeffs(i + 1, j, k)
1477 v2 = coarse_coeffs(i + 2, j, k)
1478 fi = fi + 1
1479 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1480 wj*(w_border1(1)*v0 + w_border1(2)*v1 + &
1481 w_border1(3)*v2)
1482 ELSE
1483 i = coarse_bo(1, 1)
1484 v0 = coarse_coeffs(i, j, k)
1485 v1 = coarse_coeffs(i + 1, j, k)
1486 v2 = coarse_coeffs(i + 2, j, k)
1487 fi = 2*i + f_shift(1) + 1
1488 IF (.NOT. (fi + 1 == fine_bo(1, 1) .OR. &
1489 fi + 2 == fine_bo(1, 1))) THEN
1490 CALL cp_abort(__location__, &
1491 "unexpected start index "// &
1492 trim(cp_to_string(coarse_bo(1, 1)))//" "// &
1493 trim(cp_to_string(fi)))
1494 END IF
1495 END IF
1496 fi = fi + 1
1497 IF (fi >= fine_bo(1, 1)) THEN
1498 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1499 ww0(1)*v0 + ww0(2)*v1 + &
1500 ww0(3)*v2
1501 ELSE
1502 cpassert(fi + 1 == fine_bo(1, 1))
1503 END IF
1504! CALL timestop(handle2)
1505! CALL timeset(routineN//"_core",handle2)
1506 DO i = coarse_bo(1, 1) + 3, floor((fine_bo(2, 1) - f_shift(1))/2.) - 3, 4
1507 v3 = coarse_coeffs(i, j, k)
1508 fi = fi + 1
1509 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1510 (ww1(1)*v0 + ww1(2)*v1 + &
1511 ww1(3)*v2 + ww1(4)*v3)
1512 fi = fi + 1
1513 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1514 (ww0(1)*v1 + ww0(2)*v2 + &
1515 ww0(3)*v3)
1516 v0 = coarse_coeffs(i + 1, j, k)
1517 fi = fi + 1
1518 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1519 (ww1(4)*v0 + ww1(1)*v1 + &
1520 ww1(2)*v2 + ww1(3)*v3)
1521 fi = fi + 1
1522 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1523 (ww0(1)*v2 + ww0(2)*v3 + &
1524 ww0(3)*v0)
1525 v1 = coarse_coeffs(i + 2, j, k)
1526 fi = fi + 1
1527 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1528 (ww1(3)*v0 + ww1(4)*v1 + &
1529 ww1(1)*v2 + ww1(2)*v3)
1530 fi = fi + 1
1531 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1532 (ww0(1)*v3 + ww0(2)*v0 + &
1533 ww0(3)*v1)
1534 v2 = coarse_coeffs(i + 3, j, k)
1535 fi = fi + 1
1536 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1537 (ww1(2)*v0 + ww1(3)*v1 + &
1538 ww1(4)*v2 + ww1(1)*v3)
1539 fi = fi + 1
1540 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1541 (ww0(1)*v0 + ww0(2)*v1 + &
1542 ww0(3)*v2)
1543 END DO
1544! CALL timestop(handle2)
1545! CALL timeset(routineN//"_clean",handle2)
1546 rest_b = modulo(floor((fine_bo(2, 1) - f_shift(1))/2.) - coarse_bo(1, 1) - 3 + 1, 4)
1547 IF (rest_b > 0) THEN
1548 i = floor((fine_bo(2, 1) - f_shift(1))/2.) - rest_b + 1
1549 v3 = coarse_coeffs(i, j, k)
1550 fi = fi + 1
1551 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1552 (ww1(1)*v0 + ww1(2)*v1 + &
1553 ww1(3)*v2 + ww1(4)*v3)
1554 fi = fi + 1
1555 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1556 (ww0(1)*v1 + ww0(2)*v2 + &
1557 ww0(3)*v3)
1558 IF (rest_b > 1) THEN
1559 v0 = coarse_coeffs(i + 1, j, k)
1560 fi = fi + 1
1561 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1562 (ww1(4)*v0 + ww1(1)*v1 + &
1563 ww1(2)*v2 + ww1(3)*v3)
1564 fi = fi + 1
1565 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1566 (ww0(1)*v2 + ww0(2)*v3 + &
1567 ww0(3)*v0)
1568 IF (rest_b > 2) THEN
1569 v1 = coarse_coeffs(i + 2, j, k)
1570 fi = fi + 1
1571 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1572 (ww1(3)*v0 + ww1(4)*v1 + &
1573 ww1(1)*v2 + ww1(2)*v3)
1574 fi = fi + 1
1575 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1576 (ww0(1)*v3 + ww0(2)*v0 + &
1577 ww0(3)*v1)
1578 IF (pbc .AND. .NOT. is_split) THEN
1579 v2 = coarse_coeffs(coarse_bo(1, 1), j, k)
1580 fi = fi + 1
1581 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1582 ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1583 fi = fi + 1
1584 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1585 ww0(1)*v0 + ww0(2)*v1 + ww0(3)*v2
1586 v3 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1587 fi = fi + 1
1588 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1589 ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1590 ELSE IF (has_i_ubound) THEN
1591 v2 = coarse_coeffs(i + 3, j, k)
1592 fi = fi + 1
1593 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1594 ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1595 fi = fi + 1
1596 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1597 ww0(1)*v0 + ww0(2)*v1 + ww0(3)*v2
1598 IF (fi + 1 == fine_bo(2, 1)) THEN
1599 v3 = coarse_coeffs(i + 4, j, k)
1600 fi = fi + 1
1601 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1602 ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1603 END IF
1604 ELSE
1605 fi = fi + 1
1606 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1607 wj*(w_border1(3)*v3 + w_border1(2)*v0 + &
1608 w_border1(1)*v1)
1609 fi = fi + 1
1610 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1611 w_border0*wj*v1
1612 END IF
1613 ELSE IF (pbc .AND. .NOT. is_split) THEN
1614 v1 = coarse_coeffs(coarse_bo(1, 1), j, k)
1615 fi = fi + 1
1616 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1617 ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1618 fi = fi + 1
1619 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1620 ww0(1)*v3 + ww0(2)*v0 + ww0(3)*v1
1621 v2 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1622 fi = fi + 1
1623 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1624 ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1625 ELSE IF (has_i_ubound) THEN
1626 v1 = coarse_coeffs(i + 2, j, k)
1627 fi = fi + 1
1628 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1629 ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1630 fi = fi + 1
1631 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1632 ww0(1)*v3 + ww0(2)*v0 + ww0(3)*v1
1633 IF (fi + 1 == fine_bo(2, 1)) THEN
1634 v2 = coarse_coeffs(i + 3, j, k)
1635 fi = fi + 1
1636 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1637 ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1638 END IF
1639 ELSE
1640 fi = fi + 1
1641 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1642 wj*(w_border1(3)*v2 + w_border1(2)*v3 + &
1643 w_border1(1)*v0)
1644 fi = fi + 1
1645 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1646 w_border0*wj*v0
1647 END IF
1648 ELSE IF (pbc .AND. .NOT. is_split) THEN
1649 v0 = coarse_coeffs(coarse_bo(1, 1), j, k)
1650 fi = fi + 1
1651 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1652 ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1653 fi = fi + 1
1654 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1655 ww0(1)*v2 + ww0(2)*v3 + ww0(3)*v0
1656 v1 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1657 fi = fi + 1
1658 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1659 ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1660 ELSE IF (has_i_ubound) THEN
1661 v0 = coarse_coeffs(i + 1, j, k)
1662 fi = fi + 1
1663 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1664 ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1665 fi = fi + 1
1666 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1667 ww0(1)*v2 + ww0(2)*v3 + ww0(3)*v0
1668 IF (fi + 1 == fine_bo(2, 1)) THEN
1669 v1 = coarse_coeffs(i + 2, j, k)
1670 fi = fi + 1
1671 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1672 ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1673 END IF
1674 ELSE
1675 fi = fi + 1
1676 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1677 wj*(w_border1(3)*v1 + w_border1(2)*v2 + &
1678 w_border1(1)*v3)
1679 fi = fi + 1
1680 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1681 w_border0*wj*v3
1682 END IF
1683 ELSE IF (pbc .AND. .NOT. is_split) THEN
1684 v3 = coarse_coeffs(coarse_bo(1, 1), j, k)
1685 fi = fi + 1
1686 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1687 ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1688 fi = fi + 1
1689 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1690 ww0(1)*v1 + ww0(2)*v2 + ww0(3)*v3
1691 v0 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1692 fi = fi + 1
1693 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1694 ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1695 ELSE IF (has_i_ubound) THEN
1696 v3 = coarse_coeffs(i, j, k)
1697 fi = fi + 1
1698 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1699 ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1700 fi = fi + 1
1701 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1702 ww0(1)*v1 + ww0(2)*v2 + ww0(3)*v3
1703 IF (fi + 1 == fine_bo(2, 1)) THEN
1704 v0 = coarse_coeffs(i + 1, j, k)
1705 fi = fi + 1
1706 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1707 ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1708 END IF
1709 ELSE
1710 fi = fi + 1
1711 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1712 wj*(w_border1(3)*v0 + w_border1(2)*v1 + &
1713 w_border1(1)*v2)
1714 fi = fi + 1
1715 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1716 w_border0*wj*v2
1717 END IF
1718 cpassert(fi == fine_bo(2, 1))
1719 END IF
1720! CALL timestop(handle2)
1721 END DO
1722 END DO
1723 END DO
1724 END DO
1725
1726 IF (is_split) THEN
1727 DEALLOCATE (coarse_coeffs)
1728 END IF
1729 CALL timestop(handle)
1730 END SUBROUTINE add_coarse2fine
1731
1732! **************************************************************************************************
1733!> \brief low level function that adds a coarse grid (without boundary)
1734!> to a fine grid.
1735!>
1736!> It will add to
1737!>
1738!> coarse_coeffs(coarse_bounds(1,1):coarse_bounds(2,1),
1739!> coarse_bounds(1,2):coarse_bounds(2,2),
1740!> coarse_bounds(1,3):coarse_bounds(2,3))
1741!>
1742!> using
1743!>
1744!> fine_values(2*coarse_bounds(1,1):2*coarse_bounds(2,1),
1745!> 2*coarse_bounds(1,2):2*coarse_bounds(2,2),
1746!> 2*coarse_bounds(1,3):2*coarse_bounds(2,3))
1747!>
1748!> composed with the weights obtained by the direct product of the
1749!> 1d coefficients weights:
1750!>
1751!> for i,j,k in -3..3
1752!> w(i,j,k)=weights_1d(abs(i)+1)*weights_1d(abs(j)+1)*
1753!> weights_1d(abs(k)+1)
1754!> \param fine_values_pw 3d array where to add the values due to the
1755!> coarse coeffs
1756!> \param coarse_coeffs_pw 3d array with boundary of size 1 with the values of the
1757!> coefficients
1758!> \param weights_1d the weights of the 1d smearing
1759!> \param w_border0 the 1d weight at the border
1760!> \param w_border1 the 1d weights for a point one off the border
1761!> (w_border1(1) is the weight of the coefficent at the border)
1762!> \param pbc ...
1763!> \param safe_computation ...
1764!> \author fawzi
1765!> \note
1766!> see coarse2fine for some relevant notes
1767! **************************************************************************************************
1768 SUBROUTINE add_fine2coarse(fine_values_pw, coarse_coeffs_pw, &
1769 weights_1d, w_border0, w_border1, pbc, safe_computation)
1770 TYPE(pw_r3d_rs_type), INTENT(IN) :: fine_values_pw, coarse_coeffs_pw
1771 REAL(kind=dp), DIMENSION(4), INTENT(in) :: weights_1d
1772 REAL(kind=dp), INTENT(in) :: w_border0
1773 REAL(kind=dp), DIMENSION(3), INTENT(in) :: w_border1
1774 LOGICAL, INTENT(in) :: pbc
1775 LOGICAL, INTENT(in), OPTIONAL :: safe_computation
1776
1777 CHARACTER(len=*), PARAMETER :: routinen = 'add_fine2coarse'
1778
1779 INTEGER :: coarse_slice_size, f_shift(3), fi, fj, fk, handle, handle2, i, ii, ij, ik, ip, j, &
1780 k, n_procs, p, p_old, rcv_tot_size, rest_b, s(3), send_tot_size, ss, x, x_att
1781 INTEGER, ALLOCATABLE, DIMENSION(:) :: pp_lb, pp_ub, rcv_offset, rcv_size, &
1782 real_rcv_size, send_offset, send_size, &
1783 sent_size
1784 INTEGER, DIMENSION(2, 3) :: coarse_bo, coarse_gbo, fine_bo, &
1785 fine_gbo, my_coarse_bo
1786 INTEGER, DIMENSION(:), POINTER :: pos_of_x
1787 LOGICAL :: has_i_lbound, has_i_ubound, is_split, &
1788 local_data, safe_calc
1789 REAL(kind=dp) :: vv0, vv1, vv2, vv3, vv4, vv5, vv6, vv7, &
1790 wi, wj, wk
1791 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rcv_buf, send_buf
1792 REAL(kind=dp), DIMENSION(3) :: w_0, ww0
1793 REAL(kind=dp), DIMENSION(4) :: w_1, ww1
1794 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: coarse_coeffs, fine_values
1795
1796 CALL timeset(routinen, handle)
1797
1798 safe_calc = .false.
1799 IF (PRESENT(safe_computation)) safe_calc = safe_computation
1800
1801 my_coarse_bo = coarse_coeffs_pw%pw_grid%bounds_local
1802 coarse_gbo = coarse_coeffs_pw%pw_grid%bounds
1803 fine_bo = fine_values_pw%pw_grid%bounds_local
1804 fine_gbo = fine_values_pw%pw_grid%bounds
1805 f_shift = fine_gbo(1, :) - 2*coarse_gbo(1, :)
1806 is_split = any(coarse_gbo(:, 1) /= my_coarse_bo(:, 1))
1807 coarse_bo = my_coarse_bo
1808 IF (fine_bo(1, 1) <= fine_bo(2, 1)) THEN
1809 coarse_bo(1, 1) = floor(real(fine_bo(1, 1) - f_shift(1), dp)/2._dp) - 1
1810 coarse_bo(2, 1) = floor(real(fine_bo(2, 1) + 1 - f_shift(1), dp)/2._dp) + 1
1811 ELSE
1812 coarse_bo(1, 1) = coarse_gbo(2, 1)
1813 coarse_bo(2, 1) = coarse_gbo(2, 1) - 1
1814 END IF
1815 IF (.NOT. is_split .OR. .NOT. pbc) THEN
1816 coarse_bo(1, 1) = max(coarse_gbo(1, 1), coarse_bo(1, 1))
1817 coarse_bo(2, 1) = min(coarse_gbo(2, 1), coarse_bo(2, 1))
1818 END IF
1819 has_i_ubound = (fine_gbo(2, 1) /= fine_bo(2, 1)) .OR. pbc .AND. is_split
1820 has_i_lbound = (fine_gbo(1, 1) /= fine_bo(1, 1)) .OR. pbc .AND. is_split
1821
1822 IF (pbc) THEN
1823 cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1824 cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + f_shift + 1))
1825 ELSE
1826 cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + f_shift))
1827 cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1828 END IF
1829 cpassert(coarse_gbo(2, 1) - coarse_gbo(1, 2) > 1)
1830 local_data = is_split ! ANY(coarse_bo/=my_coarse_bo)
1831 IF (local_data) THEN
1832 ALLOCATE (coarse_coeffs(coarse_bo(1, 1):coarse_bo(2, 1), &
1833 coarse_bo(1, 2):coarse_bo(2, 2), &
1834 coarse_bo(1, 3):coarse_bo(2, 3)))
1835 coarse_coeffs = 0._dp
1836 ELSE
1837 coarse_coeffs => coarse_coeffs_pw%array
1838 END IF
1839
1840 fine_values => fine_values_pw%array
1841 w_0 = (/weights_1d(3), weights_1d(1), weights_1d(3)/)
1842 w_1 = (/weights_1d(4), weights_1d(2), weights_1d(2), weights_1d(4)/)
1843
1844 DO i = 1, 3
1845 s(i) = coarse_gbo(2, i) - coarse_gbo(1, i) + 1
1846 END DO
1847 IF (any(s < 1)) RETURN
1848
1849 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1850 DO ik = -3, 3
1851 IF (pbc) THEN
1852 wk = weights_1d(abs(ik) + 1)
1853 fk = fine_gbo(1, 3) + modulo(2*k + ik - fine_gbo(1, 3) + f_shift(3), 2*s(3))
1854 ELSE
1855 fk = 2*k + ik + f_shift(3)
1856 IF (fk <= fine_bo(1, 3) + 1 .OR. fk >= fine_bo(2, 3) - 1) THEN
1857 IF (fk < fine_bo(1, 3) .OR. fk > fine_bo(2, 3)) cycle
1858 IF (fk == fine_bo(1, 3) .OR. fk == fine_bo(2, 3)) THEN
1859 IF (ik /= 0) cycle
1860 wk = w_border0
1861 ELSE IF (fk == fine_bo(1, 3) + 1) THEN
1862 SELECT CASE (ik)
1863 CASE (1)
1864 wk = w_border1(1)
1865 CASE (-1)
1866 wk = w_border1(2)
1867 CASE (-3)
1868 wk = w_border1(3)
1869 CASE default
1870 cpabort("")
1871 cycle
1872 END SELECT
1873 ELSE
1874 SELECT CASE (ik)
1875 CASE (3)
1876 wk = w_border1(3)
1877 CASE (1)
1878 wk = w_border1(2)
1879 CASE (-1)
1880 wk = w_border1(1)
1881 CASE default
1882 cpabort("")
1883 cycle
1884 END SELECT
1885 END IF
1886 ELSE
1887 wk = weights_1d(abs(ik) + 1)
1888 END IF
1889 END IF
1890 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1891 DO ij = -3, 3
1892 IF (pbc) THEN
1893 fj = fine_gbo(1, 2) + modulo(2*j + ij - fine_gbo(1, 2) + f_shift(2), &
1894 2*s(2))
1895 wj = weights_1d(abs(ij) + 1)*wk
1896 ELSE
1897 fj = 2*j + ij + f_shift(2)
1898 IF (fj <= fine_bo(1, 2) + 1 .OR. fj >= fine_bo(2, 2) - 1) THEN
1899 IF (fj < fine_bo(1, 2) .OR. fj > fine_bo(2, 2)) cycle
1900 IF (fj == fine_bo(1, 2) .OR. fj == fine_bo(2, 2)) THEN
1901 IF (ij /= 0) cycle
1902 wj = w_border0*wk
1903 ELSE IF (fj == fine_bo(1, 2) + 1) THEN
1904 SELECT CASE (ij)
1905 CASE (1)
1906 wj = w_border1(1)*wk
1907 CASE (-1)
1908 wj = w_border1(2)*wk
1909 CASE (-3)
1910 wj = w_border1(3)*wk
1911 CASE default
1912 cpabort("")
1913 cycle
1914 END SELECT
1915 ELSE
1916 SELECT CASE (ij)
1917 CASE (-1)
1918 wj = w_border1(1)*wk
1919 CASE (1)
1920 wj = w_border1(2)*wk
1921 CASE (3)
1922 wj = w_border1(3)*wk
1923 CASE default
1924 cpabort("")
1925 cycle
1926 END SELECT
1927 END IF
1928 ELSE
1929 wj = weights_1d(abs(ij) + 1)*wk
1930 END IF
1931 END IF
1932
1933 IF (coarse_bo(2, 1) - coarse_bo(1, 1) < 7 .OR. safe_calc) THEN
1934 DO i = coarse_bo(1, 1), coarse_bo(2, 1)
1935 DO ii = -3, 3
1936 IF (pbc .AND. .NOT. is_split) THEN
1937 wi = weights_1d(abs(ii) + 1)*wj
1938 fi = fine_gbo(1, 1) + modulo(2*i + ii - fine_gbo(1, 1) + f_shift(1), 2*s(1))
1939 ELSE
1940 fi = 2*i + ii + f_shift(1)
1941 IF (fi < fine_bo(1, 1) .OR. fi > fine_bo(2, 1)) cycle
1942 IF (((.NOT. pbc) .AND. fi <= fine_gbo(1, 1) + 1) .OR. &
1943 ((.NOT. pbc) .AND. fi >= fine_gbo(2, 1) - 1)) THEN
1944 IF (fi == fine_gbo(1, 1) .OR. fi == fine_gbo(2, 1)) THEN
1945 IF (ii /= 0) cycle
1946 wi = w_border0*wj
1947 ELSE IF (fi == fine_gbo(1, 1) + 1) THEN
1948 SELECT CASE (ii)
1949 CASE (1)
1950 wi = w_border1(1)*wj
1951 CASE (-1)
1952 wi = w_border1(2)*wj
1953 CASE (-3)
1954 wi = w_border1(3)*wj
1955 CASE default
1956 cycle
1957 END SELECT
1958 ELSE
1959 SELECT CASE (ii)
1960 CASE (-1)
1961 wi = w_border1(1)*wj
1962 CASE (1)
1963 wi = w_border1(2)*wj
1964 CASE (3)
1965 wi = w_border1(3)*wj
1966 CASE default
1967 cycle
1968 END SELECT
1969 END IF
1970 ELSE
1971 wi = weights_1d(abs(ii) + 1)*wj
1972 END IF
1973 END IF
1974 coarse_coeffs(i, j, k) = &
1975 coarse_coeffs(i, j, k) + &
1976 wi*fine_values(fi, fj, fk)
1977 END DO
1978 END DO
1979 ELSE
1980 ww0 = wj*w_0
1981 ww1 = wj*w_1
1982 IF (pbc .AND. .NOT. is_split) THEN
1983 i = coarse_bo(1, 1) - 1
1984 vv2 = fine_values(fine_bo(2, 1) - 2, fj, fk)
1985 vv3 = fine_values(fine_bo(2, 1) - 1, fj, fk)
1986 vv4 = fine_values(fine_bo(2, 1), fj, fk)
1987 fi = fine_bo(1, 1)
1988 vv5 = fine_values(fi, fj, fk)
1989 fi = fi + 1
1990 vv6 = fine_values(fi, fj, fk)
1991 fi = fi + 1
1992 vv7 = fine_values(fi, fj, fk)
1993 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
1994 + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7
1995 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
1996 + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7
1997 coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
1998 + ww1(4)*vv6 + ww0(3)*vv7
1999 ELSE IF (has_i_lbound) THEN
2000 i = coarse_bo(1, 1)
2001 fi = fine_bo(1, 1) - 1
2002 IF (i + 1 == floor((fine_bo(1, 1) + 1 - f_shift(1))/2._dp)) THEN
2003 fi = fi + 1
2004 vv0 = fine_values(fi, fj, fk)
2005 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) + &
2006 vv0*ww0(3)
2007 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) + &
2008 vv0*ww0(2)
2009 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) + &
2010 vv0*ww0(1)
2011 END IF
2012 ELSE
2013 i = coarse_bo(1, 1)
2014 fi = 2*i + f_shift(1)
2015 vv0 = fine_values(fi, fj, fk)
2016 fi = fi + 1
2017 vv1 = fine_values(fi, fj, fk)
2018 fi = fi + 1
2019 vv2 = fine_values(fi, fj, fk)
2020 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) + &
2021 (vv0*w_border0 + vv1*w_border1(1))*wj + vv2*ww0(1)
2022 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) + &
2023 wj*w_border1(2)*vv1 + ww0(2)*vv2
2024 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) + &
2025 wj*w_border1(3)*vv1 + ww0(3)*vv2
2026 END IF
2027 DO i = coarse_bo(1, 1) + 3, floor((fine_bo(2, 1) - f_shift(1))/2._dp) - 3, 4
2028 fi = fi + 1
2029 vv0 = fine_values(fi, fj, fk)
2030 fi = fi + 1
2031 vv1 = fine_values(fi, fj, fk)
2032 fi = fi + 1
2033 vv2 = fine_values(fi, fj, fk)
2034 fi = fi + 1
2035 vv3 = fine_values(fi, fj, fk)
2036 fi = fi + 1
2037 vv4 = fine_values(fi, fj, fk)
2038 fi = fi + 1
2039 vv5 = fine_values(fi, fj, fk)
2040 fi = fi + 1
2041 vv6 = fine_values(fi, fj, fk)
2042 fi = fi + 1
2043 vv7 = fine_values(fi, fj, fk)
2044 coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2045 + ww1(1)*vv0
2046 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2047 + ww1(2)*vv0 + ww0(1)*vv1 + ww1(1)*vv2
2048 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2049 + ww1(3)*vv0 + ww0(2)*vv1 + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2050 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2051 + ww1(4)*vv0 + ww0(3)*vv1 + ww1(3)*vv2 + ww0(2)*vv3 + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2052 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2053 + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7
2054 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2055 + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7
2056 coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
2057 + ww1(4)*vv6 + ww0(3)*vv7
2058 END DO
2059 IF (.NOT. floor((fine_bo(2, 1) - f_shift(1))/2._dp) - coarse_bo(1, 1) >= 4) THEN
2060 cpabort("FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)-coarse_bo(1,1)>=4")
2061 END IF
2062 rest_b = modulo(floor((fine_bo(2, 1) - f_shift(1))/2._dp) - coarse_bo(1, 1) - 6, 4)
2063 i = floor((fine_bo(2, 1) - f_shift(1))/2._dp) - 3 - rest_b + 4
2064 cpassert(fi == (i - 2)*2 + f_shift(1))
2065 IF (rest_b > 0) THEN
2066 fi = fi + 1
2067 vv0 = fine_values(fi, fj, fk)
2068 fi = fi + 1
2069 vv1 = fine_values(fi, fj, fk)
2070 coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2071 + ww1(1)*vv0
2072 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2073 + ww1(2)*vv0 + ww0(1)*vv1
2074 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2075 + ww1(3)*vv0 + ww0(2)*vv1
2076 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2077 + ww1(4)*vv0 + ww0(3)*vv1
2078 IF (rest_b > 1) THEN
2079 fi = fi + 1
2080 vv2 = fine_values(fi, fj, fk)
2081 fi = fi + 1
2082 vv3 = fine_values(fi, fj, fk)
2083 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2084 + ww1(1)*vv2
2085 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2086 + ww1(2)*vv2 + ww0(1)*vv3
2087 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2088 + ww1(3)*vv2 + ww0(2)*vv3
2089 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2090 + ww1(4)*vv2 + ww0(3)*vv3
2091 IF (rest_b > 2) THEN
2092 fi = fi + 1
2093 vv4 = fine_values(fi, fj, fk)
2094 fi = fi + 1
2095 vv5 = fine_values(fi, fj, fk)
2096 fi = fi + 1
2097 vv6 = fine_values(fi, fj, fk)
2098 fi = fi + 1
2099 vv7 = fine_values(fi, fj, fk)
2100 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2101 + ww1(1)*vv4
2102 IF (has_i_ubound) THEN
2103 IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2104 fi = fi + 1
2105 vv0 = fine_values(fi, fj, fk)
2106 coarse_coeffs(i + 4, j, k) = coarse_coeffs(i + 4, j, k) &
2107 + vv0*ww1(4)
2108 ELSE
2109 vv0 = 0._dp
2110 END IF
2111 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2112 + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2113 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2114 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7 + vv0*ww1(1)
2115 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2116 + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7 + vv0*ww1(2)
2117 coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
2118 + ww1(4)*vv6 + ww0(3)*vv7 + vv0*ww1(3)
2119 ELSEIF (pbc .AND. .NOT. is_split) THEN
2120 fi = fi + 1
2121 vv0 = fine_values(fi, fj, fk)
2122 vv1 = fine_values(fine_bo(1, 1), fj, fk)
2123 vv2 = fine_values(fine_bo(1, 1) + 1, fj, fk)
2124 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2125 + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2126 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2127 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7 + vv0*ww1(1)
2128 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2129 + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7 + vv0*ww1(2) &
2130 + vv1*ww0(1) + vv2*ww1(1)
2131 ELSE
2132 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2133 + ww1(2)*vv4 + ww0(1)*vv5 + wj*w_border1(3)*vv6
2134 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2135 + ww1(3)*vv4 + ww0(2)*vv5 + wj*w_border1(2)*vv6
2136 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2137 + ww1(4)*vv4 + ww0(3)*vv5 + wj*w_border1(1)*vv6 + w_border0*wj*vv7
2138 END IF
2139 ELSE
2140 fi = fi + 1
2141 vv4 = fine_values(fi, fj, fk)
2142 fi = fi + 1
2143 vv5 = fine_values(fi, fj, fk)
2144 IF (has_i_ubound) THEN
2145 IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2146 fi = fi + 1
2147 vv6 = fine_values(fi, fj, fk)
2148 coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
2149 + ww1(4)*vv6
2150 ELSE
2151 vv6 = 0._dp
2152 END IF
2153 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2154 + ww1(1)*vv4
2155 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2156 + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2157 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2158 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6
2159 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2160 + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6
2161 ELSEIF (pbc .AND. .NOT. is_split) THEN
2162 fi = fi + 1
2163 vv6 = fine_values(fi, fj, fk)
2164 vv7 = fine_values(fine_bo(1, 1), fj, fk)
2165 vv0 = fine_values(fine_bo(1, 1) + 1, fj, fk)
2166 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2167 + ww1(1)*vv4
2168 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2169 + ww1(4)*vv0 + ww0(3)*vv1 + ww1(3)*vv2 + ww0(2)*vv3 + &
2170 ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2171 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2172 + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 &
2173 + ww0(1)*vv7 + ww1(1)*vv0
2174 ELSE
2175 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2176 + wj*w_border1(3)*vv4
2177 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2178 + wj*w_border1(2)*vv4
2179 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2180 + wj*(w_border1(1)*vv4 + w_border0*vv5)
2181 END IF
2182 END IF
2183 ELSE
2184 fi = fi + 1
2185 vv2 = fine_values(fi, fj, fk)
2186 fi = fi + 1
2187 vv3 = fine_values(fi, fj, fk)
2188 IF (has_i_ubound) THEN
2189 IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2190 fi = fi + 1
2191 vv4 = fine_values(fi, fj, fk)
2192 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2193 + ww1(4)*vv4
2194 ELSE
2195 vv4 = 0._dp
2196 END IF
2197 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2198 + ww1(1)*vv2
2199 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2200 + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2201 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2202 + ww1(3)*vv2 + ww0(2)*vv3 + ww1(2)*vv4
2203 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2204 + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4
2205 ELSEIF (pbc .AND. .NOT. is_split) THEN
2206 fi = fi + 1
2207 vv4 = fine_values(fi, fj, fk)
2208 vv5 = fine_values(fine_bo(1, 1), fj, fk)
2209 vv6 = fine_values(fine_bo(1, 1) + 1, fj, fk)
2210 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2211 + ww1(1)*vv2
2212 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2213 + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2214 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2215 + ww1(3)*vv2 + ww0(2)*vv3 + ww1(2)*vv4 + vv5*ww0(1) + ww1(1)*vv6
2216 ELSE
2217 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2218 + wj*w_border1(3)*vv2
2219 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2220 + wj*w_border1(2)*vv2
2221 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2222 + wj*(w_border1(1)*vv2 + w_border0*vv3)
2223 END IF
2224 END IF
2225 ELSE
2226 fi = fi + 1
2227 vv0 = fine_values(fi, fj, fk)
2228 fi = fi + 1
2229 vv1 = fine_values(fi, fj, fk)
2230 IF (has_i_ubound) THEN
2231 IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2232 fi = fi + 1
2233 vv2 = fine_values(fi, fj, fk)
2234 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2235 + ww1(4)*vv2
2236 ELSE
2237 vv2 = 0._dp
2238 END IF
2239 coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2240 + ww1(1)*vv0
2241 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2242 + ww1(2)*vv0 + ww0(1)*vv1 + ww1(1)*vv2
2243 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2244 + ww1(3)*vv0 + ww0(2)*vv1 + ww1(2)*vv2
2245 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2246 + ww1(4)*vv0 + ww0(3)*vv1 + ww1(3)*vv2
2247 ELSEIF (pbc .AND. .NOT. is_split) THEN
2248 fi = fi + 1
2249 vv2 = fine_values(fi, fj, fk)
2250 vv3 = fine_values(fine_bo(1, 1), fk, fk)
2251 vv4 = fine_values(fine_bo(1, 1) + 1, fk, fk)
2252 coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2253 + ww1(1)*vv0
2254 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2255 + ww1(2)*vv0 + ww0(1)*vv1 + ww1(1)*vv2
2256 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2257 + ww1(3)*vv0 + ww0(2)*vv1 + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2258 ELSE
2259 coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2260 + wj*w_border1(3)*vv0
2261 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2262 + wj*w_border1(2)*vv0
2263 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2264 + wj*(w_border1(1)*vv0 + w_border0*vv1)
2265 END IF
2266 END IF
2267 cpassert(fi == fine_bo(2, 1))
2268 END IF
2269 END DO
2270 END DO
2271 END DO
2272 END DO
2273
2274 ! *** parallel case
2275 IF (is_split) THEN
2276 CALL timeset(routinen//"_comm", handle2)
2277 coarse_slice_size = (coarse_bo(2, 2) - coarse_bo(1, 2) + 1)* &
2278 (coarse_bo(2, 3) - coarse_bo(1, 3) + 1)
2279 n_procs = coarse_coeffs_pw%pw_grid%para%group%num_pe
2280 ALLOCATE (send_size(0:n_procs - 1), send_offset(0:n_procs - 1), &
2281 sent_size(0:n_procs - 1), rcv_size(0:n_procs - 1), &
2282 rcv_offset(0:n_procs - 1), pp_lb(0:n_procs - 1), &
2283 pp_ub(0:n_procs - 1), real_rcv_size(0:n_procs - 1))
2284
2285 ! ** send size count
2286
2287 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
2288 send_size = 0
2289 DO x = coarse_bo(1, 1), coarse_bo(2, 1)
2290 p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
2291 send_size(p) = send_size(p) + coarse_slice_size
2292 END DO
2293
2294 ! ** rcv size count
2295
2296 pos_of_x => fine_values_pw%pw_grid%para%pos_of_x
2297 p_old = pos_of_x(fine_gbo(1, 1))
2298 pp_lb = fine_gbo(2, 1)
2299 pp_ub = fine_gbo(2, 1) - 1
2300 pp_lb(p_old) = fine_gbo(1, 1)
2301 DO x = fine_gbo(1, 1), fine_gbo(2, 1)
2302 p = pos_of_x(x)
2303 IF (p /= p_old) THEN
2304 pp_ub(p_old) = x - 1
2305 pp_lb(p) = x
2306 p_old = p
2307 END IF
2308 END DO
2309 pp_ub(p_old) = fine_gbo(2, 1)
2310
2311 DO ip = 0, n_procs - 1
2312 IF (pp_lb(ip) <= pp_ub(ip)) THEN
2313 pp_lb(ip) = floor(real(pp_lb(ip) - f_shift(1), dp)/2._dp) - 1
2314 pp_ub(ip) = floor(real(pp_ub(ip) + 1 - f_shift(1), dp)/2._dp) + 1
2315 ELSE
2316 pp_lb(ip) = coarse_gbo(2, 1)
2317 pp_ub(ip) = coarse_gbo(2, 1) - 1
2318 END IF
2319 IF (.NOT. is_split .OR. .NOT. pbc) THEN
2320 pp_lb(ip) = max(pp_lb(ip), coarse_gbo(1, 1))
2321 pp_ub(ip) = min(pp_ub(ip), coarse_gbo(2, 1))
2322 END IF
2323 END DO
2324
2325 rcv_size = 0
2326 DO ip = 0, n_procs - 1
2327 DO x = pp_lb(ip), coarse_gbo(1, 1) - 1
2328 x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2329 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2330 rcv_size(ip) = rcv_size(ip) + coarse_slice_size
2331 END IF
2332 END DO
2333 rcv_size(ip) = rcv_size(ip) + coarse_slice_size* &
2334 max(0, &
2335 min(pp_ub(ip), my_coarse_bo(2, 1)) - max(pp_lb(ip), my_coarse_bo(1, 1)) + 1)
2336 DO x = coarse_gbo(2, 1) + 1, pp_ub(ip)
2337 x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2338 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2339 rcv_size(ip) = rcv_size(ip) + coarse_slice_size
2340 END IF
2341 END DO
2342 END DO
2343
2344 ! ** offsets & alloc send-rcv
2345
2346 send_tot_size = 0
2347 DO ip = 0, n_procs - 1
2348 send_offset(ip) = send_tot_size
2349 send_tot_size = send_tot_size + send_size(ip)
2350 END DO
2351 IF (send_tot_size /= (coarse_bo(2, 1) - coarse_bo(1, 1) + 1)*coarse_slice_size) &
2352 cpabort("Error calculating send_tot_size")
2353 ALLOCATE (send_buf(0:send_tot_size - 1))
2354
2355 rcv_tot_size = 0
2356 DO ip = 0, n_procs - 1
2357 rcv_offset(ip) = rcv_tot_size
2358 rcv_tot_size = rcv_tot_size + rcv_size(ip)
2359 END DO
2360 ALLOCATE (rcv_buf(0:rcv_tot_size - 1))
2361
2362 ! ** fill send buffer
2363
2364 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
2365 p_old = pos_of_x(coarse_gbo(1, 1) &
2366 + modulo(coarse_bo(1, 1) - coarse_gbo(1, 1), s(1)))
2367 sent_size(:) = send_offset
2368 ss = coarse_bo(2, 1) - coarse_bo(1, 1) + 1
2369 DO x = coarse_bo(1, 1), coarse_bo(2, 1)
2370 p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
2371 CALL dcopy(coarse_slice_size, &
2372 coarse_coeffs(x, coarse_bo(1, 2), &
2373 coarse_bo(1, 3)), ss, send_buf(sent_size(p)), 1)
2374 sent_size(p) = sent_size(p) + coarse_slice_size
2375 END DO
2376
2377 IF (any(sent_size(0:n_procs - 2) /= send_offset(1:n_procs - 1))) &
2378 cpabort("error 1 filling send buffer")
2379 IF (sent_size(n_procs - 1) /= send_tot_size) &
2380 cpabort("error 2 filling send buffer")
2381
2382 IF (local_data) THEN
2383 DEALLOCATE (coarse_coeffs)
2384 ELSE
2385 NULLIFY (coarse_coeffs)
2386 END IF
2387
2388 cpassert(all(sent_size(:n_procs - 2) == send_offset(1:)))
2389 cpassert(sent_size(n_procs - 1) == send_tot_size)
2390 ! test send/rcv sizes
2391 CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(send_size, real_rcv_size, 1)
2392
2393 cpassert(all(real_rcv_size == rcv_size))
2394 ! all2all
2395 CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(sb=send_buf, scount=send_size, sdispl=send_offset, &
2396 rb=rcv_buf, rcount=rcv_size, rdispl=rcv_offset)
2397
2398 ! ** sum & reorder rcv buffer
2399
2400 sent_size(:) = rcv_offset
2401 DO ip = 0, n_procs - 1
2402
2403 DO x = pp_lb(ip), coarse_gbo(1, 1) - 1
2404 x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2405 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2406 ii = sent_size(ip)
2407 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
2408 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
2409 coarse_coeffs_pw%array(x_att, j, k) = coarse_coeffs_pw%array(x_att, j, k) + rcv_buf(ii)
2410 ii = ii + 1
2411 END DO
2412 END DO
2413 sent_size(ip) = ii
2414 END IF
2415 END DO
2416
2417 ii = sent_size(ip)
2418 DO x_att = max(pp_lb(ip), my_coarse_bo(1, 1)), min(pp_ub(ip), my_coarse_bo(2, 1))
2419 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
2420 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
2421 coarse_coeffs_pw%array(x_att, j, k) = coarse_coeffs_pw%array(x_att, j, k) + rcv_buf(ii)
2422 ii = ii + 1
2423 END DO
2424 END DO
2425 END DO
2426 sent_size(ip) = ii
2427
2428 DO x = coarse_gbo(2, 1) + 1, pp_ub(ip)
2429 x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2430 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2431 ii = sent_size(ip)
2432 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
2433 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
2434 coarse_coeffs_pw%array(x_att, j, k) = coarse_coeffs_pw%array(x_att, j, k) + rcv_buf(ii)
2435 ii = ii + 1
2436 END DO
2437 END DO
2438 sent_size(ip) = ii
2439 END IF
2440 END DO
2441
2442 END DO
2443
2444 IF (any(sent_size(0:n_procs - 2) /= rcv_offset(1:n_procs - 1))) &
2445 cpabort("error 1 handling the rcv buffer")
2446 IF (sent_size(n_procs - 1) /= rcv_tot_size) &
2447 cpabort("error 2 handling the rcv buffer")
2448
2449 ! dealloc
2450 DEALLOCATE (send_size, send_offset, rcv_size, rcv_offset)
2451 DEALLOCATE (send_buf, rcv_buf, real_rcv_size)
2452 DEALLOCATE (pp_ub, pp_lb)
2453 CALL timestop(handle2)
2454 ELSE
2455 cpassert(.NOT. local_data)
2456 END IF
2457
2458 CALL timestop(handle)
2459 END SUBROUTINE add_fine2coarse
2460
2461! **************************************************************************************************
2462!> \brief ...
2463!> \param preconditioner the preconditioner to create
2464!> \param precond_kind the kind of preconditioner to use
2465!> \param pool a pool with grids of the same type as the elements to
2466!> precondition
2467!> \param pbc if periodic boundary conditions should be applied
2468!> \param transpose ...
2469!> \author fawzi
2470! **************************************************************************************************
2471 SUBROUTINE pw_spline_precond_create(preconditioner, precond_kind, &
2472 pool, pbc, transpose)
2473 TYPE(pw_spline_precond_type), INTENT(OUT) :: preconditioner
2474 INTEGER, INTENT(in) :: precond_kind
2475 TYPE(pw_pool_type), INTENT(IN), POINTER :: pool
2476 LOGICAL, INTENT(in) :: pbc, transpose
2477
2478 preconditioner%kind = no_precond
2479 preconditioner%pool => pool
2480 preconditioner%pbc = pbc
2481 preconditioner%transpose = transpose
2482 CALL pool%retain()
2483 CALL pw_spline_precond_set_kind(preconditioner, precond_kind)
2484 END SUBROUTINE pw_spline_precond_create
2485
2486! **************************************************************************************************
2487!> \brief switches the types of precoditioner to use
2488!> \param preconditioner the preconditioner to be changed
2489!> \param precond_kind the new kind of preconditioner to use
2490!> \param pbc ...
2491!> \param transpose ...
2492!> \author fawzi
2493! **************************************************************************************************
2494 SUBROUTINE pw_spline_precond_set_kind(preconditioner, precond_kind, pbc, &
2495 transpose)
2496 TYPE(pw_spline_precond_type), INTENT(INOUT) :: preconditioner
2497 INTEGER, INTENT(in) :: precond_kind
2498 LOGICAL, INTENT(in), OPTIONAL :: pbc, transpose
2499
2500 LOGICAL :: do_3d_coeff
2501 REAL(kind=dp) :: s
2502
2503 IF (PRESENT(transpose)) preconditioner%transpose = transpose
2504 do_3d_coeff = .false.
2505 preconditioner%kind = precond_kind
2506 IF (PRESENT(pbc)) preconditioner%pbc = pbc
2507 SELECT CASE (precond_kind)
2508 CASE (no_precond)
2509 CASE (precond_spl3_aint2)
2510 preconditioner%coeffs_1d = (/-1.66_dp*0.25_dp, 1.66_dp, -1.66_dp*0.25_dp/)
2511 preconditioner%sharpen = .false.
2512 preconditioner%normalize = .false.
2513 do_3d_coeff = .true.
2514 CASE (precond_spl3_3)
2515 preconditioner%coeffs_1d(1) = -0.25_dp*1.6_dp
2516 preconditioner%coeffs_1d(2) = 1.6_dp
2517 preconditioner%coeffs_1d(3) = -0.25_dp*1.6_dp
2518 preconditioner%sharpen = .false.
2519 preconditioner%normalize = .false.
2520 do_3d_coeff = .true.
2521 CASE (precond_spl3_2)
2522 preconditioner%coeffs_1d(1) = -0.26_dp*1.76_dp
2523 preconditioner%coeffs_1d(2) = 1.76_dp
2524 preconditioner%coeffs_1d(3) = -0.26_dp*1.76_dp
2525 preconditioner%sharpen = .false.
2526 preconditioner%normalize = .false.
2527 do_3d_coeff = .true.
2528 CASE (precond_spl3_aint)
2529 preconditioner%coeffs_1d = spl3_1d_coeffs0
2530 preconditioner%sharpen = .true.
2531 preconditioner%normalize = .true.
2532 do_3d_coeff = .true.
2533 CASE (precond_spl3_1)
2534 preconditioner%coeffs_1d(1) = 0.5_dp/3._dp**(1._dp/3._dp)
2535 preconditioner%coeffs_1d(2) = 4._dp/3._dp**(1._dp/3._dp)
2536 preconditioner%coeffs_1d(3) = 0.5_dp/3._dp**(1._dp/3._dp)
2537 preconditioner%sharpen = .true.
2538 preconditioner%normalize = .false.
2539 do_3d_coeff = .true.
2540 CASE default
2541 cpabort("")
2542 END SELECT
2543 IF (do_3d_coeff) THEN
2544 s = 1._dp
2545 IF (preconditioner%sharpen) s = -1._dp
2546 preconditioner%coeffs(1) = &
2547 s*preconditioner%coeffs_1d(2)* &
2548 preconditioner%coeffs_1d(2)* &
2549 preconditioner%coeffs_1d(2)
2550 preconditioner%coeffs(2) = &
2551 s*preconditioner%coeffs_1d(1)* &
2552 preconditioner%coeffs_1d(2)* &
2553 preconditioner%coeffs_1d(2)
2554 preconditioner%coeffs(3) = &
2555 s*preconditioner%coeffs_1d(1)* &
2556 preconditioner%coeffs_1d(1)* &
2557 preconditioner%coeffs_1d(2)
2558 preconditioner%coeffs(4) = &
2559 s*preconditioner%coeffs_1d(1)* &
2560 preconditioner%coeffs_1d(1)* &
2561 preconditioner%coeffs_1d(1)
2562 IF (preconditioner%sharpen) THEN
2563 IF (preconditioner%normalize) THEN
2564 preconditioner%coeffs(1) = 2._dp + &
2565 preconditioner%coeffs(1)
2566 ELSE
2567 preconditioner%coeffs(1) = -preconditioner%coeffs(1)
2568 END IF
2569 END IF
2570 END IF
2571 END SUBROUTINE pw_spline_precond_set_kind
2572
2573! **************************************************************************************************
2574!> \brief releases the preconditioner
2575!> \param preconditioner the preconditioner to release
2576!> \author fawzi
2577! **************************************************************************************************
2578 SUBROUTINE pw_spline_precond_release(preconditioner)
2579 TYPE(pw_spline_precond_type), INTENT(INOUT) :: preconditioner
2580
2581 CALL pw_pool_release(preconditioner%pool)
2582 END SUBROUTINE pw_spline_precond_release
2583
2584! **************************************************************************************************
2585!> \brief applies the preconditioner to the system of equations to find the
2586!> coefficients of the spline
2587!> \param preconditioner the preconditioner to apply
2588!> \param in_v the grid on which the preconditioner should be applied
2589!> \param out_v place to store the preconditioner applied on v_out
2590!> \author fawzi
2591! **************************************************************************************************
2592 SUBROUTINE pw_spline_do_precond(preconditioner, in_v, out_v)
2593 TYPE(pw_spline_precond_type), INTENT(IN) :: preconditioner
2594 TYPE(pw_r3d_rs_type), INTENT(IN) :: in_v
2595 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: out_v
2596
2597 SELECT CASE (preconditioner%kind)
2598 CASE (no_precond)
2599 CALL pw_copy(in_v, out_v)
2600 CASE (precond_spl3_aint, precond_spl3_1)
2601 CALL pw_zero(out_v)
2602 IF (preconditioner%pbc) THEN
2603 CALL pw_nn_smear_r(pw_in=in_v, pw_out=out_v, &
2604 coeffs=preconditioner%coeffs)
2605 ELSE
2606 CALL pw_nn_compose_r_no_pbc(weights_1d=preconditioner%coeffs_1d, &
2607 pw_in=in_v, pw_out=out_v, sharpen=preconditioner%sharpen, &
2608 normalize=preconditioner%normalize, &
2609 transpose=preconditioner%transpose)
2610 END IF
2611 CASE (precond_spl3_3, precond_spl3_2, precond_spl3_aint2)
2612 CALL pw_zero(out_v)
2613 IF (preconditioner%pbc) THEN
2614 CALL pw_nn_smear_r(pw_in=in_v, pw_out=out_v, &
2615 coeffs=preconditioner%coeffs)
2616 ELSE
2617 CALL pw_nn_compose_r_no_pbc(weights_1d=preconditioner%coeffs_1d, &
2618 pw_in=in_v, pw_out=out_v, sharpen=preconditioner%sharpen, &
2619 normalize=preconditioner%normalize, smooth_boundary=.true., &
2620 transpose=preconditioner%transpose)
2621 END IF
2622 CASE default
2623 cpabort("")
2624 END SELECT
2625 END SUBROUTINE pw_spline_do_precond
2626
2627! **************************************************************************************************
2628!> \brief solves iteratively (CG) a systmes of linear equations
2629!> linOp(coeffs)=values
2630!> (for example those needed to find the coefficients of a spline)
2631!> Returns true if the it succeeded to achieve the requested accuracy
2632!> \param values the right hand side of the system
2633!> \param coeffs will contain the solution of the system (and on entry
2634!> it contains the starting point)
2635!> \param linOp the linear operator to be inverted
2636!> \param preconditioner the preconditioner to apply
2637!> \param pool a pool of grids (for the temporary objects)
2638!> \param eps_r the requested precision on the residual
2639!> \param eps_x the requested precision on the solution
2640!> \param max_iter maximum number of iteration allowed
2641!> \param sumtype ...
2642!> \return ...
2643!> \author fawzi
2644! **************************************************************************************************
2645 FUNCTION find_coeffs(values, coeffs, linOp, preconditioner, pool, &
2646 eps_r, eps_x, max_iter, sumtype) RESULT(res)
2647
2648 TYPE(pw_r3d_rs_type), INTENT(IN) :: values
2649 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: coeffs
2650 INTERFACE
2651 SUBROUTINE linop(pw_in, pw_out)
2652 USE pw_types, ONLY: pw_r3d_rs_type
2653 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
2654 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
2655 END SUBROUTINE linop
2656 END INTERFACE
2657 TYPE(pw_spline_precond_type), INTENT(IN) :: preconditioner
2658 TYPE(pw_pool_type), POINTER :: pool
2659 REAL(kind=dp), INTENT(in) :: eps_r, eps_x
2660 INTEGER, INTENT(in) :: max_iter
2661 INTEGER, INTENT(in), OPTIONAL :: sumtype
2662 LOGICAL :: res
2663
2664 INTEGER :: i, iiter, iter, j, k
2665 INTEGER, DIMENSION(2, 3) :: bo
2666 LOGICAL :: last
2667 REAL(kind=dp) :: alpha, beta, eps_r_att, eps_x_att, r_z, &
2668 r_z_new
2669 TYPE(cp_logger_type), POINTER :: logger
2670 TYPE(pw_r3d_rs_type) :: ap, p, r, z
2671
2672 last = .false.
2673
2674 res = .false.
2675 logger => cp_get_default_logger()
2676 CALL pool%create_pw(r)
2677 CALL pool%create_pw(z)
2678 CALL pool%create_pw(p)
2679 CALL pool%create_pw(ap)
2680
2681 !CALL cp_add_iter_level(logger%iter_info,level_name="SPLINE_FIND_COEFFS")
2682 ext_do: DO iiter = 1, max_iter, 10
2683 CALL pw_zero(r)
2684 CALL linop(pw_in=coeffs, pw_out=r)
2685 r%array = -r%array
2686 CALL pw_axpy(values, r)
2687 CALL pw_spline_do_precond(preconditioner, in_v=r, out_v=z)
2688 CALL pw_copy(z, p)
2689 r_z = pw_integral_ab(r, z, sumtype)
2690
2691 DO iter = iiter, min(iiter + 9, max_iter)
2692 eps_r_att = sqrt(pw_integral_ab(r, r, sumtype))
2693 IF (eps_r_att == 0._dp) THEN
2694 eps_x_att = 0._dp
2695 last = .true.
2696 ELSE
2697 CALL pw_zero(ap)
2698 CALL linop(pw_in=p, pw_out=ap)
2699 alpha = r_z/pw_integral_ab(ap, p, sumtype)
2700
2701 CALL pw_axpy(p, coeffs, alpha=alpha)
2702
2703 eps_x_att = alpha*sqrt(pw_integral_ab(p, p, sumtype)) ! try to spare if unneeded?
2704 IF (eps_r_att < eps_r .AND. eps_x_att < eps_x) last = .true.
2705 END IF
2706 !CALL cp_iterate(logger%iter_info,last=last)
2707 IF (last) THEN
2708 res = .true.
2709 EXIT ext_do
2710 END IF
2711
2712 CALL pw_axpy(ap, r, alpha=-alpha)
2713
2714 CALL pw_spline_do_precond(preconditioner, in_v=r, out_v=z)
2715
2716 r_z_new = pw_integral_ab(r, z, sumtype)
2717 beta = r_z_new/r_z
2718 r_z = r_z_new
2719
2720 bo = p%pw_grid%bounds_local
2721 DO k = bo(1, 3), bo(2, 3)
2722 DO j = bo(1, 2), bo(2, 2)
2723 DO i = bo(1, 1), bo(2, 1)
2724 p%array(i, j, k) = z%array(i, j, k) + beta*p%array(i, j, k)
2725 END DO
2726 END DO
2727 END DO
2728
2729 END DO
2730 END DO ext_do
2731 !CALL cp_rm_iter_level(logger%iter_info,level_name="SPLINE_FIND_COEFFS")
2732
2733 CALL pool%give_back_pw(r)
2734 CALL pool%give_back_pw(z)
2735 CALL pool%give_back_pw(p)
2736 CALL pool%give_back_pw(ap)
2737
2738 END FUNCTION find_coeffs
2739
2740! **************************************************************************************************
2741!> \brief adds to pw_out pw_in composed with the weights
2742!> pw_out%array(i,j,k)=pw_out%array(i,j,k)+sum(pw_in%array(i+l,j+m,k+n)*
2743!> weights_1d(abs(l)+1)*weights_1d(abs(m)+1)*weights_1d(abs(n)+1),
2744!> l=-1..1,m=-1..1,n=-1..1)
2745!> \param weights_1d ...
2746!> \param pw_in ...
2747!> \param pw_out ...
2748!> \param sharpen ...
2749!> \param normalize ...
2750!> \param transpose ...
2751!> \param smooth_boundary ...
2752!> \author fawzi
2753! **************************************************************************************************
2754 SUBROUTINE pw_nn_compose_r_no_pbc(weights_1d, pw_in, pw_out, &
2755 sharpen, normalize, transpose, smooth_boundary)
2756 REAL(kind=dp), DIMENSION(-1:1) :: weights_1d
2757 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in, pw_out
2758 LOGICAL, INTENT(in), OPTIONAL :: sharpen, normalize, transpose, &
2759 smooth_boundary
2760
2761 INTEGER :: first_index, i, j, jw, k, kw, &
2762 last_index, myj, myk, n_els
2763 INTEGER, DIMENSION(2, 3) :: bo, gbo
2764 INTEGER, DIMENSION(3) :: s
2765 LOGICAL :: has_l_boundary, has_u_boundary, &
2766 is_split, my_normalize, my_sharpen, &
2767 my_smooth_boundary, my_transpose
2768 REAL(kind=dp) :: in_val_f, in_val_l, in_val_tmp, w_j, w_k
2769 REAL(kind=dp), DIMENSION(-1:1) :: w
2770 REAL(kind=dp), DIMENSION(:, :), POINTER :: l_boundary, tmp, u_boundary
2771 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: in_val, out_val
2772
2773 bo = pw_in%pw_grid%bounds_local
2774 gbo = pw_in%pw_grid%bounds
2775 in_val => pw_in%array
2776 out_val => pw_out%array
2777 my_sharpen = .false.
2778 IF (PRESENT(sharpen)) my_sharpen = sharpen
2779 my_normalize = .false.
2780 IF (PRESENT(normalize)) my_normalize = normalize
2781 my_transpose = .false.
2782 IF (PRESENT(transpose)) my_transpose = transpose
2783 my_smooth_boundary = .false.
2784 IF (PRESENT(smooth_boundary)) my_smooth_boundary = smooth_boundary
2785 cpassert(.NOT. my_normalize .OR. my_sharpen)
2786 cpassert(.NOT. my_smooth_boundary .OR. .NOT. my_sharpen)
2787 DO i = 1, 3
2788 s(i) = bo(2, i) - bo(1, i) + 1
2789 END DO
2790 IF (any(s < 1)) RETURN
2791 is_split = any(pw_in%pw_grid%bounds_local(:, 1) /= &
2792 pw_in%pw_grid%bounds(:, 1))
2793 has_l_boundary = (gbo(1, 1) == bo(1, 1))
2794 has_u_boundary = (gbo(2, 1) == bo(2, 1))
2795 IF (is_split) THEN
2796 ALLOCATE (l_boundary(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)), &
2797 u_boundary(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)), &
2798 tmp(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
2799 tmp(:, :) = pw_in%array(bo(2, 1), :, :)
2800 CALL pw_in%pw_grid%para%group%sendrecv(tmp, pw_in%pw_grid%para%pos_of_x( &
2801 gbo(1, 1) + modulo(bo(2, 1) + 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)), &
2802 l_boundary, pw_in%pw_grid%para%pos_of_x( &
2803 gbo(1, 1) + modulo(bo(1, 1) - 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)))
2804 tmp(:, :) = pw_in%array(bo(1, 1), :, :)
2805 CALL pw_in%pw_grid%para%group%sendrecv(tmp, pw_in%pw_grid%para%pos_of_x( &
2806 gbo(1, 1) + modulo(bo(1, 1) - 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)), &
2807 u_boundary, pw_in%pw_grid%para%pos_of_x( &
2808 gbo(1, 1) + modulo(bo(2, 1) + 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)))
2809 DEALLOCATE (tmp)
2810 END IF
2811
2812 n_els = s(1)
2813 IF (has_l_boundary) THEN
2814 n_els = n_els - 1
2815 first_index = bo(1, 1) + 1
2816 ELSE
2817 first_index = bo(1, 1)
2818 END IF
2819 IF (has_u_boundary) THEN
2820 n_els = n_els - 1
2821 last_index = bo(2, 1) - 1
2822 ELSE
2823 last_index = bo(2, 1)
2824 END IF
2825!$OMP PARALLEL DO DEFAULT(NONE) &
2826!$OMP PRIVATE(k, kw, myk, j, jw, myj, in_val_f, in_val_l, w_k, w_j, in_val_tmp, w) &
2827!$OMP SHARED(bo, in_val, out_val, s, l_boundary, u_boundary, weights_1d, is_split, &
2828!$OMP my_transpose, gbo, my_smooth_boundary, has_l_boundary, has_u_boundary, &
2829!$OMP my_sharpen, last_index, first_index, my_normalize, n_els)
2830 DO k = bo(1, 3), bo(2, 3)
2831 DO kw = -1, 1
2832 myk = k + kw
2833 IF (my_transpose) THEN
2834 IF (k >= gbo(2, 3) - 1 .OR. k <= gbo(1, 3) + 1) THEN
2835 IF (k == gbo(2, 3) .OR. k == gbo(1, 3)) THEN
2836 IF (myk < gbo(2, 3) .AND. myk > gbo(1, 3)) THEN
2837 w_k = weights_1d(kw)
2838 IF (my_smooth_boundary) THEN
2839 w_k = weights_1d(kw)/weights_1d(0)
2840 END IF
2841 ELSE IF (kw == 0) THEN
2842 w_k = 1._dp
2843 ELSE
2844 cycle
2845 END IF
2846 ELSE
2847 IF (myk == gbo(2, 3) .OR. myk == gbo(1, 3)) cycle
2848 w_k = weights_1d(kw)
2849 END IF
2850 ELSE
2851 w_k = weights_1d(kw)
2852 END IF
2853 ELSE
2854 IF (k >= gbo(2, 3) - 1 .OR. k <= gbo(1, 3) + 1) THEN
2855 IF (k == gbo(2, 3) .OR. k == gbo(1, 3)) THEN
2856 IF (kw /= 0) cycle
2857 w_k = 1._dp
2858 ELSE
2859 IF (my_smooth_boundary .AND. ((k == gbo(1, 3) + 1 .AND. myk == gbo(1, 3)) .OR. &
2860 (k == gbo(2, 3) - 1 .AND. myk == gbo(2, 3)))) THEN
2861 w_k = weights_1d(kw)/weights_1d(0)
2862 ELSE
2863 w_k = weights_1d(kw)
2864 END IF
2865 END IF
2866 ELSE
2867 w_k = weights_1d(kw)
2868 END IF
2869 END IF
2870 DO j = bo(1, 2), bo(2, 2)
2871 DO jw = -1, 1
2872 myj = j + jw
2873 IF (j < gbo(2, 2) - 1 .AND. j > gbo(1, 2) + 1) THEN
2874 w_j = w_k*weights_1d(jw)
2875 ELSE
2876 IF (my_transpose) THEN
2877 IF (j == gbo(2, 2) .OR. j == gbo(1, 2)) THEN
2878 IF (myj < gbo(2, 2) .AND. myj > gbo(1, 2)) THEN
2879 w_j = weights_1d(jw)*w_k
2880 IF (my_smooth_boundary) THEN
2881 w_j = weights_1d(jw)/weights_1d(0)*w_k
2882 END IF
2883 ELSE IF (jw == 0) THEN
2884 w_j = w_k
2885 ELSE
2886 cycle
2887 END IF
2888 ELSE
2889 IF (myj == gbo(2, 2) .OR. myj == gbo(1, 2)) cycle
2890 w_j = w_k*weights_1d(jw)
2891 END IF
2892 ELSE
2893 IF (j == gbo(2, 2) .OR. j == gbo(1, 2)) THEN
2894 IF (jw /= 0) cycle
2895 w_j = w_k
2896 ELSE IF (my_smooth_boundary .AND. ((j == gbo(1, 2) + 1 .AND. myj == gbo(1, 2)) .OR. &
2897 (j == gbo(2, 2) - 1 .AND. myj == gbo(2, 2)))) THEN
2898 w_j = w_k*weights_1d(jw)/weights_1d(0)
2899 ELSE
2900 w_j = w_k*weights_1d(jw)
2901 END IF
2902 END IF
2903 END IF
2904
2905 IF (has_l_boundary) THEN
2906 IF (my_transpose) THEN
2907 IF (s(1) == 1) THEN
2908 cpassert(.NOT. has_u_boundary)
2909 in_val_tmp = u_boundary(myj, myk)
2910 ELSE
2911 in_val_tmp = in_val(bo(1, 1) + 1, myj, myk)
2912 END IF
2913 IF (my_sharpen) THEN
2914 IF (kw == 0 .AND. jw == 0) THEN
2915 IF (my_normalize) THEN
2916 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2917 (2.0_dp - w_j)*in_val(bo(1, 1), myj, myk) - &
2918 in_val_tmp*weights_1d(1)*w_j
2919 ELSE
2920 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2921 in_val(bo(1, 1), myj, myk)*w_j - &
2922 in_val_tmp*weights_1d(1)*w_j
2923 END IF
2924 ELSE
2925 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) - &
2926 in_val(bo(1, 1), myj, myk)*w_j - &
2927 in_val_tmp*weights_1d(1)*w_j
2928 END IF
2929 ELSE IF (my_smooth_boundary) THEN
2930 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2931 w_j*(in_val(bo(1, 1), myj, myk) + &
2932 in_val_tmp*weights_1d(1)/weights_1d(0))
2933 ELSE
2934 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2935 w_j*(in_val(bo(1, 1), myj, myk) + &
2936 in_val_tmp*weights_1d(1))
2937 END IF
2938 in_val_f = 0.0_dp
2939 ELSE
2940 in_val_f = in_val(bo(1, 1), myj, myk)
2941 IF (my_sharpen) THEN
2942 IF (kw == 0 .AND. jw == 0) THEN
2943 IF (my_normalize) THEN
2944 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2945 (2.0_dp - w_j)*in_val_f
2946 ELSE
2947 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2948 in_val_f*w_j
2949 END IF
2950 ELSE
2951 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) - &
2952 in_val_f*w_j
2953 END IF
2954 ELSE
2955 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2956 in_val_f*w_j
2957 END IF
2958 END IF
2959 ELSE
2960 in_val_f = l_boundary(myj, myk)
2961 END IF
2962 IF (has_u_boundary) THEN
2963 IF (my_transpose) THEN
2964 in_val_l = in_val(bo(2, 1), myj, myk)
2965 IF (s(1) == 1) THEN
2966 cpassert(.NOT. has_l_boundary)
2967 in_val_tmp = l_boundary(myj, myk)
2968 ELSE
2969 in_val_tmp = in_val(bo(2, 1) - 1, myj, myk)
2970 END IF
2971 IF (my_sharpen) THEN
2972 IF (kw == 0 .AND. jw == 0) THEN
2973 IF (my_normalize) THEN
2974 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2975 in_val_l*(2._dp - w_j) - &
2976 in_val_tmp*weights_1d(1)*w_j
2977 ELSE
2978 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2979 in_val_l*w_j - &
2980 in_val_tmp*weights_1d(1)*w_j
2981 END IF
2982 ELSE
2983 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) - &
2984 w_j*in_val_l - &
2985 in_val_tmp*weights_1d(1)*w_j
2986 END IF
2987 ELSE IF (my_smooth_boundary) THEN
2988 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2989 w_j*(in_val_l + in_val_tmp*weights_1d(1)/weights_1d(0))
2990 ELSE
2991 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2992 w_j*(in_val_l + in_val_tmp*weights_1d(1))
2993 END IF
2994 in_val_l = 0._dp
2995 ELSE
2996 in_val_l = in_val(bo(2, 1), myj, myk)
2997 IF (my_sharpen) THEN
2998 IF (kw == 0 .AND. jw == 0) THEN
2999 IF (my_normalize) THEN
3000 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
3001 in_val_l*(2._dp - w_j)
3002 ELSE
3003 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
3004 in_val_l*w_j
3005 END IF
3006 ELSE
3007 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) - &
3008 w_j*in_val_l
3009 END IF
3010 ELSE
3011 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
3012 w_j*in_val_l
3013 END IF
3014 END IF
3015 ELSE
3016 in_val_l = u_boundary(myj, myk)
3017 END IF
3018 IF (last_index >= first_index) THEN
3019 IF (my_transpose) THEN
3020 IF (bo(1, 1) - 1 == gbo(1, 1)) THEN
3021 in_val_f = 0._dp
3022 ELSE IF (bo(2, 1) + 1 == gbo(2, 1)) THEN
3023 in_val_l = 0._dp
3024 END IF
3025 END IF
3026 IF (my_sharpen) THEN
3027 w = -weights_1d*w_j
3028 IF (kw == 0 .AND. jw == 0) THEN
3029 IF (my_normalize) THEN
3030 w(0) = w(0) + 2._dp
3031 ELSE
3032 w(0) = -w(0)
3033 END IF
3034 END IF
3035 ELSE
3036 w = weights_1d*w_j
3037 END IF
3038 IF (my_smooth_boundary .AND. (.NOT. my_transpose)) THEN
3039 IF (gbo(1, 1) + 1 >= bo(1, 1) .AND. &
3040 gbo(1, 1) + 1 <= bo(2, 1) .AND. gbo(2, 1) - gbo(1, 1) > 2) THEN
3041 IF (gbo(1, 1) >= bo(1, 1)) THEN
3042 out_val(gbo(1, 1) + 1, j, k) = out_val(gbo(1, 1) + 1, j, k) + &
3043 in_val(gbo(1, 1), myj, myk)*w_j*weights_1d(-1)* &
3044 (1._dp/weights_1d(0) - 1._dp)
3045 ELSE
3046 out_val(gbo(1, 1) + 1, j, k) = out_val(gbo(1, 1) + 1, j, k) + &
3047 l_boundary(myj, myk)*w_j*weights_1d(-1)* &
3048 (1._dp/weights_1d(0) - 1._dp)
3049 END IF
3050 END IF
3051 END IF
3052 CALL pw_compose_stripe(weights=w, &
3053 in_val=in_val(first_index:last_index, myj, myk), &
3054 in_val_first=in_val_f, in_val_last=in_val_l, &
3055 out_val=out_val(first_index:last_index, j, k), &
3056 n_el=n_els)
3057!FM call pw_compose_stripe2(weights=w,&
3058!FM in_val=in_val,&
3059!FM in_val_first=in_val_f,in_val_last=in_val_l,&
3060!FM out_val=out_val,&
3061!FM first_val=first_index,last_val=last_index,&
3062!FM myj=myj,myk=myk,j=j,k=k)
3063 IF (my_smooth_boundary .AND. (.NOT. my_transpose)) THEN
3064 IF (gbo(2, 1) - 1 >= bo(1, 1) .AND. &
3065 gbo(2, 1) - 1 <= bo(2, 1) .AND. gbo(2, 1) - gbo(1, 1) > 2) THEN
3066 IF (gbo(2, 1) <= bo(2, 1)) THEN
3067 out_val(gbo(2, 1) - 1, j, k) = out_val(gbo(2, 1) - 1, j, k) + &
3068 in_val(gbo(2, 1), myj, myk)*w_j*weights_1d(1)* &
3069 (1._dp/weights_1d(0) - 1._dp)
3070 ELSE
3071 out_val(gbo(2, 1) - 1, j, k) = out_val(gbo(2, 1) - 1, j, k) + &
3072 u_boundary(myj, myk)*w_j*weights_1d(1)* &
3073 (1._dp/weights_1d(0) - 1._dp)
3074 END IF
3075 END IF
3076 END IF
3077
3078 END IF
3079 END DO
3080 END DO
3081 END DO
3082 END DO
3083
3084 IF (is_split) THEN
3085 DEALLOCATE (l_boundary, u_boundary)
3086 END IF
3087 END SUBROUTINE pw_nn_compose_r_no_pbc
3088
3089! **************************************************************************************************
3090!> \brief ...
3091!> \param pw_in ...
3092!> \param pw_out ...
3093! **************************************************************************************************
3094 SUBROUTINE spl3_nopbc(pw_in, pw_out)
3095
3096 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
3097 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
3098
3099 CALL pw_zero(pw_out)
3100 CALL pw_nn_compose_r_no_pbc(weights_1d=spl3_1d_coeffs0, pw_in=pw_in, &
3101 pw_out=pw_out, sharpen=.false., normalize=.false.)
3102
3103 END SUBROUTINE spl3_nopbc
3104
3105! **************************************************************************************************
3106!> \brief ...
3107!> \param pw_in ...
3108!> \param pw_out ...
3109! **************************************************************************************************
3110 SUBROUTINE spl3_nopbct(pw_in, pw_out)
3111
3112 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
3113 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
3114
3115 CALL pw_zero(pw_out)
3116 CALL pw_nn_compose_r_no_pbc(weights_1d=spl3_1d_coeffs0, pw_in=pw_in, &
3117 pw_out=pw_out, sharpen=.false., normalize=.false., transpose=.true.)
3118
3119 END SUBROUTINE spl3_nopbct
3120
3121! **************************************************************************************************
3122!> \brief ...
3123!> \param pw_in ...
3124!> \param pw_out ...
3125! **************************************************************************************************
3126 SUBROUTINE spl3_pbc(pw_in, pw_out)
3127
3128 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
3129 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
3130
3131 CALL pw_zero(pw_out)
3132 CALL pw_nn_smear_r(pw_in, pw_out, coeffs=spline3_coeffs)
3133
3134 END SUBROUTINE spl3_pbc
3135
3136! **************************************************************************************************
3137!> \brief Evaluates the PBC interpolated Spline (pw) function on the generic
3138!> input vector (vec)
3139!> \param vec ...
3140!> \param pw ...
3141!> \return ...
3142!> \par History
3143!> 12.2007 Adapted for use with distributed grids [rdeclerck]
3144!> \author Teodoro Laino 12/2005 [tlaino]
3145!> \note
3146!> Requires the Spline coefficients to be computed with PBC
3147! **************************************************************************************************
3148 FUNCTION eval_interp_spl3_pbc(vec, pw) RESULT(val)
3149 REAL(kind=dp), DIMENSION(3), INTENT(in) :: vec
3150 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
3151 REAL(kind=dp) :: val
3152
3153 INTEGER :: i, ivec(3), j, k, npts(3)
3154 INTEGER, DIMENSION(2, 3) :: bo, bo_l
3155 INTEGER, DIMENSION(4) :: ii, ij, ik
3156 LOGICAL :: my_mpsum
3157 REAL(kind=dp) :: a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, dr1, dr2, dr3, e1, e2, e3, &
3158 f1, f2, f3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, r1, r2, r3, s1, s2, s3, s4, &
3159 t1, t2, t3, t4, u1, u2, u3, v1, v2, v3, v4, xd1, xd2, xd3
3160 REAL(kind=dp), DIMENSION(4, 4, 4) :: box
3161 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
3162
3163 NULLIFY (grid)
3164 my_mpsum = (pw%pw_grid%para%mode /= pw_mode_local)
3165 npts = pw%pw_grid%npts
3166 ivec = floor(vec/pw%pw_grid%dr)
3167 dr1 = pw%pw_grid%dr(1)
3168 dr2 = pw%pw_grid%dr(2)
3169 dr3 = pw%pw_grid%dr(3)
3170
3171 xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
3172 xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
3173 xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
3174 grid => pw%array(:, :, :)
3175 bo = pw%pw_grid%bounds
3176 bo_l = pw%pw_grid%bounds_local
3177
3178 ik(1) = modulo(ivec(3) - 1, npts(3)) + bo(1, 3)
3179 ik(2) = modulo(ivec(3), npts(3)) + bo(1, 3)
3180 ik(3) = modulo(ivec(3) + 1, npts(3)) + bo(1, 3)
3181 ik(4) = modulo(ivec(3) + 2, npts(3)) + bo(1, 3)
3182
3183 ij(1) = modulo(ivec(2) - 1, npts(2)) + bo(1, 2)
3184 ij(2) = modulo(ivec(2), npts(2)) + bo(1, 2)
3185 ij(3) = modulo(ivec(2) + 1, npts(2)) + bo(1, 2)
3186 ij(4) = modulo(ivec(2) + 2, npts(2)) + bo(1, 2)
3187
3188 ii(1) = modulo(ivec(1) - 1, npts(1)) + bo(1, 1)
3189 ii(2) = modulo(ivec(1), npts(1)) + bo(1, 1)
3190 ii(3) = modulo(ivec(1) + 1, npts(1)) + bo(1, 1)
3191 ii(4) = modulo(ivec(1) + 2, npts(1)) + bo(1, 1)
3192
3193 DO k = 1, 4
3194 DO j = 1, 4
3195 DO i = 1, 4
3196 IF ( &
3197 ii(i) >= bo_l(1, 1) .AND. &
3198 ii(i) <= bo_l(2, 1) .AND. &
3199 ij(j) >= bo_l(1, 2) .AND. &
3200 ij(j) <= bo_l(2, 2) .AND. &
3201 ik(k) >= bo_l(1, 3) .AND. &
3202 ik(k) <= bo_l(2, 3) &
3203 ) THEN
3204 box(i, j, k) = grid(ii(i) + 1 - bo_l(1, 1), &
3205 ij(j) + 1 - bo_l(1, 2), &
3206 ik(k) + 1 - bo_l(1, 3))
3207 ELSE
3208 box(i, j, k) = 0.0_dp
3209 END IF
3210 END DO
3211 END DO
3212 END DO
3213
3214 a1 = 3.0_dp + xd1
3215 a2 = a1*a1
3216 a3 = a2*a1
3217 b1 = 2.0_dp + xd1
3218 b2 = b1*b1
3219 b3 = b2*b1
3220 c1 = 1.0_dp + xd1
3221 c2 = c1*c1
3222 c3 = c2*c1
3223 d1 = xd1
3224 d2 = d1*d1
3225 d3 = d2*d1
3226 e1 = 3.0_dp + xd2
3227 e2 = e1*e1
3228 e3 = e2*e1
3229 f1 = 2.0_dp + xd2
3230 f2 = f1*f1
3231 f3 = f2*f1
3232 g1 = 1.0_dp + xd2
3233 g2 = g1*g1
3234 g3 = g2*g1
3235 h1 = xd2
3236 h2 = h1*h1
3237 h3 = h2*h1
3238 p1 = 3.0_dp + xd3
3239 p2 = p1*p1
3240 p3 = p2*p1
3241 q1 = 2.0_dp + xd3
3242 q2 = q1*q1
3243 q3 = q2*q1
3244 r1 = 1.0_dp + xd3
3245 r2 = r1*r1
3246 r3 = r2*r1
3247 u1 = xd3
3248 u2 = u1*u1
3249 u3 = u2*u1
3250
3251 t1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
3252 t2 = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
3253 t3 = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
3254 t4 = 1.0_dp/6.0_dp*d3
3255 s1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
3256 s2 = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
3257 s3 = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
3258 s4 = 1.0_dp/6.0_dp*h3
3259 v1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
3260 v2 = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
3261 v3 = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
3262 v4 = 1.0_dp/6.0_dp*u3
3263
3264 val = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3265 (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3266 (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3267 (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3268 ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3269 (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3270 (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3271 (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3272 ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3273 (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3274 (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3275 (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3276 ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3277 (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3278 (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3279 (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3280
3281 IF (my_mpsum) CALL pw%pw_grid%para%group%sum(val)
3282
3283 END FUNCTION eval_interp_spl3_pbc
3284
3285! **************************************************************************************************
3286!> \brief Evaluates the derivatives of the PBC interpolated Spline (pw)
3287!> function on the generic input vector (vec)
3288!> \param vec ...
3289!> \param pw ...
3290!> \return ...
3291!> \par History
3292!> 12.2007 Adapted for use with distributed grids [rdeclerck]
3293!> \author Teodoro Laino 12/2005 [tlaino]
3294!> \note
3295!> Requires the Spline coefficients to be computed with PBC
3296! **************************************************************************************************
3297 FUNCTION eval_d_interp_spl3_pbc(vec, pw) RESULT(val)
3298 REAL(kind=dp), DIMENSION(3), INTENT(in) :: vec
3299 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
3300 REAL(kind=dp) :: val(3)
3301
3302 INTEGER :: i, ivec(3), j, k, npts(3)
3303 INTEGER, DIMENSION(2, 3) :: bo, bo_l
3304 INTEGER, DIMENSION(4) :: ii, ij, ik
3305 LOGICAL :: my_mpsum
3306 REAL(kind=dp) :: a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, dr1, dr1i, dr2, dr2i, dr3, &
3307 dr3i, e1, e2, e3, f1, f2, f3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, r1, r2, r3, &
3308 s1, s1d, s1o, s2, s2d, s2o, s3, s3d, s3o, s4, s4d, s4o, t1, t1d, t1o, t2, t2d, t2o, t3, &
3309 t3d, t3o, t4, t4d, t4o, u1, u2, u3, v1, v1d, v1o, v2, v2d, v2o, v3, v3d, v3o, v4, v4d, &
3310 v4o, xd1, xd2, xd3
3311 REAL(kind=dp), DIMENSION(4, 4, 4) :: box
3312 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
3313
3314 NULLIFY (grid)
3315 my_mpsum = (pw%pw_grid%para%mode /= pw_mode_local)
3316 npts = pw%pw_grid%npts
3317 ivec = floor(vec/pw%pw_grid%dr)
3318 dr1 = pw%pw_grid%dr(1)
3319 dr2 = pw%pw_grid%dr(2)
3320 dr3 = pw%pw_grid%dr(3)
3321 dr1i = 1.0_dp/dr1
3322 dr2i = 1.0_dp/dr2
3323 dr3i = 1.0_dp/dr3
3324 xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
3325 xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
3326 xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
3327 grid => pw%array(:, :, :)
3328 bo = pw%pw_grid%bounds
3329 bo_l = pw%pw_grid%bounds_local
3330
3331 ik(1) = modulo(ivec(3) - 1, npts(3)) + bo(1, 3)
3332 ik(2) = modulo(ivec(3), npts(3)) + bo(1, 3)
3333 ik(3) = modulo(ivec(3) + 1, npts(3)) + bo(1, 3)
3334 ik(4) = modulo(ivec(3) + 2, npts(3)) + bo(1, 3)
3335
3336 ij(1) = modulo(ivec(2) - 1, npts(2)) + bo(1, 2)
3337 ij(2) = modulo(ivec(2), npts(2)) + bo(1, 2)
3338 ij(3) = modulo(ivec(2) + 1, npts(2)) + bo(1, 2)
3339 ij(4) = modulo(ivec(2) + 2, npts(2)) + bo(1, 2)
3340
3341 ii(1) = modulo(ivec(1) - 1, npts(1)) + bo(1, 1)
3342 ii(2) = modulo(ivec(1), npts(1)) + bo(1, 1)
3343 ii(3) = modulo(ivec(1) + 1, npts(1)) + bo(1, 1)
3344 ii(4) = modulo(ivec(1) + 2, npts(1)) + bo(1, 1)
3345
3346 DO k = 1, 4
3347 DO j = 1, 4
3348 DO i = 1, 4
3349 IF ( &
3350 ii(i) >= bo_l(1, 1) .AND. &
3351 ii(i) <= bo_l(2, 1) .AND. &
3352 ij(j) >= bo_l(1, 2) .AND. &
3353 ij(j) <= bo_l(2, 2) .AND. &
3354 ik(k) >= bo_l(1, 3) .AND. &
3355 ik(k) <= bo_l(2, 3) &
3356 ) THEN
3357 box(i, j, k) = grid(ii(i) + 1 - bo_l(1, 1), &
3358 ij(j) + 1 - bo_l(1, 2), &
3359 ik(k) + 1 - bo_l(1, 3))
3360 ELSE
3361 box(i, j, k) = 0.0_dp
3362 END IF
3363 END DO
3364 END DO
3365 END DO
3366
3367 a1 = 3.0_dp + xd1
3368 a2 = a1*a1
3369 a3 = a2*a1
3370 b1 = 2.0_dp + xd1
3371 b2 = b1*b1
3372 b3 = b2*b1
3373 c1 = 1.0_dp + xd1
3374 c2 = c1*c1
3375 c3 = c2*c1
3376 d1 = xd1
3377 d2 = d1*d1
3378 d3 = d2*d1
3379 e1 = 3.0_dp + xd2
3380 e2 = e1*e1
3381 e3 = e2*e1
3382 f1 = 2.0_dp + xd2
3383 f2 = f1*f1
3384 f3 = f2*f1
3385 g1 = 1.0_dp + xd2
3386 g2 = g1*g1
3387 g3 = g2*g1
3388 h1 = xd2
3389 h2 = h1*h1
3390 h3 = h2*h1
3391 p1 = 3.0_dp + xd3
3392 p2 = p1*p1
3393 p3 = p2*p1
3394 q1 = 2.0_dp + xd3
3395 q2 = q1*q1
3396 q3 = q2*q1
3397 r1 = 1.0_dp + xd3
3398 r2 = r1*r1
3399 r3 = r2*r1
3400 u1 = xd3
3401 u2 = u1*u1
3402 u3 = u2*u1
3403
3404 t1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
3405 t2o = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
3406 t3o = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
3407 t4o = 1.0_dp/6.0_dp*d3
3408 s1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
3409 s2o = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
3410 s3o = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
3411 s4o = 1.0_dp/6.0_dp*h3
3412 v1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
3413 v2o = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
3414 v3o = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
3415 v4o = 1.0_dp/6.0_dp*u3
3416
3417 t1d = -8.0_dp + 4.0_dp*a1 - 0.5_dp*a2
3418 t2d = 10.0_dp - 8.0_dp*b1 + 1.5_dp*b2
3419 t3d = -2.0_dp + 4.0_dp*c1 - 1.5_dp*c2
3420 t4d = 0.5_dp*d2
3421 s1d = -8.0_dp + 4.0_dp*e1 - 0.5_dp*e2
3422 s2d = 10.0_dp - 8.0_dp*f1 + 1.5_dp*f2
3423 s3d = -2.0_dp + 4.0_dp*g1 - 1.5_dp*g2
3424 s4d = 0.5_dp*h2
3425 v1d = -8.0_dp + 4.0_dp*p1 - 0.5_dp*p2
3426 v2d = 10.0_dp - 8.0_dp*q1 + 1.5_dp*q2
3427 v3d = -2.0_dp + 4.0_dp*r1 - 1.5_dp*r2
3428 v4d = 0.5_dp*u2
3429
3430 t1 = t1d*dr1i
3431 t2 = t2d*dr1i
3432 t3 = t3d*dr1i
3433 t4 = t4d*dr1i
3434 s1 = s1o
3435 s2 = s2o
3436 s3 = s3o
3437 s4 = s4o
3438 v1 = v1o
3439 v2 = v2o
3440 v3 = v3o
3441 v4 = v4o
3442 val(1) = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3443 (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3444 (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3445 (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3446 ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3447 (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3448 (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3449 (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3450 ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3451 (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3452 (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3453 (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3454 ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3455 (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3456 (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3457 (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3458
3459 t1 = t1o
3460 t2 = t2o
3461 t3 = t3o
3462 t4 = t4o
3463 s1 = s1d*dr2i
3464 s2 = s2d*dr2i
3465 s3 = s3d*dr2i
3466 s4 = s4d*dr2i
3467 v1 = v1o
3468 v2 = v2o
3469 v3 = v3o
3470 v4 = v4o
3471 val(2) = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3472 (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3473 (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3474 (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3475 ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3476 (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3477 (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3478 (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3479 ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3480 (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3481 (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3482 (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3483 ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3484 (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3485 (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3486 (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3487
3488 t1 = t1o
3489 t2 = t2o
3490 t3 = t3o
3491 t4 = t4o
3492 s1 = s1o
3493 s2 = s2o
3494 s3 = s3o
3495 s4 = s4o
3496 v1 = v1d*dr3i
3497 v2 = v2d*dr3i
3498 v3 = v3d*dr3i
3499 v4 = v4d*dr3i
3500 val(3) = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3501 (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3502 (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3503 (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3504 ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3505 (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3506 (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3507 (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3508 ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3509 (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3510 (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3511 (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3512 ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3513 (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3514 (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3515 (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3516
3517 IF (my_mpsum) CALL pw%pw_grid%para%group%sum(val)
3518
3519 END FUNCTION eval_d_interp_spl3_pbc
3520
3521END MODULE pw_spline_utils
3522
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_copy
Definition pw_methods.F:145
pw_methods::pw_integral_ab
Definition pw_methods.F:221
pw_methods::pw_zero
Definition pw_methods.F:82
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_comm_congruent
integer, parameter, public mp_comm_congruent
Definition message_passing.F:827
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner.F:15
pw_grid_types
Definition pw_grid_types.F:13
pw_grid_types::pw_mode_local
integer, parameter, public pw_mode_local
Definition pw_grid_types.F:27
pw_grid_types::fullspace
integer, parameter, public fullspace
Definition pw_grid_types.F:26
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_pool_types::pw_pool_release
subroutine, public pw_pool_release(pool)
releases the given pool (see cp2k/doc/ReferenceCounting.html)
Definition pw_pool_types.F:259
pw_spline_utils
different utils that are useful to manipulate splines on the regular grid of a pw
Definition pw_spline_utils.F:18
pw_spline_utils::add_fine2coarse
subroutine, public add_fine2coarse(fine_values_pw, coarse_coeffs_pw, weights_1d, w_border0, w_border1, pbc, safe_computation)
low level function that adds a coarse grid (without boundary) to a fine grid.
Definition pw_spline_utils.F:1772
pw_spline_utils::precond_spl3_3
integer, parameter, public precond_spl3_3
Definition pw_spline_utils.F:44
pw_spline_utils::pw_spline_precond_release
subroutine, public pw_spline_precond_release(preconditioner)
releases the preconditioner
Definition pw_spline_utils.F:2581
pw_spline_utils::pw_spline_precond_create
subroutine, public pw_spline_precond_create(preconditioner, precond_kind, pool, pbc, transpose)
...
Definition pw_spline_utils.F:2475
pw_spline_utils::pw_nn_deriv_r
subroutine, public pw_nn_deriv_r(pw_in, pw_out, coeffs, idir)
calculates a nearest neighbor central derivative. for the x dir: pw_outarray(i,j,k)=( pw_in(i+1,...
Definition pw_spline_utils.F:899
pw_spline_utils::pw_spline_do_precond
subroutine, public pw_spline_do_precond(preconditioner, in_v, out_v)
applies the preconditioner to the system of equations to find the coefficients of the spline
Definition pw_spline_utils.F:2595
pw_spline_utils::pw_spline3_deriv_g
subroutine, public pw_spline3_deriv_g(spline_g, idir)
calculates the FFT of the values of the x,y,z (idir=1,2,3) derivative of the cubic spline
Definition pw_spline_utils.F:366
pw_spline_utils::pw_spline_precond_set_kind
subroutine, public pw_spline_precond_set_kind(preconditioner, precond_kind, pbc, transpose)
switches the types of precoditioner to use
Definition pw_spline_utils.F:2498
pw_spline_utils::spl3_1d_transf_coeffs
real(kind=dp), dimension(4), parameter, public spl3_1d_transf_coeffs
Definition pw_spline_utils.F:51
pw_spline_utils::spl3_1d_transf_border1
real(kind=dp), dimension(3), parameter, public spl3_1d_transf_border1
Definition pw_spline_utils.F:70
pw_spline_utils::spline2_coeffs
real(kind=dp), dimension(4), parameter, public spline2_coeffs
Definition pw_spline_utils.F:51
pw_spline_utils::add_coarse2fine
subroutine, public add_coarse2fine(coarse_coeffs_pw, fine_values_pw, weights_1d, w_border0, w_border1, pbc, safe_computation)
low level function that adds a coarse grid to a fine grid. If pbc is true periodic boundary condition...
Definition pw_spline_utils.F:977
pw_spline_utils::eval_d_interp_spl3_pbc
real(kind=dp) function, dimension(3), public eval_d_interp_spl3_pbc(vec, pw)
Evaluates the derivatives of the PBC interpolated Spline (pw) function on the generic input vector (v...
Definition pw_spline_utils.F:3300
pw_spline_utils::eval_interp_spl3_pbc
real(kind=dp) function, public eval_interp_spl3_pbc(vec, pw)
Evaluates the PBC interpolated Spline (pw) function on the generic input vector (vec)
Definition pw_spline_utils.F:3151
pw_spline_utils::spline3_deriv_coeffs
real(kind=dp), dimension(3), parameter, public spline3_deriv_coeffs
Definition pw_spline_utils.F:70
pw_spline_utils::precond_spl3_aint
integer, parameter, public precond_spl3_aint
Definition pw_spline_utils.F:44
pw_spline_utils::pw_spline_scale_deriv
subroutine, public pw_spline_scale_deriv(deriv_vals_r, transpose, scale)
rescales the derivatives from gridspacing=1 to the real derivatives
Definition pw_spline_utils.F:275
pw_spline_utils::spline2_deriv_coeffs
real(kind=dp), dimension(3), parameter, public spline2_deriv_coeffs
Definition pw_spline_utils.F:70
pw_spline_utils::pw_nn_smear_r
subroutine, public pw_nn_smear_r(pw_in, pw_out, coeffs)
calculates the values of a nearest neighbor smearing
Definition pw_spline_utils.F:856
pw_spline_utils::spline3_coeffs
real(kind=dp), dimension(4), parameter, public spline3_coeffs
Definition pw_spline_utils.F:51
pw_spline_utils::spl3_nopbc
subroutine, public spl3_nopbc(pw_in, pw_out)
...
Definition pw_spline_utils.F:3097
pw_spline_utils::find_coeffs
logical function, public find_coeffs(values, coeffs, linop, preconditioner, pool, eps_r, eps_x, max_iter, sumtype)
solves iteratively (CG) a systmes of linear equations linOp(coeffs)=values (for example those needed ...
Definition pw_spline_utils.F:2649
pw_spline_utils::nn50_deriv_coeffs
real(kind=dp), dimension(3), parameter, public nn50_deriv_coeffs
Definition pw_spline_utils.F:70
pw_spline_utils::pw_spline3_interpolate_values_g
subroutine, public pw_spline3_interpolate_values_g(spline_g)
calculates the FFT of the coefficients of the2 cubic spline that interpolates the given values
Definition pw_spline_utils.F:200
pw_spline_utils::pw_spline2_interpolate_values_g
subroutine, public pw_spline2_interpolate_values_g(spline_g)
calculates the FFT of the coefficients of the quadratic spline that interpolates the given values
Definition pw_spline_utils.F:129
pw_spline_utils::spl3_nopbct
subroutine, public spl3_nopbct(pw_in, pw_out)
...
Definition pw_spline_utils.F:3113
pw_spline_utils::nn10_deriv_coeffs
real(kind=dp), dimension(3), parameter, public nn10_deriv_coeffs
Definition pw_spline_utils.F:70
pw_spline_utils::spl3_precond1_coeff
real(kind=dp), dimension(4), parameter, public spl3_precond1_coeff
Definition pw_spline_utils.F:51
pw_spline_utils::spl3_1d_coeffs0
real(kind=dp), dimension(3), parameter, public spl3_1d_coeffs0
Definition pw_spline_utils.F:70
pw_spline_utils::pw_spline2_deriv_g
subroutine, public pw_spline2_deriv_g(spline_g, idir)
calculates the FFT of the values of the x,y,z (idir=1,2,3) derivative of the quadratic spline
Definition pw_spline_utils.F:501
pw_spline_utils::spl3_aint_coeff
real(kind=dp), dimension(4), parameter, public spl3_aint_coeff
Definition pw_spline_utils.F:51
pw_spline_utils::spl3_pbc
subroutine, public spl3_pbc(pw_in, pw_out)
...
Definition pw_spline_utils.F:3129
pw_spline_utils::no_precond
integer, parameter, public no_precond
Definition pw_spline_utils.F:44
pw_spline_utils::precond_spl3_2
integer, parameter, public precond_spl3_2
Definition pw_spline_utils.F:44
pw_spline_utils::nn10_coeffs
real(kind=dp), dimension(4), parameter, public nn10_coeffs
Definition pw_spline_utils.F:51
pw_spline_utils::precond_spl3_aint2
integer, parameter, public precond_spl3_aint2
Definition pw_spline_utils.F:44
pw_spline_utils::nn50_coeffs
real(kind=dp), dimension(4), parameter, public nn50_coeffs
Definition pw_spline_utils.F:51
pw_spline_utils::precond_spl3_1
integer, parameter, public precond_spl3_1
Definition pw_spline_utils.F:44
pw_types
Definition pw_types.F:24
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_spline_utils::pw_spline_precond_type
stores information for the preconditioner used to calculate the coeffs of splines
Definition pw_spline_utils.F:107
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62