Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2026 CP2K developers group <https://cp2k.org> */
4 : /* */
5 : /* SPDX-License-Identifier: BSD-3-Clause */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : #include <assert.h>
9 : #include <limits.h>
10 : #include <math.h>
11 : #include <stdio.h>
12 : #include <stdlib.h>
13 : #include <string.h>
14 :
15 : #if defined(__AVX2__) && defined(__FMA__)
16 : #include <immintrin.h>
17 : #endif
18 :
19 : #include "../common/grid_common.h"
20 : #include "../common/grid_library.h"
21 : #include "../common/grid_sphere_cache.h"
22 :
23 : #define GRID_MAX_LP_OPTIMIZED 9
24 :
25 : #if (GRID_DO_COLLOCATE)
26 : #define GRID_CONST_WHEN_COLLOCATE const
27 : #define GRID_CONST_WHEN_INTEGRATE
28 : #else
29 : #define GRID_CONST_WHEN_COLLOCATE
30 : #define GRID_CONST_WHEN_INTEGRATE const
31 : #endif
32 :
33 : /*******************************************************************************
34 : * \brief Simple loop body for ortho_cx_to_grid using plain C.
35 : * \author Ole Schuett
36 : ******************************************************************************/
37 : static inline void __attribute__((always_inline))
38 >32811*10^7 : ortho_cx_to_grid_scalar(const int lp, const int cmax, const int i,
39 >32811*10^7 : const double pol[3][lp + 1][2 * cmax + 1],
40 : GRID_CONST_WHEN_COLLOCATE double *cx,
41 : GRID_CONST_WHEN_INTEGRATE double *grid_0,
42 : GRID_CONST_WHEN_INTEGRATE double *grid_1,
43 : GRID_CONST_WHEN_INTEGRATE double *grid_2,
44 : GRID_CONST_WHEN_INTEGRATE double *grid_3) {
45 :
46 : #if (GRID_DO_COLLOCATE)
47 : // collocate
48 >17232*10^7 : double reg[4] = {0.0, 0.0, 0.0, 0.0};
49 >13785*10^8 : #pragma omp simd reduction(+ : reg)
50 : for (int lxp = 0; lxp <= lp; lxp++) {
51 >47347*10^7 : const double p = pol[0][lxp][i + cmax];
52 >47347*10^7 : reg[0] += cx[lxp * 4 + 0] * p;
53 >47347*10^7 : reg[1] += cx[lxp * 4 + 1] * p;
54 >47347*10^7 : reg[2] += cx[lxp * 4 + 2] * p;
55 >47347*10^7 : reg[3] += cx[lxp * 4 + 3] * p;
56 : }
57 >17232*10^7 : *grid_0 += reg[0];
58 >17232*10^7 : *grid_1 += reg[1];
59 >17232*10^7 : *grid_2 += reg[2];
60 >17232*10^7 : *grid_3 += reg[3];
61 :
62 : #else
63 : // integrate
64 >15579*10^7 : const double reg[4] = {*grid_0, *grid_1, *grid_2, *grid_3};
65 : #pragma omp simd
66 >15579*10^7 : for (int lxp = 0; lxp <= lp; lxp++) {
67 >43568*10^7 : const double p = pol[0][lxp][i + cmax];
68 >43568*10^7 : cx[lxp * 4 + 0] += reg[0] * p;
69 >43568*10^7 : cx[lxp * 4 + 1] += reg[1] * p;
70 >43568*10^7 : cx[lxp * 4 + 2] += reg[2] * p;
71 >43568*10^7 : cx[lxp * 4 + 3] += reg[3] * p;
72 : }
73 : #endif
74 : }
75 :
76 : /*******************************************************************************
77 : * \brief Optimized loop body for ortho_cx_to_grid using AVX2 Intel Intrinsics.
78 : * This routine always processes four consecutive grid elements at once.
79 : * \author Ole Schuett
80 : ******************************************************************************/
81 : #if defined(__AVX2__) && defined(__FMA__)
82 : static inline void __attribute__((always_inline))
83 : ortho_cx_to_grid_avx2(const int lp, const int cmax, const int i,
84 : const double pol[3][lp + 1][2 * cmax + 1],
85 : GRID_CONST_WHEN_COLLOCATE double *cx,
86 : GRID_CONST_WHEN_INTEGRATE double *grid_0,
87 : GRID_CONST_WHEN_INTEGRATE double *grid_1,
88 : GRID_CONST_WHEN_INTEGRATE double *grid_2,
89 : GRID_CONST_WHEN_INTEGRATE double *grid_3) {
90 :
91 : const int icmax = i + cmax;
92 :
93 : #if (GRID_DO_COLLOCATE)
94 : // collocate
95 : // First iteration for lxp == 0 does not need add instructions.
96 : __m256d p_vec = _mm256_loadu_pd(&pol[0][0][icmax]);
97 : __m256d r_vec_0 = _mm256_mul_pd(p_vec, _mm256_set1_pd(cx[0]));
98 : __m256d r_vec_1 = _mm256_mul_pd(p_vec, _mm256_set1_pd(cx[1]));
99 : __m256d r_vec_2 = _mm256_mul_pd(p_vec, _mm256_set1_pd(cx[2]));
100 : __m256d r_vec_3 = _mm256_mul_pd(p_vec, _mm256_set1_pd(cx[3]));
101 :
102 : // Remaining iterations for lxp > 0 use fused multiply adds.
103 : GRID_PRAGMA_UNROLL_UP_TO(GRID_MAX_LP_OPTIMIZED)
104 : for (int lxp = 1; lxp <= lp; lxp++) {
105 : const double *cx_base = &cx[lxp * 4];
106 : p_vec = _mm256_loadu_pd(&pol[0][lxp][icmax]);
107 : r_vec_0 = _mm256_fmadd_pd(p_vec, _mm256_set1_pd(cx_base[0]), r_vec_0);
108 : r_vec_1 = _mm256_fmadd_pd(p_vec, _mm256_set1_pd(cx_base[1]), r_vec_1);
109 : r_vec_2 = _mm256_fmadd_pd(p_vec, _mm256_set1_pd(cx_base[2]), r_vec_2);
110 : r_vec_3 = _mm256_fmadd_pd(p_vec, _mm256_set1_pd(cx_base[3]), r_vec_3);
111 : }
112 :
113 : // Add vectors to grid one at a time, because they can aliase when cube wraps.
114 : _mm256_storeu_pd(grid_0, _mm256_add_pd(_mm256_loadu_pd(grid_0), r_vec_0));
115 : _mm256_storeu_pd(grid_1, _mm256_add_pd(_mm256_loadu_pd(grid_1), r_vec_1));
116 : _mm256_storeu_pd(grid_2, _mm256_add_pd(_mm256_loadu_pd(grid_2), r_vec_2));
117 : _mm256_storeu_pd(grid_3, _mm256_add_pd(_mm256_loadu_pd(grid_3), r_vec_3));
118 :
119 : #else
120 : // integrate
121 : __m256d grid_vec_0 = _mm256_loadu_pd(grid_0);
122 : __m256d grid_vec_1 = _mm256_loadu_pd(grid_1);
123 : __m256d grid_vec_2 = _mm256_loadu_pd(grid_2);
124 : __m256d grid_vec_3 = _mm256_loadu_pd(grid_3);
125 :
126 : GRID_PRAGMA_UNROLL_UP_TO(GRID_MAX_LP_OPTIMIZED + 1)
127 : for (int lxp = 0; lxp <= lp; lxp++) {
128 : __m256d p_vec = _mm256_loadu_pd(&pol[0][lxp][icmax]);
129 :
130 : // Do 4 dot products at once. https://stackoverflow.com/a/10454420
131 : __m256d xy0 = _mm256_mul_pd(p_vec, grid_vec_0);
132 : __m256d xy1 = _mm256_mul_pd(p_vec, grid_vec_1);
133 : __m256d xy2 = _mm256_mul_pd(p_vec, grid_vec_2);
134 : __m256d xy3 = _mm256_mul_pd(p_vec, grid_vec_3);
135 :
136 : // low to high: xy00+xy01 xy10+xy11 xy02+xy03 xy12+xy13
137 : __m256d temp01 = _mm256_hadd_pd(xy0, xy1);
138 :
139 : // low to high: xy20+xy21 xy30+xy31 xy22+xy23 xy32+xy33
140 : __m256d temp23 = _mm256_hadd_pd(xy2, xy3);
141 :
142 : // low to high: xy02+xy03 xy12+xy13 xy20+xy21 xy30+xy31
143 : __m256d swapped = _mm256_permute2f128_pd(temp01, temp23, 0x21);
144 :
145 : // low to high: xy00+xy01 xy10+xy11 xy22+xy23 xy32+xy33
146 : __m256d blended = _mm256_blend_pd(temp01, temp23, 0b1100);
147 :
148 : __m256d r_vec = _mm256_add_pd(swapped, blended);
149 :
150 : // cx += r_vec
151 : double *cx_base = &cx[lxp * 4];
152 : _mm256_storeu_pd(cx_base, _mm256_add_pd(r_vec, _mm256_loadu_pd(cx_base)));
153 : }
154 : #endif
155 : }
156 : #endif // __AVX2__ && __FMA__
157 :
158 : /*******************************************************************************
159 : * \brief Collocates coefficients C_x onto the grid for orthorhombic case.
160 : * \author Ole Schuett
161 : ******************************************************************************/
162 : static inline void __attribute__((always_inline))
163 17608027260 : ortho_cx_to_grid(const int lp, const int kg1, const int kg2, const int jg1,
164 : const int jg2, const int cmax,
165 17608027260 : const double pol[3][lp + 1][2 * cmax + 1],
166 17608027260 : const int map[3][2 * cmax + 1],
167 17608027260 : const int sections[3][2 * cmax + 1], const int npts_local[3],
168 : int **sphere_bounds_iter, GRID_CONST_WHEN_COLLOCATE double *cx,
169 : GRID_CONST_WHEN_INTEGRATE double *grid) {
170 :
171 : // Lower and upper sphere bounds relative to center, ie. in cube coordinates.
172 17608027260 : const int lb = *((*sphere_bounds_iter)++);
173 17608027260 : const int ub = 1 - lb;
174 :
175 : // AVX instructions can only load/store from evenly spaced memory locations.
176 : // Since the sphere bounds can wrap around due to the grid's periodicity,
177 : // the inner loop runs over sections with homogeneous cube to grid mapping.
178 50740777890 : for (int istart = lb; istart <= ub; istart++) {
179 33132750630 : const int istop = imin(ub, istart + sections[0][istart + cmax]);
180 33132750630 : const int cube2grid = map[0][istart + cmax] - istart;
181 :
182 33132750630 : const int stride = npts_local[1] * npts_local[0];
183 33132750630 : const int grid_index_0 = kg1 * stride + jg1 * npts_local[0];
184 33132750630 : const int grid_index_1 = kg2 * stride + jg1 * npts_local[0];
185 33132750630 : const int grid_index_2 = kg1 * stride + jg2 * npts_local[0];
186 33132750630 : const int grid_index_3 = kg2 * stride + jg2 * npts_local[0];
187 33132750630 : GRID_CONST_WHEN_INTEGRATE double *grid_base_0 = &grid[grid_index_0];
188 33132750630 : GRID_CONST_WHEN_INTEGRATE double *grid_base_1 = &grid[grid_index_1];
189 33132750630 : GRID_CONST_WHEN_INTEGRATE double *grid_base_2 = &grid[grid_index_2];
190 33132750630 : GRID_CONST_WHEN_INTEGRATE double *grid_base_3 = &grid[grid_index_3];
191 :
192 : // Use AVX2 to process grid points in chunks of four, ie. 256 bit vectors.
193 : #if defined(__AVX2__) && defined(__FMA__)
194 : const int istop_vec = istart + 4 * ((istop - istart + 1) / 4) - 1;
195 : for (int i = istart; i <= istop_vec; i += 4) {
196 : const int ig = i + cube2grid;
197 : ortho_cx_to_grid_avx2(lp, cmax, i, pol, cx, &grid_base_0[ig],
198 : &grid_base_1[ig], &grid_base_2[ig],
199 : &grid_base_3[ig]);
200 : }
201 : istart = istop_vec + 1;
202 : #endif
203 :
204 : // Process up to 3 remaining points - or everything if AVX2 isn't available.
205 >36125*10^7 : for (int i = istart; i <= istop; i++) {
206 >32811*10^7 : const int ig = i + cube2grid;
207 >32811*10^7 : ortho_cx_to_grid_scalar(lp, cmax, i, pol, cx, &grid_base_0[ig],
208 : &grid_base_1[ig], &grid_base_2[ig],
209 >32811*10^7 : &grid_base_3[ig]);
210 : }
211 33132750630 : istart = istop;
212 : }
213 : }
214 :
215 : /*******************************************************************************
216 : * \brief Transforms coefficients C_xy into C_x by fixing grid index j.
217 : * \author Ole Schuett
218 : ******************************************************************************/
219 : static inline void __attribute__((always_inline))
220 17608027260 : ortho_cxy_to_cx(const int lp, const int j1, const int j2, const int cmax,
221 17608027260 : const double pol[3][lp + 1][2 * cmax + 1],
222 : GRID_CONST_WHEN_COLLOCATE double *cxy,
223 : GRID_CONST_WHEN_INTEGRATE double *cx) {
224 :
225 66383873777 : for (int lyp = 0; lyp <= lp; lyp++) {
226 >15299*10^7 : for (int lxp = 0; lxp <= lp - lyp; lxp++) {
227 >10422*10^7 : const double p1 = pol[1][lyp][j1 + cmax];
228 >10422*10^7 : const double p2 = pol[1][lyp][j2 + cmax];
229 >10422*10^7 : const int cxy_index = lyp * (lp + 1) * 2 + lxp * 2; // [lyp, lxp, 0]
230 :
231 : #if (GRID_DO_COLLOCATE)
232 : // collocate
233 53598473151 : cx[lxp * 4 + 0] += cxy[cxy_index + 0] * p1;
234 53598473151 : cx[lxp * 4 + 1] += cxy[cxy_index + 1] * p1;
235 53598473151 : cx[lxp * 4 + 2] += cxy[cxy_index + 0] * p2;
236 53598473151 : cx[lxp * 4 + 3] += cxy[cxy_index + 1] * p2;
237 : #else
238 : // integrate
239 50624561796 : cxy[cxy_index + 0] += cx[lxp * 4 + 0] * p1;
240 50624561796 : cxy[cxy_index + 1] += cx[lxp * 4 + 1] * p1;
241 50624561796 : cxy[cxy_index + 0] += cx[lxp * 4 + 2] * p2;
242 50624561796 : cxy[cxy_index + 1] += cx[lxp * 4 + 3] * p2;
243 : #endif
244 : }
245 : }
246 : }
247 :
248 : /*******************************************************************************
249 : * \brief Loop body of ortho_cxy_to_grid to be inlined for low values of lp.
250 : * \author Ole Schuett
251 : ******************************************************************************/
252 17608027260 : static inline void __attribute__((always_inline)) ortho_cxy_to_grid_low(
253 : const int lp, const int j1, const int j2, const int kg1, const int kg2,
254 : const int jg1, const int jg2, const int cmax,
255 17608027260 : const double pol[3][lp + 1][2 * cmax + 1], const int map[3][2 * cmax + 1],
256 17608027260 : const int sections[3][2 * cmax + 1], const int npts_local[3],
257 : int **sphere_bounds_iter, double *cx, GRID_CONST_WHEN_COLLOCATE double *cxy,
258 : GRID_CONST_WHEN_INTEGRATE double *grid) {
259 :
260 : #if (GRID_DO_COLLOCATE)
261 : // collocate
262 18372391086 : ortho_cxy_to_cx(lp, j1, j2, cmax, pol, cxy, cx);
263 9186195543 : ortho_cx_to_grid(lp, kg1, kg2, jg1, jg2, cmax, pol, map, sections, npts_local,
264 : sphere_bounds_iter, cx, grid);
265 : #else
266 : // integrate
267 16843663434 : ortho_cx_to_grid(lp, kg1, kg2, jg1, jg2, cmax, pol, map, sections, npts_local,
268 : sphere_bounds_iter, cx, grid);
269 8421831717 : ortho_cxy_to_cx(lp, j1, j2, cmax, pol, cxy, cx);
270 : #endif
271 : }
272 :
273 : /*******************************************************************************
274 : * \brief Collocates coefficients C_xy onto the grid for orthorhombic case.
275 : * \author Ole Schuett
276 : ******************************************************************************/
277 1933958092 : static inline void ortho_cxy_to_grid(
278 : const int lp, const int kg1, const int kg2, const int cmax,
279 1933958092 : const double pol[3][lp + 1][2 * cmax + 1], const int map[3][2 * cmax + 1],
280 1933958092 : const int sections[3][2 * cmax + 1], const int npts_local[3],
281 : int **sphere_bounds_iter, GRID_CONST_WHEN_COLLOCATE double *cxy,
282 1933958092 : GRID_CONST_WHEN_INTEGRATE double *grid) {
283 :
284 : // The cube contains an even number of grid points in each direction and
285 : // collocation is always performed on a pair of two opposing grid points.
286 : // Hence, the points with index 0 and 1 are both assigned distance zero via
287 : // the formular distance=(2*index-1)/2.
288 :
289 1933958092 : const int jstart = *((*sphere_bounds_iter)++);
290 1933958092 : const size_t cx_size = (lp + 1) * 4;
291 1933958092 : double cx[cx_size];
292 19541985352 : for (int j1 = jstart; j1 <= 0; j1++) {
293 17608027260 : const int j2 = 1 - j1;
294 17608027260 : const int jg1 = map[1][j1 + cmax];
295 17608027260 : const int jg2 = map[1][j2 + cmax];
296 :
297 17608027260 : memset(cx, 0, cx_size * sizeof(double));
298 :
299 : // Generate separate branches for low values of lp gives up to 30% speedup.
300 17608027260 : if (lp <= GRID_MAX_LP_OPTIMIZED) {
301 : GRID_PRAGMA_UNROLL(GRID_MAX_LP_OPTIMIZED + 1)
302 >19368*10^7 : for (int ilp = 0; ilp <= GRID_MAX_LP_OPTIMIZED; ilp++) {
303 >17608*10^7 : if (lp == ilp) {
304 >19368*10^7 : ortho_cxy_to_grid_low(ilp, j1, j2, kg1, kg2, jg1, jg2, cmax, pol, map,
305 : sections, npts_local, sphere_bounds_iter, cx,
306 : cxy, grid);
307 : }
308 : }
309 : } else {
310 17608030878 : ortho_cxy_to_grid_low(lp, j1, j2, kg1, kg2, jg1, jg2, cmax, pol, map,
311 : sections, npts_local, sphere_bounds_iter, cx, cxy,
312 : grid);
313 : }
314 : }
315 1933958092 : }
316 :
317 : /*******************************************************************************
318 : * \brief Transforms coefficients C_xyz into C_xz by fixing grid index k.
319 : * \author Ole Schuett
320 : ******************************************************************************/
321 1933958092 : static inline void ortho_cxyz_to_cxy(const int lp, const int k1, const int k2,
322 : const int cmax,
323 1933958092 : const double pol[3][lp + 1][2 * cmax + 1],
324 : GRID_CONST_WHEN_COLLOCATE double *cxyz,
325 : GRID_CONST_WHEN_INTEGRATE double *cxy) {
326 :
327 7224325341 : for (int lzp = 0; lzp <= lp; lzp++) {
328 16501116777 : for (int lyp = 0; lyp <= lp - lzp; lyp++) {
329 31922543123 : for (int lxp = 0; lxp <= lp - lzp - lyp; lxp++) {
330 20711793595 : const double p1 = pol[2][lzp][k1 + cmax];
331 20711793595 : const double p2 = pol[2][lzp][k2 + cmax];
332 20711793595 : const int cxyz_index =
333 20711793595 : lzp * (lp + 1) * (lp + 1) + lyp * (lp + 1) + lxp; // [lzp, lyp, lxp]
334 20711793595 : const int cxy_index = lyp * (lp + 1) * 2 + lxp * 2; // [lyp, lxp, 0]
335 :
336 : #if (GRID_DO_COLLOCATE)
337 : // collocate
338 10445325220 : cxy[cxy_index + 0] += cxyz[cxyz_index] * p1;
339 10445325220 : cxy[cxy_index + 1] += cxyz[cxyz_index] * p2;
340 : #else
341 : // integrate
342 10266468375 : cxyz[cxyz_index] += cxy[cxy_index + 0] * p1;
343 10266468375 : cxyz[cxyz_index] += cxy[cxy_index + 1] * p2;
344 : #endif
345 : }
346 : }
347 : }
348 1933958092 : }
349 :
350 : /*******************************************************************************
351 : * \brief Collocates coefficients C_xyz onto the grid for orthorhombic case.
352 : * \author Ole Schuett
353 : ******************************************************************************/
354 : static inline void
355 181657159 : ortho_cxyz_to_grid(const int lp, const double zetp, const double dh[3][3],
356 : const double dh_inv[3][3], const double rp[3],
357 : const int npts_global[3], const int npts_local[3],
358 : const int shift_local[3], const double radius,
359 : GRID_CONST_WHEN_COLLOCATE double *cxyz,
360 181657159 : GRID_CONST_WHEN_INTEGRATE double *grid) {
361 :
362 : // *** position of the gaussian product
363 : //
364 : // this is the actual definition of the position on the grid
365 : // i.e. a point rp(:) gets here grid coordinates
366 : // MODULO(rp(:)/dr(:),npts_global(:))+1
367 : // hence (0.0,0.0,0.0) in real space is rsgrid%lb on the rsgrid in Fortran
368 : // and (1,1,1) on grid here in C.
369 :
370 : // cubecenter(:) = FLOOR(MATMUL(dh_inv, rp))
371 181657159 : int cubecenter[3];
372 726628636 : for (int i = 0; i < 3; i++) {
373 : double dh_inv_rp = 0.0;
374 2179885908 : for (int j = 0; j < 3; j++) {
375 1634914431 : dh_inv_rp += dh_inv[j][i] * rp[j];
376 : }
377 544971477 : cubecenter[i] = (int)floor(dh_inv_rp);
378 : }
379 :
380 : double roffset[3];
381 726628636 : for (int i = 0; i < 3; i++) {
382 544971477 : roffset[i] = rp[i] - ((double)cubecenter[i]) * dh[i][i];
383 : }
384 :
385 : // Lookup loop bounds for spherical cutoff.
386 181657159 : int *sphere_bounds;
387 181657159 : double disr_radius;
388 181657159 : grid_sphere_cache_lookup(radius, dh, dh_inv, &sphere_bounds, &disr_radius);
389 181657159 : int **sphere_bounds_iter = &sphere_bounds;
390 :
391 : // Cube bounds.
392 181657159 : int lb_cube[3], ub_cube[3];
393 726628636 : for (int i = 0; i < 3; i++) {
394 544971477 : lb_cube[i] = (int)ceil(-1e-8 - disr_radius * dh_inv[i][i]);
395 544971477 : ub_cube[i] = 1 - lb_cube[i];
396 : // If grid is not period check that cube fits without wrapping.
397 544971477 : if (npts_global[i] != npts_local[i]) {
398 20802 : const int offset =
399 20802 : modulo(cubecenter[i] + lb_cube[i] - shift_local[i], npts_global[i]) -
400 : lb_cube[i];
401 20802 : assert(offset + ub_cube[i] < npts_local[i]);
402 20802 : assert(offset + lb_cube[i] >= 0);
403 : }
404 : }
405 :
406 : // cmax = MAXVAL(ub_cube)
407 181657159 : const int cmax = imax(imax(ub_cube[0], ub_cube[1]), ub_cube[2]);
408 :
409 : // Precompute (x-xp)**lp*exp(..) for each direction.
410 181657159 : double pol_mutable[3][lp + 1][2 * cmax + 1];
411 726628636 : for (int idir = 0; idir < 3; idir++) {
412 544971477 : const double dr = dh[idir][idir];
413 544971477 : const double ro = roffset[idir];
414 : // Reuse the result from the previous gridpoint to avoid to many exps:
415 : // exp( -a*(x+d)**2) = exp(-a*x**2)*exp(-2*a*x*d)*exp(-a*d**2)
416 : // exp(-2*a*(x+d)*d) = exp(-2*a*x*d)*exp(-2*a*d**2)
417 544971477 : const double t_exp_1 = exp(-zetp * pow(dr, 2));
418 544971477 : const double t_exp_2 = pow(t_exp_1, 2);
419 544971477 : double t_exp_min_1 = exp(-zetp * pow(+dr - ro, 2));
420 544971477 : double t_exp_min_2 = exp(-2 * zetp * (+dr - ro) * (-dr));
421 6356587719 : for (int ig = 0; ig >= lb_cube[idir]; ig--) {
422 5811616242 : const double rpg = ig * dr - ro;
423 5811616242 : t_exp_min_1 *= t_exp_min_2 * t_exp_1;
424 5811616242 : t_exp_min_2 *= t_exp_2;
425 5811616242 : double pg = t_exp_min_1;
426 21710827601 : for (int icoef = 0; icoef <= lp; icoef++) {
427 15899211359 : pol_mutable[idir][icoef][ig + cmax] = pg; // exp(-zetp*rpg**2)
428 15899211359 : pg *= rpg;
429 : }
430 : }
431 544971477 : double t_exp_plus_1 = exp(-zetp * pow(-ro, 2));
432 544971477 : double t_exp_plus_2 = exp(-2 * zetp * (-ro) * (+dr));
433 6356587719 : for (int ig = 0; ig >= lb_cube[idir]; ig--) {
434 5811616242 : const double rpg = (1 - ig) * dr - ro;
435 5811616242 : t_exp_plus_1 *= t_exp_plus_2 * t_exp_1;
436 5811616242 : t_exp_plus_2 *= t_exp_2;
437 5811616242 : double pg = t_exp_plus_1;
438 21710827601 : for (int icoef = 0; icoef <= lp; icoef++) {
439 15899211359 : pol_mutable[idir][icoef][1 - ig + cmax] = pg; // exp(-zetp*rpg**2)
440 15899211359 : pg *= rpg;
441 : }
442 : }
443 : }
444 181657159 : const double (*pol)[lp + 1][2 * cmax + 1] =
445 181657159 : (const double (*)[lp + 1][2 * cmax + 1]) pol_mutable;
446 :
447 : // Precompute mapping from cube to grid indices for each direction
448 181657159 : int map_mutable[3][2 * cmax + 1];
449 726628636 : for (int i = 0; i < 3; i++) {
450 12797724180 : for (int k = -cmax; k <= +cmax; k++) {
451 12252752703 : map_mutable[i][k + cmax] =
452 12252752703 : modulo(cubecenter[i] + k - shift_local[i], npts_global[i]);
453 : }
454 : }
455 181657159 : const int (*map)[2 * cmax + 1] = (const int (*)[2 * cmax + 1]) map_mutable;
456 :
457 : // Precompute length of sections with homogeneous cube to grid mapping.
458 181657159 : int sections_mutable[3][2 * cmax + 1];
459 726628636 : for (int i = 0; i < 3; i++) {
460 12797724180 : for (int kg = 2 * cmax; kg >= 0; kg--) {
461 12252752703 : if (kg == 2 * cmax || map[i][kg] != map[i][kg + 1] - 1) {
462 1173000143 : sections_mutable[i][kg] = 0;
463 : } else {
464 11079752560 : sections_mutable[i][kg] = sections_mutable[i][kg + 1] + 1;
465 : }
466 : }
467 : }
468 181657159 : const int (*sections)[2 * cmax + 1] =
469 181657159 : (const int (*)[2 * cmax + 1]) sections_mutable;
470 :
471 : // Loop over k dimension of the cube.
472 181657159 : const int kstart = *((*sphere_bounds_iter)++);
473 181657159 : const size_t cxy_size = (lp + 1) * (lp + 1) * 2;
474 181657159 : double cxy[cxy_size];
475 2115615251 : for (int k1 = kstart; k1 <= 0; k1++) {
476 1933958092 : const int k2 = 1 - k1;
477 1933958092 : const int kg1 = map[2][k1 + cmax];
478 1933958092 : const int kg2 = map[2][k2 + cmax];
479 :
480 1933958092 : memset(cxy, 0, cxy_size * sizeof(double));
481 :
482 : #if (GRID_DO_COLLOCATE)
483 : // collocate
484 996603654 : ortho_cxyz_to_cxy(lp, k1, k2, cmax, pol, cxyz, cxy);
485 996603654 : ortho_cxy_to_grid(lp, kg1, kg2, cmax, pol, map, sections, npts_local,
486 : sphere_bounds_iter, cxy, grid);
487 : #else
488 : // integrate
489 937354438 : ortho_cxy_to_grid(lp, kg1, kg2, cmax, pol, map, sections, npts_local,
490 : sphere_bounds_iter, cxy, grid);
491 937354438 : ortho_cxyz_to_cxy(lp, k1, k2, cmax, pol, cxyz, cxy);
492 : #endif
493 : }
494 181657159 : }
495 :
496 : /*******************************************************************************
497 : * \brief Collocates coefficients C_i onto the grid for general case.
498 : * \author Ole Schuett
499 : ******************************************************************************/
500 5859092566 : static inline void __attribute__((always_inline)) general_ci_to_grid(
501 : const int lp, const int jg, const int kg, const int ismin, const int ismax,
502 : const int npts_local[3], const int index_min[3], const int index_max[3],
503 : const int map_i[], const int sections_i[], const double gp[3], const int k,
504 : const int j, const double exp_ij[], const double exp_jk[],
505 : const double exp_ik[], GRID_CONST_WHEN_COLLOCATE double *ci,
506 : GRID_CONST_WHEN_INTEGRATE double *grid) {
507 :
508 5859092566 : const int base = kg * npts_local[1] * npts_local[0] + jg * npts_local[0];
509 5859092566 : const int stride_i = index_max[0] - index_min[0] + 1;
510 5859092566 : const int stride_j = index_max[1] - index_min[1] + 1;
511 5859092566 : const int idx_jk = (k - index_min[2]) * stride_j + j - index_min[1];
512 :
513 : #if defined(__AVX2__) && defined(__FMA__)
514 : const __m256d v_gp0 = _mm256_set1_pd(gp[0]);
515 : const __m256d v_four = _mm256_set1_pd(4.0);
516 : const __m256d v_ejk = _mm256_set1_pd(exp_jk[idx_jk]);
517 : __m256d v_ci[lp + 1];
518 : GRID_PRAGMA_UNROLL_UP_TO(GRID_MAX_LP_OPTIMIZED)
519 : for (int il = 0; il <= lp; il++) {
520 : #if (GRID_DO_COLLOCATE)
521 : v_ci[il] = _mm256_set1_pd(ci[il]);
522 : #else
523 : v_ci[il] = _mm256_setzero_pd();
524 : #endif
525 : }
526 : #endif
527 :
528 : // AVX instructions can only load/store from evenly spaced memory locations.
529 : // Since the cube can wrap around due to the grid's periodicity,
530 : // the inner loop runs over sections with homogeneous cube to grid mapping.
531 18568204154 : for (int istart = ismin; istart <= ismax; istart++) {
532 12709111588 : const int istop = imin(ismax, istart + sections_i[istart - index_min[0]]);
533 12709111588 : if (map_i[istart - index_min[0]] < 0) {
534 0 : istart = istop; // skip over out-of-bounds indicies
535 0 : continue;
536 : }
537 :
538 12709111588 : const int cube2grid = map_i[istart - index_min[0]] - istart;
539 :
540 : #if defined(__AVX2__) && defined(__FMA__)
541 : const int istart_tail = istart + ((istop - istart + 1) / 4) * 4;
542 : __m256d v_i = _mm256_set_pd(istart + 3, istart + 2, istart + 1, istart);
543 : for (int i = istart; i < istart_tail; i += 4) {
544 : const int idx_ij = (j - index_min[1]) * stride_i + i - index_min[0];
545 : const int idx_ik = (k - index_min[2]) * stride_i + i - index_min[0];
546 :
547 : const __m256d v_di = _mm256_sub_pd(v_i, v_gp0);
548 : v_i = _mm256_add_pd(v_i, v_four);
549 : const __m256d v_eij = _mm256_loadu_pd(&exp_ij[idx_ij]);
550 : const __m256d v_eik = _mm256_loadu_pd(&exp_ik[idx_ik]);
551 : __m256d v_dip = _mm256_mul_pd(_mm256_mul_pd(v_eij, v_ejk), v_eik);
552 :
553 : const int grid_index = base + i + cube2grid;
554 : #if (GRID_DO_COLLOCATE)
555 : __m256d v_reg = _mm256_mul_pd(v_ci[0], v_dip);
556 : GRID_PRAGMA_UNROLL_UP_TO(GRID_MAX_LP_OPTIMIZED)
557 : for (int il = 1; il <= lp; il++) {
558 : v_dip = _mm256_mul_pd(v_dip, v_di);
559 : v_reg = _mm256_fmadd_pd(v_ci[il], v_dip, v_reg);
560 : }
561 : __m256d v_grid = _mm256_loadu_pd(&grid[grid_index]);
562 : _mm256_storeu_pd(&grid[grid_index], _mm256_add_pd(v_grid, v_reg));
563 : #else
564 : const __m256d v_reg = _mm256_loadu_pd(&grid[grid_index]);
565 : v_ci[0] = _mm256_fmadd_pd(v_reg, v_dip, v_ci[0]);
566 : GRID_PRAGMA_UNROLL_UP_TO(GRID_MAX_LP_OPTIMIZED)
567 : for (int il = 1; il <= lp; il++) {
568 : v_dip = _mm256_mul_pd(v_dip, v_di);
569 : v_ci[il] = _mm256_fmadd_pd(v_reg, v_dip, v_ci[il]);
570 : }
571 : #endif
572 : }
573 : #else
574 12709111588 : const int istart_tail = istart;
575 : #endif
576 :
577 >11679*10^7 : for (int i = istart_tail; i <= istop; i++) {
578 >10408*10^7 : const int ig = i + cube2grid;
579 >10408*10^7 : const double di = i - gp[0];
580 :
581 >10408*10^7 : const int idx_ij = (j - index_min[1]) * stride_i + i - index_min[0];
582 >10408*10^7 : const int idx_ik = (k - index_min[2]) * stride_i + i - index_min[0];
583 :
584 : // Mathieu's trick: Calculate 3D Gaussian from three precomputed 2D tables
585 : //
586 : // r = (i-gp[0])*dh[0,:] + (j-gp[1])*dh[1,:] + (k-gp[2])*dh[2,:]
587 : // = a + b + c
588 : //
589 : // r**2 = (a + b + c)**2 = a**2 + b**2 + c**2 + 2ab + 2bc + 2ca
590 : //
591 : // exp(-r**2) = exp(-0.5*a**2 - 2ab - 0.5*b**2)
592 : // * exp(-0.5*b**2 - 2bc - 0.5*c**2)
593 : // * exp(-0.5*a**2 - 2ac - 0.5*c**2)
594 : //
595 >10408*10^7 : const double gaussian = exp_ij[idx_ij] * exp_jk[idx_jk] * exp_ik[idx_ik];
596 :
597 >10408*10^7 : const int grid_index = base + ig; // [kg, jg, ig]
598 >10408*10^7 : double dip = gaussian;
599 :
600 : #if (GRID_DO_COLLOCATE)
601 : // collocate
602 57697642572 : double reg = 0.0;
603 >24032*10^7 : for (int il = 0; il <= lp; il++) {
604 >18262*10^7 : reg += ci[il] * dip;
605 >18262*10^7 : dip *= di;
606 : }
607 57697642572 : grid[grid_index] += reg;
608 : #else
609 : // integrate
610 46389547583 : const double reg = grid[grid_index];
611 >19259*10^7 : for (int il = 0; il <= lp; il++) {
612 >14621*10^7 : ci[il] += reg * dip;
613 >14621*10^7 : dip *= di;
614 : }
615 : #endif
616 : }
617 : istart = istop;
618 : }
619 :
620 : #if defined(__AVX2__) && defined(__FMA__) && (!GRID_DO_COLLOCATE)
621 : GRID_PRAGMA_UNROLL_UP_TO(GRID_MAX_LP_OPTIMIZED)
622 : for (int il = 0; il <= lp; il++) {
623 : __m128d v_low = _mm256_castpd256_pd128(v_ci[il]);
624 : __m128d v_high = _mm256_extractf128_pd(v_ci[il], 1);
625 : __m128d v_sum2 = _mm_add_pd(v_low, v_high);
626 : __m128d v_shuf = _mm_unpackhi_pd(v_sum2, v_sum2);
627 : ci[il] += _mm_cvtsd_f64(_mm_add_sd(v_sum2, v_shuf));
628 : }
629 : #endif
630 : }
631 :
632 : /*******************************************************************************
633 : * \brief Transforms coefficients C_ij into C_i by fixing grid index j.
634 : * \author Ole Schuett
635 : ******************************************************************************/
636 : static inline void __attribute__((always_inline))
637 3259667752 : general_cij_to_ci(const int lp, const double dj,
638 : GRID_CONST_WHEN_COLLOCATE double *cij,
639 : GRID_CONST_WHEN_INTEGRATE double *ci) {
640 3259667752 : double djp = 1.0;
641 32460877539 : for (int jl = 0; jl <= lp; jl++) {
642 61580354583 : for (int il = 0; il <= lp - jl; il++) {
643 43143699969 : const int cij_index = jl * (lp + 1) + il; // [jl, il]
644 : #if (GRID_DO_COLLOCATE)
645 23818594426 : ci[il] += cij[cij_index] * djp; // collocate
646 : #else
647 19325105543 : cij[cij_index] += ci[il] * djp; // integrate
648 : #endif
649 : }
650 18436654614 : djp *= dj;
651 : }
652 : }
653 :
654 : /*******************************************************************************
655 : * \brief Loop body of general_cij_to_grid to be inlined for low values of lp.
656 : * \author Ole Schuett
657 : ******************************************************************************/
658 5859092566 : static inline void __attribute__((always_inline)) general_cij_to_grid_low(
659 : const int lp, const int jg, const int kg, const int ismin, const int ismax,
660 : const int npts_local[3], const int index_min[3], const int index_max[3],
661 : const int map_i[], const int sections_i[], const double gp[3], const int k,
662 : const int j, const double exp_ij[], const double exp_jk[],
663 : const double exp_ik[], const double dj, double *ci,
664 : GRID_CONST_WHEN_COLLOCATE double *cij,
665 : GRID_CONST_WHEN_INTEGRATE double *grid) {
666 :
667 : #if (GRID_DO_COLLOCATE)
668 : // collocate
669 6519335504 : general_cij_to_ci(lp, dj, cij, ci);
670 3259667752 : general_ci_to_grid(lp, jg, kg, ismin, ismax, npts_local, index_min, index_max,
671 : map_i, sections_i, gp, k, j, exp_ij, exp_jk, exp_ik, ci,
672 : grid);
673 : #else
674 : // integrate
675 10764555173 : general_ci_to_grid(lp, jg, kg, ismin, ismax, npts_local, index_min, index_max,
676 : map_i, sections_i, gp, k, j, exp_ij, exp_jk, exp_ik, ci,
677 : grid);
678 : general_cij_to_ci(lp, dj, cij, ci);
679 : #endif
680 : }
681 :
682 : /*******************************************************************************
683 : * \brief Collocates coefficients C_ij onto the grid for general case.
684 : * \author Ole Schuett
685 : ******************************************************************************/
686 312771970 : static inline void __attribute__((always_inline)) general_cij_to_grid(
687 : const int lp, const int k, const int kg, const int npts_local[3],
688 : const int index_min[3], const int index_max[3], const int map_i[],
689 : const int map_j[], const int sections_i[], const int sections_j[],
690 : const double dh[3][3], const double gp[3], const double radius,
691 : const double exp_ij[], const double exp_jk[], const double exp_ik[],
692 : GRID_CONST_WHEN_COLLOCATE double *cij,
693 : GRID_CONST_WHEN_INTEGRATE double *grid) {
694 :
695 312771970 : const double a =
696 312771970 : dh[0][0] * dh[0][0] + dh[0][1] * dh[0][1] + dh[0][2] * dh[0][2];
697 312771970 : const double inv_2a = 1.0 / (2.0 * a);
698 312771970 : const double radius2 = radius * radius;
699 312771970 : const double gp0_dh0[3] = {(0 - gp[0]) * dh[0][0], (0 - gp[0]) * dh[0][1],
700 312771970 : (0 - gp[0]) * dh[0][2]};
701 312771970 : const double dk_dh2[3] = {(k - gp[2]) * dh[2][0], (k - gp[2]) * dh[2][1],
702 312771970 : (k - gp[2]) * dh[2][2]};
703 :
704 9161576994 : for (int j = index_min[1]; j <= index_max[1]; j++) {
705 8848805024 : const int jg = map_j[j - index_min[1]];
706 8848805024 : if (jg < 0) {
707 0 : j += sections_j[j - index_min[1]]; // skip over out-of-bounds indicies
708 0 : continue;
709 : }
710 :
711 : //--------------------------------------------------------------------
712 : // Find bounds for the inner loop based on a quadratic equation in i.
713 : //
714 : // The real-space vector from the center of the gaussian to the
715 : // grid point i,j,k is given by:
716 : // r = (i-gp[0])*dh[0,:] + (j-gp[1])*dh[1,:] + (k-gp[2])*dh[2,:]
717 : //
718 : // Separating the term that depends on i:
719 : // r = i*dh[0,:] - gp[0]*dh[0,:] + (j-gp[1])*dh[1,:] + (k-gp[2])*dh[2,:]
720 : // = i*dh[0,:] + v
721 : //
722 : // The squared distance works out to:
723 : // r**2 = dh[0,:]**2 * i**2 + 2 * v * dh[0,:] * i + v**2
724 : // = a * i**2 + b * i + c
725 : //
726 : // Solving r**2==radius**2 for i yields:
727 : // d = b**2 - 4 * a * (c - radius**2)
728 : // i = (-b \pm sqrt(d)) / (2*a)
729 : //
730 8848805024 : const double dj = j - gp[1];
731 8848805024 : double b = 0.0, c = 0.0;
732 35395220096 : for (int i = 0; i < 3; i++) {
733 26546415072 : const double v = (gp0_dh0[i] + dj * dh[1][i]) + dk_dh2[i];
734 26546415072 : b += 2.0 * v * dh[0][i];
735 26546415072 : c += v * v;
736 : }
737 8848805024 : const double d = b * b - 4.0 * a * (c - radius2);
738 :
739 8848805024 : if (0.0 < d) {
740 5859092566 : const double sqrt_d = sqrt(d);
741 5859092566 : const int ismin = (int)ceil((-b - sqrt_d) * inv_2a);
742 5859092566 : const int ismax = (int)floor((-b + sqrt_d) * inv_2a);
743 :
744 5859092566 : double ci[lp + 1];
745 5859092566 : memset(ci, 0, sizeof(ci));
746 :
747 : // Separate branches for low values of lp are generated in
748 : // general_cijk_to_grid.
749 11718185132 : general_cij_to_grid_low(lp, jg, kg, ismin, ismax, npts_local, index_min,
750 : index_max, map_i, sections_i, gp, k, j, exp_ij,
751 : exp_jk, exp_ik, dj, ci, cij, grid);
752 : }
753 : }
754 173276029 : }
755 :
756 : /*******************************************************************************
757 : * \brief Transforms coefficients C_ijk into C_ij by fixing grid index k.
758 : * \author Ole Schuett
759 : ******************************************************************************/
760 312771970 : static inline void general_cijk_to_cij(const int lp, const double dk,
761 : GRID_CONST_WHEN_COLLOCATE double *cijk,
762 : GRID_CONST_WHEN_INTEGRATE double *cij) {
763 312771970 : double dkp = 1.0;
764 1283962679 : for (int kl = 0; kl <= lp; kl++) {
765 3230227009 : for (int jl = 0; jl <= lp - kl; jl++) {
766 6769013057 : for (int il = 0; il <= lp - kl - jl; il++) {
767 4509976757 : const int cij_index = jl * (lp + 1) + il; // [jl, il]
768 4509976757 : const int cijk_index =
769 4509976757 : kl * (lp + 1) * (lp + 1) + jl * (lp + 1) + il; // [kl, jl, il]
770 : #if (GRID_DO_COLLOCATE)
771 2431053971 : cij[cij_index] += cijk[cijk_index] * dkp; // collocate
772 : #else
773 2078922786 : cijk[cijk_index] += cij[cij_index] * dkp; // integrate
774 : #endif
775 : }
776 : }
777 971190709 : dkp *= dk;
778 : }
779 312771970 : }
780 :
781 : /*******************************************************************************
782 : * \brief Precompute mapping of grid indices and its homogeneous sections.
783 : * Returns false if no index in [index_min, index_max] maps into bounds.
784 : * \author Ole Schuett
785 : ******************************************************************************/
786 : static inline bool
787 39011790 : general_precompute_mapping(const int index_min, const int index_max,
788 : const int shift_local, const int npts_global,
789 : const int bounds[2], int map[], int sections[]) {
790 :
791 : // Precompute mapping from continous grid indices to pbc wraped.
792 39011790 : bool has_valid = false;
793 1017394642 : for (int k = index_min; k <= index_max; k++) {
794 978382852 : const int kg = modulo(k - shift_local, npts_global);
795 978382852 : if (bounds[0] <= kg && kg <= bounds[1]) {
796 978377645 : map[k - index_min] = kg;
797 978377645 : has_valid = true;
798 : } else {
799 5207 : map[k - index_min] = INT_MIN; // out of bounds - not mapped
800 : }
801 : }
802 :
803 39011790 : if (!has_valid) {
804 : return false;
805 : }
806 :
807 : // Precompute length of sections with homogeneous cube to grid mapping.
808 39011740 : const int range = index_max - index_min + 1;
809 1017393718 : for (int kg = range - 1; kg >= 0; kg--) {
810 978381978 : if (kg == range - 1 || map[kg] != map[kg + 1] - 1) {
811 98693701 : sections[kg] = 0;
812 : } else {
813 879688277 : sections[kg] = sections[kg + 1] + 1;
814 : }
815 : }
816 : return true;
817 : }
818 :
819 : /*******************************************************************************
820 : * \brief Fill one of the 2D tables that speedup 3D Gaussian (Mathieu's trick).
821 : * \author Ole Schuett
822 : ******************************************************************************/
823 : static inline void
824 39011643 : general_fill_exp_table(const int idir, const int jdir, const int index_min[3],
825 : const int index_max[3], const double zetp,
826 : const double dh[3][3], const double gp[3],
827 39011643 : double exp_table[]) {
828 :
829 39011643 : const int stride_i = index_max[idir] - index_min[idir] + 1;
830 39011643 : const double h_ii = dh[idir][0] * dh[idir][0] + dh[idir][1] * dh[idir][1] +
831 39011643 : dh[idir][2] * dh[idir][2];
832 39011643 : const double h_jj = dh[jdir][0] * dh[jdir][0] + dh[jdir][1] * dh[jdir][1] +
833 39011643 : dh[jdir][2] * dh[jdir][2];
834 39011643 : const double h_ij = dh[idir][0] * dh[jdir][0] + dh[idir][1] * dh[jdir][1] +
835 39011643 : dh[idir][2] * dh[jdir][2];
836 :
837 39011643 : double exp_half_i[stride_i];
838 1038571661 : for (int i = index_min[idir]; i <= index_max[idir]; i++) {
839 999560018 : const double di = i - gp[idir];
840 999560018 : exp_half_i[i - index_min[idir]] = exp(-0.5 * zetp * di * di * h_ii);
841 : }
842 :
843 996212127 : for (int j = index_min[jdir]; j <= index_max[jdir]; j++) {
844 957200484 : const double dj = j - gp[jdir];
845 957200484 : const double rjj = 0.5 * dj * dj * h_jj;
846 957200484 : const double rij_unit = dj * h_ij;
847 957200484 : const double exp_ij_unit = exp(-zetp * 2.0 * rij_unit);
848 :
849 : // compute exponentials symmetrically around cube center
850 957200484 : const int i_center = (int)gp[idir];
851 957200484 : const double di_center = i_center - gp[idir];
852 957200484 : const double rij_center = di_center * rij_unit;
853 957200484 : const double exp_ij_center = exp(-zetp * (rjj + 2.0 * rij_center));
854 :
855 : // above center
856 957200484 : double exp_ij = exp_ij_center;
857 15076790989 : for (int i = i_center; i <= index_max[idir]; i++) {
858 14119590505 : const int idx = (j - index_min[jdir]) * stride_i + i - index_min[idir];
859 14119590505 : exp_table[idx] = exp_ij * exp_half_i[i - index_min[idir]];
860 14119590505 : exp_ij *= exp_ij_unit;
861 : }
862 :
863 : // below center
864 957200484 : const double exp_ij_unit_inv = 1.0 / exp_ij_unit;
865 957200484 : exp_ij = exp_ij_center * exp_ij_unit_inv;
866 14050863775 : for (int i = i_center - 1; i >= index_min[idir]; i--) {
867 13093663291 : const int idx = (j - index_min[jdir]) * stride_i + i - index_min[idir];
868 13093663291 : exp_table[idx] = exp_ij * exp_half_i[i - index_min[idir]];
869 13093663291 : exp_ij *= exp_ij_unit_inv;
870 : }
871 : }
872 39011643 : }
873 :
874 : /*******************************************************************************
875 : * \brief K-loop body of general_cijk_to_grid to be inlined for low values of
876 : * lp.
877 : ******************************************************************************/
878 13003881 : static inline void __attribute__((always_inline)) general_cijk_kloop_low(
879 : const int lp, const int npts_local[3], const int index_min[3],
880 : const int index_max[3], const int map_i[], const int map_j[],
881 : const int map_k[], const int sections_i[], const int sections_j[],
882 : const int sections_k[], const double dh[3][3], const double gp[3],
883 : const double radius, const double exp_ij[], const double exp_jk[],
884 : const double exp_ik[], double *cij, GRID_CONST_WHEN_COLLOCATE double *cijk,
885 : GRID_CONST_WHEN_INTEGRATE double *grid) {
886 :
887 13003881 : const int cij_size = (lp + 1) * (lp + 1);
888 325780184 : for (int k = index_min[2]; k <= index_max[2]; k++) {
889 312776303 : const int kg = map_k[k - index_min[2]];
890 312776303 : if (kg < 0) {
891 4333 : k += sections_k[k - index_min[2]]; // skip over out-of-bounds indicies
892 4333 : continue;
893 : }
894 :
895 : // zero coef_xyt
896 312771970 : memset(cij, 0, cij_size * sizeof(double));
897 :
898 : #if (GRID_DO_COLLOCATE)
899 : // collocate
900 173276029 : general_cijk_to_cij(lp, (double)k - gp[2], cijk, cij);
901 346552058 : general_cij_to_grid(lp, k, kg, npts_local, index_min, index_max, map_i,
902 : map_j, sections_i, sections_j, dh, gp, radius, exp_ij,
903 : exp_jk, exp_ik, cij, grid);
904 : #else
905 : // integrate
906 139495941 : general_cij_to_grid(lp, k, kg, npts_local, index_min, index_max, map_i,
907 : map_j, sections_i, sections_j, dh, gp, radius, exp_ij,
908 : exp_jk, exp_ik, cij, grid);
909 139495941 : general_cijk_to_cij(lp, (double)k - gp[2], cijk, cij);
910 : #endif
911 : }
912 : }
913 :
914 : /*******************************************************************************
915 : * \brief Collocates coefficients C_ijk onto the grid for general case.
916 : * \author Ole Schuett
917 : ******************************************************************************/
918 : // Keep out-of-line (called once per task) so GCC optimizes the 11 unrolled
919 : // general k-loop lp branches as a standalone function and does not exhaust
920 : // the inlining budget in grid_cpu_*_pgf_product (ortho_cxy_to_grid).
921 : static void __attribute__((noinline))
922 13003931 : general_cijk_to_grid(const int border_mask, const int lp, const double zetp,
923 : const double dh[3][3], const double dh_inv[3][3],
924 : const double rp[3], const int npts_global[3],
925 : const int npts_local[3], const int shift_local[3],
926 : const int border_width[3], const double radius,
927 : GRID_CONST_WHEN_COLLOCATE double *cijk,
928 13003931 : GRID_CONST_WHEN_INTEGRATE double *grid) {
929 :
930 : // Default for border_mask == 0.
931 13003931 : int bounds_i[2] = {0, npts_local[0] - 1};
932 13003931 : int bounds_j[2] = {0, npts_local[1] - 1};
933 13003931 : int bounds_k[2] = {0, npts_local[2] - 1};
934 :
935 : // See also rs_find_node() in task_list_methods.F.
936 : // If the bit is set then we need to exclude the border in that direction.
937 13003931 : if (border_mask & (1 << 0))
938 152246 : bounds_i[0] += border_width[0];
939 13003931 : if (border_mask & (1 << 1))
940 152246 : bounds_i[1] -= border_width[0];
941 13003931 : if (border_mask & (1 << 2))
942 152246 : bounds_j[0] += border_width[1];
943 13003931 : if (border_mask & (1 << 3))
944 152246 : bounds_j[1] -= border_width[1];
945 13003931 : if (border_mask & (1 << 4))
946 152238 : bounds_k[0] += border_width[2];
947 13003931 : if (border_mask & (1 << 5))
948 152246 : bounds_k[1] -= border_width[2];
949 :
950 13003931 : if (bounds_i[0] > bounds_i[1] || bounds_j[0] > bounds_j[1] ||
951 13003931 : bounds_k[0] > bounds_k[1]) {
952 50 : return;
953 : }
954 :
955 : // center in grid coords
956 : // gp = MATMUL(dh_inv, rp)
957 13003931 : double gp[3] = {0.0, 0.0, 0.0};
958 52015724 : for (int i = 0; i < 3; i++) {
959 156047172 : for (int j = 0; j < 3; j++) {
960 117035379 : gp[i] += dh_inv[j][i] * rp[j];
961 : }
962 : }
963 :
964 : // Get the min max indices that contain the sphere around rp of radius radius.
965 : // Project the sphere onto each grid axis via the reciprocal-lattice vector
966 : // norms ||dh_inv[:,idir]||_2 so non-orthogonal cells avoid excess box
967 : // corners.
968 : // Guard against floating-point rounding in sqrt/norm when gp +/- rad_idx hits
969 : // an int.
970 : const double sphere_bound_eps = 1.0 + 1e-12;
971 : int index_min[3];
972 : int index_max[3];
973 52015724 : for (int idir = 0; idir < 3; idir++) {
974 39011793 : const double norm = sqrt(dh_inv[0][idir] * dh_inv[0][idir] +
975 39011793 : dh_inv[1][idir] * dh_inv[1][idir] +
976 39011793 : dh_inv[2][idir] * dh_inv[2][idir]);
977 39011793 : const double rad_idx = radius * norm * sphere_bound_eps;
978 39011793 : index_min[idir] = (int)floor(gp[idir] - rad_idx);
979 39011793 : index_max[idir] = (int)ceil(gp[idir] + rad_idx);
980 : }
981 :
982 : // Mappings and exponential tables live in thread local heap scratch instead
983 : // of VLAs, because large radii on fine grids can exceed the default stack.
984 13003931 : const int range_i = index_max[0] - index_min[0] + 1;
985 13003931 : const int range_j = index_max[1] - index_min[1] + 1;
986 13003931 : const int range_k = index_max[2] - index_min[2] + 1;
987 13003931 : const size_t exp_count = (size_t)range_i * range_j +
988 : (size_t)range_j * range_k +
989 13003931 : (size_t)range_k * range_i;
990 13003931 : const size_t map_count = 2 * ((size_t)range_i + range_j + range_k);
991 26007862 : double *const scratch = (double *)grid_library_get_scratch(
992 13003931 : exp_count * sizeof(double) + map_count * sizeof(int));
993 13003931 : double *restrict const exp_ij = scratch;
994 13003931 : double *restrict const exp_jk = exp_ij + (size_t)range_i * range_j;
995 13003931 : double *restrict const exp_ik = exp_jk + (size_t)range_j * range_k;
996 13003931 : int *restrict const map_i = (int *)(exp_ik + (size_t)range_k * range_i);
997 13003931 : int *restrict const sections_i = map_i + range_i;
998 13003931 : int *restrict const map_j = sections_i + range_i;
999 13003931 : int *restrict const sections_j = map_j + range_j;
1000 13003931 : int *restrict const map_k = sections_j + range_j;
1001 13003931 : int *restrict const sections_k = map_k + range_k;
1002 13003931 : if (!general_precompute_mapping(index_min[0], index_max[0], shift_local[0],
1003 : npts_global[0], bounds_i, map_i,
1004 : sections_i)) {
1005 : return;
1006 : }
1007 13003931 : if (!general_precompute_mapping(index_min[1], index_max[1], shift_local[1],
1008 : npts_global[1], bounds_j, map_j,
1009 : sections_j)) {
1010 : return;
1011 : }
1012 13003928 : if (!general_precompute_mapping(index_min[2], index_max[2], shift_local[2],
1013 : npts_global[2], bounds_k, map_k,
1014 : sections_k)) {
1015 : return;
1016 : }
1017 :
1018 : // Precompute exponentials
1019 13003881 : general_fill_exp_table(0, 1, index_min, index_max, zetp, dh, gp, exp_ij);
1020 13003881 : general_fill_exp_table(1, 2, index_min, index_max, zetp, dh, gp, exp_jk);
1021 13003881 : general_fill_exp_table(0, 2, index_min, index_max, zetp, dh, gp, exp_ik);
1022 :
1023 : // go over the grid, but cycle if the point is not within the radius
1024 13003881 : const int cij_size = (lp + 1) * (lp + 1);
1025 13003881 : double cij[cij_size];
1026 : // Generate separate branches for low values of lp.
1027 13003881 : if (lp <= GRID_MAX_LP_OPTIMIZED) {
1028 : GRID_PRAGMA_UNROLL(GRID_MAX_LP_OPTIMIZED + 1)
1029 143042691 : for (int ilp = 0; ilp <= GRID_MAX_LP_OPTIMIZED; ilp++) {
1030 130038810 : if (lp == ilp) {
1031 143042691 : general_cijk_kloop_low(ilp, npts_local, index_min, index_max, map_i,
1032 : map_j, map_k, sections_i, sections_j, sections_k,
1033 : dh, gp, radius, exp_ij, exp_jk, exp_ik, cij,
1034 : cijk, grid);
1035 : }
1036 : }
1037 : } else {
1038 13003881 : general_cijk_kloop_low(lp, npts_local, index_min, index_max, map_i, map_j,
1039 : map_k, sections_i, sections_j, sections_k, dh, gp,
1040 : radius, exp_ij, exp_jk, exp_ik, cij, cijk, grid);
1041 : }
1042 : }
1043 :
1044 : /*******************************************************************************
1045 : * \brief Transforms coefficients C_xyz into C_ijk.
1046 : * \author Ole Schuett
1047 : ******************************************************************************/
1048 : static inline void
1049 13003931 : general_cxyz_to_cijk(const int lp, const double dh[3][3],
1050 : GRID_CONST_WHEN_COLLOCATE double *cxyz,
1051 13003931 : GRID_CONST_WHEN_INTEGRATE double *cijk) {
1052 :
1053 : // transform P_{lxp,lyp,lzp} into a P_{lip,ljp,lkp} such that
1054 : // sum_{lxp,lyp,lzp} P_{lxp,lyp,lzp} (x-x_p)**lxp (y-y_p)**lyp (z-z_p)**lzp =
1055 : // sum_{lip,ljp,lkp} P_{lip,ljp,lkp} (i-i_p)**lip (j-j_p)**ljp (k-k_p)**lkp
1056 :
1057 : // transform using multinomials
1058 13003931 : const int lp1 = lp + 1;
1059 13003931 : double hmatgridp[lp1][3][3];
1060 52015724 : for (int i = 0; i < 3; i++) {
1061 156047172 : for (int j = 0; j < 3; j++) {
1062 117035379 : hmatgridp[0][j][i] = 1.0;
1063 358479945 : for (int k = 1; k <= lp; k++) {
1064 241444566 : hmatgridp[k][j][i] = hmatgridp[k - 1][j][i] * dh[j][i];
1065 : }
1066 : }
1067 : }
1068 :
1069 13003931 : double trinom[3][lp1][lp1][lp1];
1070 52015724 : for (int d = 0; d < 3; d++) {
1071 158505108 : for (int kl = 0; kl <= lp; kl++) {
1072 396722511 : for (int jl = 0; jl <= lp - kl; jl++) {
1073 831024816 : for (int il = 0; il <= lp - kl - jl; il++) {
1074 553795620 : const int l = il + jl + kl;
1075 553795620 : trinom[d][kl][jl][il] = hmatgridp[il][0][d] * hmatgridp[jl][1][d] *
1076 553795620 : hmatgridp[kl][2][d] * fac(l) /
1077 553795620 : (fac(il) * fac(jl) * fac(kl));
1078 : }
1079 : }
1080 : }
1081 : }
1082 :
1083 : const int lpx = lp;
1084 52835036 : for (int klx = 0; klx <= lpx; klx++) {
1085 132240837 : for (int jlx = 0; jlx <= lpx - klx; jlx++) {
1086 277008272 : for (int ilx = 0; ilx <= lpx - klx - jlx; ilx++) {
1087 184598540 : const int lx = ilx + jlx + klx;
1088 184598540 : const double px = trinom[0][klx][jlx][ilx];
1089 184598540 : const int lpy = lp - lx;
1090 520374971 : for (int kly = 0; kly <= lpy; kly++) {
1091 908246703 : for (int jly = 0; jly <= lpy - kly; jly++) {
1092 1502860255 : for (int ily = 0; ily <= lpy - kly - jly; ily++) {
1093 930389983 : const int ly = ily + jly + kly;
1094 930389983 : const double pxy = px * trinom[1][kly][jly][ily];
1095 930389983 : const int il_xy = ilx + ily;
1096 930389983 : const int jl_xy = jlx + jly;
1097 930389983 : const int kl_xy = klx + kly;
1098 930389983 : const int lpz = lp - lx - ly;
1099 930389983 : const int cxyz_base = ly * lp1 + lx;
1100 2387310900 : for (int klz = 0; klz <= lpz; klz++) {
1101 1456920917 : const int kl_base = (kl_xy + klz) * lp1 * lp1;
1102 3671045581 : for (int jlz = 0; jlz <= lpz - klz; jlz++) {
1103 2214124664 : const int cijk_base = kl_base + (jl_xy + jlz) * lp1 + il_xy;
1104 5496425277 : for (int ilz = 0; ilz <= lpz - klz - jlz; ilz++) {
1105 3282300613 : const int lz = ilz + jlz + klz;
1106 3282300613 : const int cijk_index = cijk_base + ilz; // [kl,jl,il]
1107 3282300613 : const int cxyz_index =
1108 3282300613 : lz * lp1 * lp1 + cxyz_base; // [lz,ly,lx]
1109 3282300613 : const double p = pxy * trinom[2][klz][jlz][ilz];
1110 : #if (GRID_DO_COLLOCATE)
1111 1442571749 : cijk[cijk_index] += cxyz[cxyz_index] * p; // collocate
1112 : #else
1113 1839728864 : cxyz[cxyz_index] += cijk[cijk_index] * p; // integrate
1114 : #endif
1115 : }
1116 : }
1117 : }
1118 : }
1119 : }
1120 : }
1121 : }
1122 : }
1123 : }
1124 13003931 : }
1125 :
1126 : /*******************************************************************************
1127 : * \brief Collocates coefficients C_xyz onto the grid for general case.
1128 : * \author Ole Schuett
1129 : ******************************************************************************/
1130 : static inline void
1131 13003931 : general_cxyz_to_grid(const int border_mask, const int lp, const double zetp,
1132 : const double dh[3][3], const double dh_inv[3][3],
1133 : const double rp[3], const int npts_global[3],
1134 : const int npts_local[3], const int shift_local[3],
1135 : const int border_width[3], const double radius,
1136 : GRID_CONST_WHEN_COLLOCATE double *cxyz,
1137 13003931 : GRID_CONST_WHEN_INTEGRATE double *grid) {
1138 :
1139 13003931 : const size_t cijk_size = (lp + 1) * (lp + 1) * (lp + 1);
1140 13003931 : double cijk[cijk_size];
1141 13003931 : memset(cijk, 0, cijk_size * sizeof(double));
1142 :
1143 : #if (GRID_DO_COLLOCATE)
1144 : // collocate
1145 7139360 : general_cxyz_to_cijk(lp, dh, cxyz, cijk);
1146 7139360 : general_cijk_to_grid(border_mask, lp, zetp, dh, dh_inv, rp, npts_global,
1147 : npts_local, shift_local, border_width, radius, cijk,
1148 : grid);
1149 : #else
1150 : // integrate
1151 5864571 : general_cijk_to_grid(border_mask, lp, zetp, dh, dh_inv, rp, npts_global,
1152 : npts_local, shift_local, border_width, radius, cijk,
1153 : grid);
1154 5864571 : general_cxyz_to_cijk(lp, dh, cxyz, cijk);
1155 : #endif
1156 13003931 : }
1157 :
1158 : /*******************************************************************************
1159 : * \brief Collocates coefficients C_xyz onto the grid.
1160 : * \author Ole Schuett
1161 : ******************************************************************************/
1162 : static inline void
1163 194661090 : cxyz_to_grid(const bool orthorhombic, const int border_mask, const int lp,
1164 : const double zetp, const double dh[3][3],
1165 : const double dh_inv[3][3], const double rp[3],
1166 : const int npts_global[3], const int npts_local[3],
1167 : const int shift_local[3], const int border_width[3],
1168 : const double radius, GRID_CONST_WHEN_COLLOCATE double *cxyz,
1169 : GRID_CONST_WHEN_INTEGRATE double *grid) {
1170 :
1171 194661090 : enum grid_library_kernel k;
1172 194661090 : if (orthorhombic && border_mask == 0) {
1173 181657159 : k = (GRID_DO_COLLOCATE) ? GRID_COLLOCATE_ORTHO : GRID_INTEGRATE_ORTHO;
1174 181657159 : ortho_cxyz_to_grid(lp, zetp, dh, dh_inv, rp, npts_global, npts_local,
1175 : shift_local, radius, cxyz, grid);
1176 : } else {
1177 13003931 : k = (GRID_DO_COLLOCATE) ? GRID_COLLOCATE_GENERAL : GRID_INTEGRATE_GENERAL;
1178 13003931 : general_cxyz_to_grid(border_mask, lp, zetp, dh, dh_inv, rp, npts_global,
1179 : npts_local, shift_local, border_width, radius, cxyz,
1180 : grid);
1181 : }
1182 194661090 : grid_library_counter_add(lp, GRID_BACKEND_CPU, k, 1);
1183 194661090 : }
1184 :
1185 : /*******************************************************************************
1186 : * \brief Transforms coefficients C_ab into C_xyz.
1187 : * \author Ole Schuett
1188 : ******************************************************************************/
1189 194661090 : static inline void cab_to_cxyz(const int la_max, const int la_min,
1190 : const int lb_max, const int lb_min,
1191 : const double prefactor, const double ra[3],
1192 : const double rb[3], const double rp[3],
1193 : GRID_CONST_WHEN_COLLOCATE double *cab,
1194 194661090 : GRID_CONST_WHEN_INTEGRATE double *cxyz) {
1195 :
1196 : // Computes the polynomial expansion coefficients:
1197 : // (x-a)**lxa (x-b)**lxb -> sum_{ls} alpha(ls,lxa,lxb,1)*(x-p)**ls
1198 194661090 : const int lp = la_max + lb_max;
1199 194661090 : const int lp1 = lp + 1;
1200 194661090 : double alpha[3][lb_max + 1][la_max + 1][lp + 1];
1201 :
1202 778644360 : for (int i = 0; i < 3; i++) {
1203 583983270 : const double drpa = rp[i] - ra[i];
1204 583983270 : const double drpb = rp[i] - rb[i];
1205 583983270 : alpha[i][0][0][0] = 1.0;
1206 1098338715 : for (int lxa = 1; lxa <= la_max; lxa++) {
1207 514355445 : alpha[i][0][lxa][0] = drpa * alpha[i][0][lxa - 1][0];
1208 627657639 : for (int k = 1; k < lxa; k++) {
1209 113302194 : alpha[i][0][lxa][k] =
1210 113302194 : alpha[i][0][lxa - 1][k - 1] + drpa * alpha[i][0][lxa - 1][k];
1211 : }
1212 514355445 : alpha[i][0][lxa][lxa] = 1.0;
1213 : }
1214 1062874938 : for (int lxb = 1; lxb <= lb_max; lxb++) {
1215 1519570572 : for (int lxa = 0; lxa <= la_max; lxa++) {
1216 1040678904 : const int n = lxa + lxb;
1217 1040678904 : alpha[i][lxb][lxa][0] = drpb * alpha[i][lxb - 1][lxa][0];
1218 1982142534 : for (int k = 1; k < n; k++) {
1219 941463630 : alpha[i][lxb][lxa][k] =
1220 941463630 : alpha[i][lxb - 1][lxa][k - 1] + drpb * alpha[i][lxb - 1][lxa][k];
1221 : }
1222 1040678904 : alpha[i][lxb][lxa][n] = 1.0;
1223 : }
1224 : }
1225 : }
1226 :
1227 : // *** initialise the coefficient matrix, we transform the sum
1228 : //
1229 : // sum_{lxa,lya,lza,lxb,lyb,lzb} P_{lxa,lya,lza,lxb,lyb,lzb} *
1230 : // (x-a_x)**lxa (y-a_y)**lya (z-a_z)**lza (x-b_x)**lxb (y-a_y)**lya
1231 : // (z-a_z)**lza
1232 : //
1233 : // into
1234 : //
1235 : // sum_{lxp,lyp,lzp} P_{lxp,lyp,lzp} (x-p_x)**lxp (y-p_y)**lyp (z-p_z)**lzp
1236 : //
1237 : // where p is center of the product gaussian, and lp = la_max + lb_max
1238 : // (current implementation is l**7)
1239 : //
1240 : // Factor the 3D transformation (lxa,lya,lza,lxb,lyb,lzb) <-> (lxp,lyp,lzp)
1241 : // into three separable 1D steps along x, y, z.
1242 194661090 : const int ncoset_a = ncoset(la_max);
1243 194661090 : double txy[lp1 * lp1];
1244 194661090 : double tx[lp1];
1245 :
1246 : #if (GRID_DO_COLLOCATE)
1247 280278531 : for (int lzb = 0; lzb <= lb_max; lzb++) {
1248 538641639 : for (int lza = 0; lza <= la_max; lza++) {
1249 357161669 : const int lz_max = lza + lzb;
1250 357161669 : const int lyp_max = lp - lz_max;
1251 357161669 : memset(txy, 0, (lyp_max + 1) * lp1 * sizeof(double));
1252 929734709 : for (int lyb = 0; lyb <= lb_max - lzb; lyb++) {
1253 1509673391 : for (int lya = 0; lya <= la_max - lza; lya++) {
1254 937100351 : const int ly_max = lya + lyb;
1255 937100351 : const int lxb_min = imax(lb_min - lzb - lyb, 0);
1256 937100351 : const int lxa_min = imax(la_min - lza - lya, 0);
1257 937100351 : const int lxb_max = lb_max - lzb - lyb;
1258 937100351 : const int lxa_max = la_max - lza - lya;
1259 937100351 : const int lxp_max = lxa_max + lxb_max;
1260 937100351 : memset(tx, 0, (lxp_max + 1) * sizeof(double));
1261 2232244544 : for (int lxb = lxb_min; lxb <= lxb_max; lxb++) {
1262 1295144193 : const int jco = coset(lxb, lyb, lzb);
1263 3127928710 : for (int lxa = lxa_min; lxa <= lxa_max; lxa++) {
1264 1832784517 : const int ico = coset(lxa, lya, lza);
1265 1832784517 : const int cab_index = jco * ncoset_a + ico; // [jco, ico]
1266 1832784517 : const double cab_val = cab[cab_index] * prefactor;
1267 5178225963 : for (int lxp = 0; lxp <= lxa + lxb; lxp++) {
1268 3345441446 : tx[lxp] += cab_val * alpha[0][lxb][lxa][lxp];
1269 : }
1270 : }
1271 : }
1272 2760865422 : for (int lyp = 0; lyp <= ly_max; lyp++) {
1273 1823765071 : const double ay = alpha[1][lyb][lya][lyp];
1274 1823765071 : double *txy_row = &txy[lyp * lp1];
1275 5176519038 : for (int lxp = 0; lxp <= lxp_max; lxp++) {
1276 3352753967 : txy_row[lxp] += tx[lxp] * ay;
1277 : }
1278 : }
1279 : }
1280 : }
1281 1143492037 : for (int lzp = 0; lzp <= lz_max; lzp++) {
1282 786330368 : const double az = alpha[2][lzb][lza][lzp];
1283 786330368 : double *cxyz_z = &cxyz[lzp * lp1 * lp1];
1284 2385713031 : for (int lyp = 0; lyp <= lyp_max; lyp++) {
1285 1599382663 : const double *txy_row = &txy[lyp * lp1];
1286 1599382663 : double *cxyz_zy = &cxyz_z[lyp * lp1];
1287 4397314240 : for (int lxp = 0; lxp <= lyp_max - lyp; lxp++) {
1288 2797931577 : cxyz_zy[lxp] += txy_row[lxp] * az;
1289 : }
1290 : }
1291 : }
1292 : }
1293 : }
1294 : #else
1295 268674205 : for (int lzb = 0; lzb <= lb_max; lzb++) {
1296 528655880 : for (int lza = 0; lza <= la_max; lza++) {
1297 355844204 : const int lz_max = lza + lzb;
1298 355844204 : const int lyp_max = lp - lz_max;
1299 355844204 : memset(txy, 0, (lyp_max + 1) * lp1 * sizeof(double));
1300 1152453100 : for (int lzp = 0; lzp <= lz_max; lzp++) {
1301 796608896 : const double az = alpha[2][lzb][lza][lzp] * prefactor;
1302 796608896 : const double *cxyz_z = &cxyz[lzp * lp1 * lp1];
1303 2444363176 : for (int lyp = 0; lyp <= lyp_max; lyp++) {
1304 1647754280 : double *txy_row = &txy[lyp * lp1];
1305 1647754280 : const double *cxyz_zy = &cxyz_z[lyp * lp1];
1306 4580397402 : for (int lxp = 0; lxp <= lyp_max - lyp; lxp++) {
1307 2932643122 : txy_row[lxp] += cxyz_zy[lxp] * az;
1308 : }
1309 : }
1310 : }
1311 916697220 : for (int lyb = 0; lyb <= lb_max - lzb; lyb++) {
1312 1514527172 : for (int lya = 0; lya <= la_max - lza; lya++) {
1313 953674156 : const int ly_max = lya + lyb;
1314 953674156 : const int lxb_min = imax(lb_min - lzb - lyb, 0);
1315 953674156 : const int lxa_min = imax(la_min - lza - lya, 0);
1316 953674156 : const int lxb_max = lb_max - lzb - lyb;
1317 953674156 : const int lxa_max = la_max - lza - lya;
1318 953674156 : const int lxp_max = lxa_max + lxb_max;
1319 953674156 : memset(tx, 0, (lxp_max + 1) * sizeof(double));
1320 2843748767 : for (int lyp = 0; lyp <= ly_max; lyp++) {
1321 1890074611 : const double ay = alpha[1][lyb][lya][lyp];
1322 1890074611 : const double *txy_row = &txy[lyp * lp1];
1323 5436714992 : for (int lxp = 0; lxp <= lxp_max; lxp++) {
1324 3546640381 : tx[lxp] += txy_row[lxp] * ay;
1325 : }
1326 : }
1327 2267932831 : for (int lxb = lxb_min; lxb <= lxb_max; lxb++) {
1328 1314258675 : const int jco = coset(lxb, lyb, lzb);
1329 3273801815 : for (int lxa = lxa_min; lxa <= lxa_max; lxa++) {
1330 1959543140 : const int ico = coset(lxa, lya, lza);
1331 1959543140 : const int cab_index = jco * ncoset_a + ico; // [jco, ico]
1332 1959543140 : double sum_x = 0.0;
1333 5608616460 : for (int lxp = 0; lxp <= lxa + lxb; lxp++) {
1334 3649073320 : sum_x += tx[lxp] * alpha[0][lxb][lxa][lxp];
1335 : }
1336 1959543140 : cab[cab_index] += sum_x;
1337 : }
1338 : }
1339 : }
1340 : }
1341 : }
1342 : }
1343 : #endif
1344 194661090 : }
1345 :
1346 : /*******************************************************************************
1347 : * \brief Collocates coefficients C_ab onto the grid.
1348 : * \author Ole Schuett
1349 : ******************************************************************************/
1350 : static inline void
1351 244687165 : cab_to_grid(const bool orthorhombic, const int border_mask, const int la_max,
1352 : const int la_min, const int lb_max, const int lb_min,
1353 : const double zeta, const double zetb, const double rscale,
1354 : const double dh[3][3], const double dh_inv[3][3],
1355 : const double ra[3], const double rab[3], const int npts_global[3],
1356 : const int npts_local[3], const int shift_local[3],
1357 : const int border_width[3], const double radius,
1358 : GRID_CONST_WHEN_COLLOCATE double *cab,
1359 244687165 : GRID_CONST_WHEN_INTEGRATE double *grid) {
1360 :
1361 : // Check if radius is too small to be mapped onto grid of given resolution.
1362 244687165 : double dh_max = 0.0;
1363 978748660 : for (int i = 0; i < 3; i++) {
1364 2936245980 : for (int j = 0; j < 3; j++) {
1365 2202184485 : dh_max = fmax(dh_max, fabs(dh[i][j]));
1366 : }
1367 : }
1368 244687165 : if (2.0 * radius < dh_max) {
1369 50026075 : return;
1370 : }
1371 :
1372 194661090 : const double zetp = zeta + zetb;
1373 194661090 : const double f = zetb / zetp;
1374 194661090 : const double rab2 = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
1375 194661090 : const double prefactor = rscale * exp(-zeta * f * rab2);
1376 194661090 : double rp[3], rb[3];
1377 778644360 : for (int i = 0; i < 3; i++) {
1378 583983270 : rp[i] = ra[i] + f * rab[i];
1379 583983270 : rb[i] = ra[i] + rab[i];
1380 : }
1381 :
1382 194661090 : const int lp = la_max + lb_max;
1383 194661090 : const size_t cxyz_size = (lp + 1) * (lp + 1) * (lp + 1);
1384 194661090 : double cxyz[cxyz_size];
1385 194661090 : memset(cxyz, 0, cxyz_size * sizeof(double));
1386 :
1387 : #if (GRID_DO_COLLOCATE)
1388 : // collocate
1389 98798561 : cab_to_cxyz(la_max, la_min, lb_max, lb_min, prefactor, ra, rb, rp, cab, cxyz);
1390 98798561 : cxyz_to_grid(orthorhombic, border_mask, lp, zetp, dh, dh_inv, rp, npts_global,
1391 : npts_local, shift_local, border_width, radius, cxyz, grid);
1392 : #else
1393 : // integrate
1394 95862529 : cxyz_to_grid(orthorhombic, border_mask, lp, zetp, dh, dh_inv, rp, npts_global,
1395 : npts_local, shift_local, border_width, radius, cxyz, grid);
1396 95862529 : cab_to_cxyz(la_max, la_min, lb_max, lb_min, prefactor, ra, rb, rp, cab, cxyz);
1397 : #endif
1398 : }
1399 :
1400 : // EOF
|
