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 >32568*10^7 : ortho_cx_to_grid_scalar(const int lp, const int cmax, const int i,
39 >32568*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 >17106*10^7 : double reg[4] = {0.0, 0.0, 0.0, 0.0};
49 >13685*10^8 : #pragma omp simd reduction(+ : reg)
50 : for (int lxp = 0; lxp <= lp; lxp++) {
51 >46959*10^7 : const double p = pol[0][lxp][i + cmax];
52 >46959*10^7 : reg[0] += cx[lxp * 4 + 0] * p;
53 >46959*10^7 : reg[1] += cx[lxp * 4 + 1] * p;
54 >46959*10^7 : reg[2] += cx[lxp * 4 + 2] * p;
55 >46959*10^7 : reg[3] += cx[lxp * 4 + 3] * p;
56 : }
57 >17106*10^7 : *grid_0 += reg[0];
58 >17106*10^7 : *grid_1 += reg[1];
59 >17106*10^7 : *grid_2 += reg[2];
60 >17106*10^7 : *grid_3 += reg[3];
61 :
62 : #else
63 : // integrate
64 >15461*10^7 : const double reg[4] = {*grid_0, *grid_1, *grid_2, *grid_3};
65 : #pragma omp simd
66 >15461*10^7 : for (int lxp = 0; lxp <= lp; lxp++) {
67 >43198*10^7 : const double p = pol[0][lxp][i + cmax];
68 >43198*10^7 : cx[lxp * 4 + 0] += reg[0] * p;
69 >43198*10^7 : cx[lxp * 4 + 1] += reg[1] * p;
70 >43198*10^7 : cx[lxp * 4 + 2] += reg[2] * p;
71 >43198*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 17475999847 : ortho_cx_to_grid(const int lp, const int kg1, const int kg2, const int jg1,
164 : const int jg2, const int cmax,
165 17475999847 : const double pol[3][lp + 1][2 * cmax + 1],
166 17475999847 : const int map[3][2 * cmax + 1],
167 17475999847 : 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 17475999847 : const int lb = *((*sphere_bounds_iter)++);
173 17475999847 : 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 50379393594 : for (int istart = lb; istart <= ub; istart++) {
179 32903393747 : const int istop = imin(ub, istart + sections[0][istart + cmax]);
180 32903393747 : const int cube2grid = map[0][istart + cmax] - istart;
181 :
182 32903393747 : const int stride = npts_local[1] * npts_local[0];
183 32903393747 : const int grid_index_0 = kg1 * stride + jg1 * npts_local[0];
184 32903393747 : const int grid_index_1 = kg2 * stride + jg1 * npts_local[0];
185 32903393747 : const int grid_index_2 = kg1 * stride + jg2 * npts_local[0];
186 32903393747 : const int grid_index_3 = kg2 * stride + jg2 * npts_local[0];
187 32903393747 : GRID_CONST_WHEN_INTEGRATE double *grid_base_0 = &grid[grid_index_0];
188 32903393747 : GRID_CONST_WHEN_INTEGRATE double *grid_base_1 = &grid[grid_index_1];
189 32903393747 : GRID_CONST_WHEN_INTEGRATE double *grid_base_2 = &grid[grid_index_2];
190 32903393747 : 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 >35858*10^7 : for (int i = istart; i <= istop; i++) {
206 >32568*10^7 : const int ig = i + cube2grid;
207 >32568*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 >32568*10^7 : &grid_base_3[ig]);
210 : }
211 32903393747 : 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 17475999847 : ortho_cxy_to_cx(const int lp, const int j1, const int j2, const int cmax,
221 17475999847 : 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 65840291708 : for (int lyp = 0; lyp <= lp; lyp++) {
226 >15165*10^7 : for (int lxp = 0; lxp <= lp - lyp; lxp++) {
227 >10329*10^7 : const double p1 = pol[1][lyp][j1 + cmax];
228 >10329*10^7 : const double p2 = pol[1][lyp][j2 + cmax];
229 >10329*10^7 : const int cxy_index = lyp * (lp + 1) * 2 + lxp * 2; // [lyp, lxp, 0]
230 :
231 : #if (GRID_DO_COLLOCATE)
232 : // collocate
233 53123537643 : cx[lxp * 4 + 0] += cxy[cxy_index + 0] * p1;
234 53123537643 : cx[lxp * 4 + 1] += cxy[cxy_index + 1] * p1;
235 53123537643 : cx[lxp * 4 + 2] += cxy[cxy_index + 0] * p2;
236 53123537643 : cx[lxp * 4 + 3] += cxy[cxy_index + 1] * p2;
237 : #else
238 : // integrate
239 50169940326 : cxy[cxy_index + 0] += cx[lxp * 4 + 0] * p1;
240 50169940326 : cxy[cxy_index + 1] += cx[lxp * 4 + 1] * p1;
241 50169940326 : cxy[cxy_index + 0] += cx[lxp * 4 + 2] * p2;
242 50169940326 : 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 17475999847 : 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 17475999847 : const double pol[3][lp + 1][2 * cmax + 1], const int map[3][2 * cmax + 1],
256 17475999847 : 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 18236637312 : ortho_cxy_to_cx(lp, j1, j2, cmax, pol, cxy, cx);
263 9118318656 : 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 16715362382 : ortho_cx_to_grid(lp, kg1, kg2, jg1, jg2, cmax, pol, map, sections, npts_local,
268 : sphere_bounds_iter, cx, grid);
269 8357681191 : 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 1920301306 : static inline void ortho_cxy_to_grid(
278 : const int lp, const int kg1, const int kg2, const int cmax,
279 1920301306 : const double pol[3][lp + 1][2 * cmax + 1], const int map[3][2 * cmax + 1],
280 1920301306 : const int sections[3][2 * cmax + 1], const int npts_local[3],
281 : int **sphere_bounds_iter, GRID_CONST_WHEN_COLLOCATE double *cxy,
282 1920301306 : 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 1920301306 : const int jstart = *((*sphere_bounds_iter)++);
290 1920301306 : const size_t cx_size = (lp + 1) * 4;
291 1920301306 : double cx[cx_size];
292 19396301153 : for (int j1 = jstart; j1 <= 0; j1++) {
293 17475999847 : const int j2 = 1 - j1;
294 17475999847 : const int jg1 = map[1][j1 + cmax];
295 17475999847 : const int jg2 = map[1][j2 + cmax];
296 :
297 17475999847 : memset(cx, 0, cx_size * sizeof(double));
298 :
299 : // Generate separate branches for low values of lp gives up to 30% speedup.
300 17475999847 : if (lp <= GRID_MAX_LP_OPTIMIZED) {
301 : GRID_PRAGMA_UNROLL(GRID_MAX_LP_OPTIMIZED + 1)
302 >19223*10^7 : for (int ilp = 0; ilp <= GRID_MAX_LP_OPTIMIZED; ilp++) {
303 >17475*10^7 : if (lp == ilp) {
304 >19223*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 17476003465 : 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 1920301306 : }
316 :
317 : /*******************************************************************************
318 : * \brief Transforms coefficients C_xyz into C_xz by fixing grid index k.
319 : * \author Ole Schuett
320 : ******************************************************************************/
321 1920301306 : static inline void ortho_cxyz_to_cxy(const int lp, const int k1, const int k2,
322 : const int cmax,
323 1920301306 : 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 7168126701 : for (int lzp = 0; lzp <= lp; lzp++) {
328 16362553870 : for (int lyp = 0; lyp <= lp - lzp; lyp++) {
329 31643631437 : for (int lxp = 0; lxp <= lp - lzp - lyp; lxp++) {
330 20528902962 : const double p1 = pol[2][lzp][k1 + cmax];
331 20528902962 : const double p2 = pol[2][lzp][k2 + cmax];
332 20528902962 : const int cxyz_index =
333 20528902962 : lzp * (lp + 1) * (lp + 1) + lyp * (lp + 1) + lxp; // [lzp, lyp, lxp]
334 20528902962 : const int cxy_index = lyp * (lp + 1) * 2 + lxp * 2; // [lyp, lxp, 0]
335 :
336 : #if (GRID_DO_COLLOCATE)
337 : // collocate
338 10351861274 : cxy[cxy_index + 0] += cxyz[cxyz_index] * p1;
339 10351861274 : cxy[cxy_index + 1] += cxyz[cxyz_index] * p2;
340 : #else
341 : // integrate
342 10177041688 : cxyz[cxyz_index] += cxy[cxy_index + 0] * p1;
343 10177041688 : cxyz[cxyz_index] += cxy[cxy_index + 1] * p2;
344 : #endif
345 : }
346 : }
347 : }
348 1920301306 : }
349 :
350 : /*******************************************************************************
351 : * \brief Collocates coefficients C_xyz onto the grid for orthorhombic case.
352 : * \author Ole Schuett
353 : ******************************************************************************/
354 : static inline void
355 180542706 : 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 180542706 : 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 180542706 : int cubecenter[3];
372 722170824 : for (int i = 0; i < 3; i++) {
373 : double dh_inv_rp = 0.0;
374 2166512472 : for (int j = 0; j < 3; j++) {
375 1624884354 : dh_inv_rp += dh_inv[j][i] * rp[j];
376 : }
377 541628118 : cubecenter[i] = (int)floor(dh_inv_rp);
378 : }
379 :
380 : double roffset[3];
381 722170824 : for (int i = 0; i < 3; i++) {
382 541628118 : roffset[i] = rp[i] - ((double)cubecenter[i]) * dh[i][i];
383 : }
384 :
385 : // Lookup loop bounds for spherical cutoff.
386 180542706 : int *sphere_bounds;
387 180542706 : double disr_radius;
388 180542706 : grid_sphere_cache_lookup(radius, dh, dh_inv, &sphere_bounds, &disr_radius);
389 180542706 : int **sphere_bounds_iter = &sphere_bounds;
390 :
391 : // Cube bounds.
392 180542706 : int lb_cube[3], ub_cube[3];
393 722170824 : for (int i = 0; i < 3; i++) {
394 541628118 : lb_cube[i] = (int)ceil(-1e-8 - disr_radius * dh_inv[i][i]);
395 541628118 : ub_cube[i] = 1 - lb_cube[i];
396 : // If grid is not period check that cube fits without wrapping.
397 541628118 : 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 180542706 : 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 180542706 : double pol_mutable[3][lp + 1][2 * cmax + 1];
411 722170824 : for (int idir = 0; idir < 3; idir++) {
412 541628118 : const double dr = dh[idir][idir];
413 541628118 : 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 541628118 : const double t_exp_1 = exp(-zetp * pow(dr, 2));
418 541628118 : const double t_exp_2 = pow(t_exp_1, 2);
419 541628118 : double t_exp_min_1 = exp(-zetp * pow(+dr - ro, 2));
420 541628118 : double t_exp_min_2 = exp(-2 * zetp * (+dr - ro) * (-dr));
421 6312273528 : for (int ig = 0; ig >= lb_cube[idir]; ig--) {
422 5770645410 : const double rpg = ig * dr - ro;
423 5770645410 : t_exp_min_1 *= t_exp_min_2 * t_exp_1;
424 5770645410 : t_exp_min_2 *= t_exp_2;
425 5770645410 : double pg = t_exp_min_1;
426 21542229481 : for (int icoef = 0; icoef <= lp; icoef++) {
427 15771584071 : pol_mutable[idir][icoef][ig + cmax] = pg; // exp(-zetp*rpg**2)
428 15771584071 : pg *= rpg;
429 : }
430 : }
431 541628118 : double t_exp_plus_1 = exp(-zetp * pow(-ro, 2));
432 541628118 : double t_exp_plus_2 = exp(-2 * zetp * (-ro) * (+dr));
433 6312273528 : for (int ig = 0; ig >= lb_cube[idir]; ig--) {
434 5770645410 : const double rpg = (1 - ig) * dr - ro;
435 5770645410 : t_exp_plus_1 *= t_exp_plus_2 * t_exp_1;
436 5770645410 : t_exp_plus_2 *= t_exp_2;
437 5770645410 : double pg = t_exp_plus_1;
438 21542229481 : for (int icoef = 0; icoef <= lp; icoef++) {
439 15771584071 : pol_mutable[idir][icoef][1 - ig + cmax] = pg; // exp(-zetp*rpg**2)
440 15771584071 : pg *= rpg;
441 : }
442 : }
443 : }
444 180542706 : const double (*pol)[lp + 1][2 * cmax + 1] =
445 180542706 : (const double (*)[lp + 1][2 * cmax + 1]) pol_mutable;
446 :
447 : // Precompute mapping from cube to grid indices for each direction
448 180542706 : int map_mutable[3][2 * cmax + 1];
449 722170824 : for (int i = 0; i < 3; i++) {
450 12709102386 : for (int k = -cmax; k <= +cmax; k++) {
451 12167474268 : map_mutable[i][k + cmax] =
452 12167474268 : modulo(cubecenter[i] + k - shift_local[i], npts_global[i]);
453 : }
454 : }
455 180542706 : 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 180542706 : int sections_mutable[3][2 * cmax + 1];
459 722170824 : for (int i = 0; i < 3; i++) {
460 12709102386 : for (int kg = 2 * cmax; kg >= 0; kg--) {
461 12167474268 : if (kg == 2 * cmax || map[i][kg] != map[i][kg + 1] - 1) {
462 1166029076 : sections_mutable[i][kg] = 0;
463 : } else {
464 11001445192 : sections_mutable[i][kg] = sections_mutable[i][kg + 1] + 1;
465 : }
466 : }
467 : }
468 180542706 : const int (*sections)[2 * cmax + 1] =
469 180542706 : (const int (*)[2 * cmax + 1]) sections_mutable;
470 :
471 : // Loop over k dimension of the cube.
472 180542706 : const int kstart = *((*sphere_bounds_iter)++);
473 180542706 : const size_t cxy_size = (lp + 1) * (lp + 1) * 2;
474 180542706 : double cxy[cxy_size];
475 2100844012 : for (int k1 = kstart; k1 <= 0; k1++) {
476 1920301306 : const int k2 = 1 - k1;
477 1920301306 : const int kg1 = map[2][k1 + cmax];
478 1920301306 : const int kg2 = map[2][k2 + cmax];
479 :
480 1920301306 : memset(cxy, 0, cxy_size * sizeof(double));
481 :
482 : #if (GRID_DO_COLLOCATE)
483 : // collocate
484 989587955 : ortho_cxyz_to_cxy(lp, k1, k2, cmax, pol, cxyz, cxy);
485 989587955 : ortho_cxy_to_grid(lp, kg1, kg2, cmax, pol, map, sections, npts_local,
486 : sphere_bounds_iter, cxy, grid);
487 : #else
488 : // integrate
489 930713351 : ortho_cxy_to_grid(lp, kg1, kg2, cmax, pol, map, sections, npts_local,
490 : sphere_bounds_iter, cxy, grid);
491 930713351 : ortho_cxyz_to_cxy(lp, k1, k2, cmax, pol, cxyz, cxy);
492 : #endif
493 : }
494 180542706 : }
495 :
496 : /*******************************************************************************
497 : * \brief Collocates coefficients C_i onto the grid for general case.
498 : * \author Ole Schuett
499 : ******************************************************************************/
500 5608353970 : 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 5608353970 : const int base = kg * npts_local[1] * npts_local[0] + jg * npts_local[0];
509 5608353970 : const int stride_i = index_max[0] - index_min[0] + 1;
510 5608353970 : const int stride_j = index_max[1] - index_min[1] + 1;
511 5608353970 : 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 17715156039 : for (int istart = ismin; istart <= ismax; istart++) {
532 12106802069 : const int istop = imin(ismax, istart + sections_i[istart - index_min[0]]);
533 12106802069 : if (map_i[istart - index_min[0]] < 0) {
534 0 : istart = istop; // skip over out-of-bounds indicies
535 0 : continue;
536 : }
537 :
538 12106802069 : 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 12106802069 : const int istart_tail = istart;
575 : #endif
576 :
577 >11104*10^7 : for (int i = istart_tail; i <= istop; i++) {
578 98940485074 : const int ig = i + cube2grid;
579 98940485074 : const double di = i - gp[0];
580 :
581 98940485074 : const int idx_ij = (j - index_min[1]) * stride_i + i - index_min[0];
582 98940485074 : 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 98940485074 : const double gaussian = exp_ij[idx_ij] * exp_jk[idx_jk] * exp_ik[idx_ik];
596 :
597 98940485074 : const int grid_index = base + ig; // [kg, jg, ig]
598 98940485074 : double dip = gaussian;
599 :
600 : #if (GRID_DO_COLLOCATE)
601 : // collocate
602 54709449705 : double reg = 0.0;
603 >22999*10^7 : for (int il = 0; il <= lp; il++) {
604 >17528*10^7 : reg += ci[il] * dip;
605 >17528*10^7 : dip *= di;
606 : }
607 54709449705 : grid[grid_index] += reg;
608 : #else
609 : // integrate
610 44231035369 : const double reg = grid[grid_index];
611 >18547*10^7 : for (int il = 0; il <= lp; il++) {
612 >14124*10^7 : ci[il] += reg * dip;
613 >14124*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 3110597089 : 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 3110597089 : double djp = 1.0;
641 31338066537 : for (int jl = 0; jl <= lp; jl++) {
642 59752741016 : for (int il = 0; il <= lp - jl; il++) {
643 41943491694 : const int cij_index = jl * (lp + 1) + il; // [jl, il]
644 : #if (GRID_DO_COLLOCATE)
645 23078338421 : ci[il] += cij[cij_index] * djp; // collocate
646 : #else
647 18865153273 : cij[cij_index] += ci[il] * djp; // integrate
648 : #endif
649 : }
650 17809249322 : 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 5608353970 : 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 6221194178 : general_cij_to_ci(lp, dj, cij, ci);
670 3110597089 : 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 10418220126 : 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 302173874 : 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 302173874 : const double a =
696 302173874 : dh[0][0] * dh[0][0] + dh[0][1] * dh[0][1] + dh[0][2] * dh[0][2];
697 302173874 : const double inv_2a = 1.0 / (2.0 * a);
698 302173874 : const double radius2 = radius * radius;
699 302173874 : const double gp0_dh0[3] = {(0 - gp[0]) * dh[0][0], (0 - gp[0]) * dh[0][1],
700 302173874 : (0 - gp[0]) * dh[0][2]};
701 302173874 : const double dk_dh2[3] = {(k - gp[2]) * dh[2][0], (k - gp[2]) * dh[2][1],
702 302173874 : (k - gp[2]) * dh[2][2]};
703 :
704 8776238670 : for (int j = index_min[1]; j <= index_max[1]; j++) {
705 8474064796 : const int jg = map_j[j - index_min[1]];
706 8474064796 : 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 8474064796 : const double dj = j - gp[1];
731 8474064796 : double b = 0.0, c = 0.0;
732 33896259184 : for (int i = 0; i < 3; i++) {
733 25422194388 : const double v = (gp0_dh0[i] + dj * dh[1][i]) + dk_dh2[i];
734 25422194388 : b += 2.0 * v * dh[0][i];
735 25422194388 : c += v * v;
736 : }
737 8474064796 : const double d = b * b - 4.0 * a * (c - radius2);
738 :
739 8474064796 : if (0.0 < d) {
740 5608353970 : const double sqrt_d = sqrt(d);
741 5608353970 : const int ismin = (int)ceil((-b - sqrt_d) * inv_2a);
742 5608353970 : const int ismax = (int)floor((-b + sqrt_d) * inv_2a);
743 :
744 5608353970 : double ci[lp + 1];
745 5608353970 : memset(ci, 0, sizeof(ci));
746 :
747 : // Separate branches for low values of lp are generated in
748 : // general_cijk_to_grid.
749 11216707940 : 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 166860491 : }
755 :
756 : /*******************************************************************************
757 : * \brief Transforms coefficients C_ijk into C_ij by fixing grid index k.
758 : * \author Ole Schuett
759 : ******************************************************************************/
760 302173874 : 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 302173874 : double dkp = 1.0;
764 1246042155 : for (int kl = 0; kl <= lp; kl++) {
765 3149919108 : for (int jl = 0; jl <= lp - kl; jl++) {
766 6627405181 : for (int il = 0; il <= lp - kl - jl; il++) {
767 4421354354 : const int cij_index = jl * (lp + 1) + il; // [jl, il]
768 4421354354 : const int cijk_index =
769 4421354354 : kl * (lp + 1) * (lp + 1) + jl * (lp + 1) + il; // [kl, jl, il]
770 : #if (GRID_DO_COLLOCATE)
771 2375344251 : cij[cij_index] += cijk[cijk_index] * dkp; // collocate
772 : #else
773 2046010103 : cijk[cijk_index] += cij[cij_index] * dkp; // integrate
774 : #endif
775 : }
776 : }
777 943868281 : dkp *= dk;
778 : }
779 302173874 : }
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 38010828 : 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 38010828 : bool has_valid = false;
793 984121358 : for (int k = index_min; k <= index_max; k++) {
794 946110530 : const int kg = modulo(k - shift_local, npts_global);
795 946110530 : if (bounds[0] <= kg && kg <= bounds[1]) {
796 946105323 : map[k - index_min] = kg;
797 946105323 : has_valid = true;
798 : } else {
799 5207 : map[k - index_min] = INT_MIN; // out of bounds - not mapped
800 : }
801 : }
802 :
803 38010828 : if (!has_valid) {
804 : return false;
805 : }
806 :
807 : // Precompute length of sections with homogeneous cube to grid mapping.
808 38010778 : const int range = index_max - index_min + 1;
809 984120434 : for (int kg = range - 1; kg >= 0; kg--) {
810 946109656 : if (kg == range - 1 || map[kg] != map[kg + 1] - 1) {
811 95509121 : sections[kg] = 0;
812 : } else {
813 850600535 : 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 38010681 : 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 38010681 : double exp_table[]) {
828 :
829 38010681 : const int stride_i = index_max[idir] - index_min[idir] + 1;
830 38010681 : const double h_ii = dh[idir][0] * dh[idir][0] + dh[idir][1] * dh[idir][1] +
831 38010681 : dh[idir][2] * dh[idir][2];
832 38010681 : const double h_jj = dh[jdir][0] * dh[jdir][0] + dh[jdir][1] * dh[jdir][1] +
833 38010681 : dh[jdir][2] * dh[jdir][2];
834 38010681 : const double h_ij = dh[idir][0] * dh[jdir][0] + dh[idir][1] * dh[jdir][1] +
835 38010681 : dh[idir][2] * dh[jdir][2];
836 :
837 38010681 : double exp_half_i[stride_i];
838 1005059299 : for (int i = index_min[idir]; i <= index_max[idir]; i++) {
839 967048618 : const double di = i - gp[idir];
840 967048618 : exp_half_i[i - index_min[idir]] = exp(-0.5 * zetp * di * di * h_ii);
841 : }
842 :
843 963177921 : for (int j = index_min[jdir]; j <= index_max[jdir]; j++) {
844 925167240 : const double dj = j - gp[jdir];
845 925167240 : const double rjj = 0.5 * dj * dj * h_jj;
846 925167240 : const double rij_unit = dj * h_ij;
847 925167240 : const double exp_ij_unit = exp(-zetp * 2.0 * rij_unit);
848 :
849 : // compute exponentials symmetrically around cube center
850 925167240 : const int i_center = (int)gp[idir];
851 925167240 : const double di_center = i_center - gp[idir];
852 925167240 : const double rij_center = di_center * rij_unit;
853 925167240 : const double exp_ij_center = exp(-zetp * (rjj + 2.0 * rij_center));
854 :
855 : // above center
856 925167240 : double exp_ij = exp_ij_center;
857 14463011537 : for (int i = i_center; i <= index_max[idir]; i++) {
858 13537844297 : const int idx = (j - index_min[jdir]) * stride_i + i - index_min[idir];
859 13537844297 : exp_table[idx] = exp_ij * exp_half_i[i - index_min[idir]];
860 13537844297 : exp_ij *= exp_ij_unit;
861 : }
862 :
863 : // below center
864 925167240 : const double exp_ij_unit_inv = 1.0 / exp_ij_unit;
865 925167240 : exp_ij = exp_ij_center * exp_ij_unit_inv;
866 13468820529 : for (int i = i_center - 1; i >= index_min[idir]; i--) {
867 12543653289 : const int idx = (j - index_min[jdir]) * stride_i + i - index_min[idir];
868 12543653289 : exp_table[idx] = exp_ij * exp_half_i[i - index_min[idir]];
869 12543653289 : exp_ij *= exp_ij_unit_inv;
870 : }
871 : }
872 38010681 : }
873 :
874 : /*******************************************************************************
875 : * \brief K-loop body of general_cijk_to_grid to be inlined for low values of
876 : * lp.
877 : ******************************************************************************/
878 12670227 : 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 12670227 : const int cij_size = (lp + 1) * (lp + 1);
888 314848434 : for (int k = index_min[2]; k <= index_max[2]; k++) {
889 302178207 : const int kg = map_k[k - index_min[2]];
890 302178207 : 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 302173874 : memset(cij, 0, cij_size * sizeof(double));
897 :
898 : #if (GRID_DO_COLLOCATE)
899 : // collocate
900 166860491 : general_cijk_to_cij(lp, (double)k - gp[2], cijk, cij);
901 333720982 : 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 135313383 : 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 135313383 : 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 12670277 : 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 12670277 : GRID_CONST_WHEN_INTEGRATE double *grid) {
929 :
930 : // Default for border_mask == 0.
931 12670277 : int bounds_i[2] = {0, npts_local[0] - 1};
932 12670277 : int bounds_j[2] = {0, npts_local[1] - 1};
933 12670277 : 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 12670277 : if (border_mask & (1 << 0))
938 150201 : bounds_i[0] += border_width[0];
939 12670277 : if (border_mask & (1 << 1))
940 150201 : bounds_i[1] -= border_width[0];
941 12670277 : if (border_mask & (1 << 2))
942 150201 : bounds_j[0] += border_width[1];
943 12670277 : if (border_mask & (1 << 3))
944 150201 : bounds_j[1] -= border_width[1];
945 12670277 : if (border_mask & (1 << 4))
946 150193 : bounds_k[0] += border_width[2];
947 12670277 : if (border_mask & (1 << 5))
948 150201 : bounds_k[1] -= border_width[2];
949 :
950 12670277 : if (bounds_i[0] > bounds_i[1] || bounds_j[0] > bounds_j[1] ||
951 12670277 : bounds_k[0] > bounds_k[1]) {
952 50 : return;
953 : }
954 :
955 : // center in grid coords
956 : // gp = MATMUL(dh_inv, rp)
957 12670277 : double gp[3] = {0.0, 0.0, 0.0};
958 50681108 : for (int i = 0; i < 3; i++) {
959 152043324 : for (int j = 0; j < 3; j++) {
960 114032493 : 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 50681108 : for (int idir = 0; idir < 3; idir++) {
974 38010831 : const double norm = sqrt(dh_inv[0][idir] * dh_inv[0][idir] +
975 38010831 : dh_inv[1][idir] * dh_inv[1][idir] +
976 38010831 : dh_inv[2][idir] * dh_inv[2][idir]);
977 38010831 : const double rad_idx = radius * norm * sphere_bound_eps;
978 38010831 : index_min[idir] = (int)floor(gp[idir] - rad_idx);
979 38010831 : index_max[idir] = (int)ceil(gp[idir] + rad_idx);
980 : }
981 :
982 : // Precompute mappings
983 12670277 : const int range_i = index_max[0] - index_min[0] + 1;
984 12670277 : int map_i[range_i], sections_i[range_i];
985 12670277 : if (!general_precompute_mapping(index_min[0], index_max[0], shift_local[0],
986 : npts_global[0], bounds_i, map_i,
987 : sections_i)) {
988 : return;
989 : }
990 12670277 : const int range_j = index_max[1] - index_min[1] + 1;
991 12670277 : int map_j[range_j], sections_j[range_j];
992 12670277 : if (!general_precompute_mapping(index_min[1], index_max[1], shift_local[1],
993 : npts_global[1], bounds_j, map_j,
994 : sections_j)) {
995 : return;
996 : }
997 12670274 : const int range_k = index_max[2] - index_min[2] + 1;
998 12670274 : int map_k[range_k], sections_k[range_k];
999 12670274 : if (!general_precompute_mapping(index_min[2], index_max[2], shift_local[2],
1000 : npts_global[2], bounds_k, map_k,
1001 : sections_k)) {
1002 : return;
1003 : }
1004 :
1005 : // Precompute exponentials
1006 12670227 : double exp_ij[range_i * range_j];
1007 12670227 : general_fill_exp_table(0, 1, index_min, index_max, zetp, dh, gp, exp_ij);
1008 12670227 : double exp_jk[range_j * range_k];
1009 12670227 : general_fill_exp_table(1, 2, index_min, index_max, zetp, dh, gp, exp_jk);
1010 12670227 : double exp_ik[range_k * range_i];
1011 12670227 : general_fill_exp_table(0, 2, index_min, index_max, zetp, dh, gp, exp_ik);
1012 :
1013 : // go over the grid, but cycle if the point is not within the radius
1014 12670227 : const int cij_size = (lp + 1) * (lp + 1);
1015 12670227 : double cij[cij_size];
1016 : // Generate separate branches for low values of lp.
1017 12670227 : if (lp <= GRID_MAX_LP_OPTIMIZED) {
1018 : GRID_PRAGMA_UNROLL(GRID_MAX_LP_OPTIMIZED + 1)
1019 139372497 : for (int ilp = 0; ilp <= GRID_MAX_LP_OPTIMIZED; ilp++) {
1020 126702270 : if (lp == ilp) {
1021 139372497 : general_cijk_kloop_low(ilp, npts_local, index_min, index_max, map_i,
1022 : map_j, map_k, sections_i, sections_j, sections_k,
1023 : dh, gp, radius, exp_ij, exp_jk, exp_ik, cij,
1024 : cijk, grid);
1025 : }
1026 : }
1027 : } else {
1028 12670227 : general_cijk_kloop_low(lp, npts_local, index_min, index_max, map_i, map_j,
1029 : map_k, sections_i, sections_j, sections_k, dh, gp,
1030 : radius, exp_ij, exp_jk, exp_ik, cij, cijk, grid);
1031 : }
1032 : }
1033 :
1034 : /*******************************************************************************
1035 : * \brief Transforms coefficients C_xyz into C_ijk.
1036 : * \author Ole Schuett
1037 : ******************************************************************************/
1038 : static inline void
1039 12670277 : general_cxyz_to_cijk(const int lp, const double dh[3][3],
1040 : GRID_CONST_WHEN_COLLOCATE double *cxyz,
1041 12670277 : GRID_CONST_WHEN_INTEGRATE double *cijk) {
1042 :
1043 : // transform P_{lxp,lyp,lzp} into a P_{lip,ljp,lkp} such that
1044 : // sum_{lxp,lyp,lzp} P_{lxp,lyp,lzp} (x-x_p)**lxp (y-y_p)**lyp (z-z_p)**lzp =
1045 : // sum_{lip,ljp,lkp} P_{lip,ljp,lkp} (i-i_p)**lip (j-j_p)**ljp (k-k_p)**lkp
1046 :
1047 : // transform using multinomials
1048 12670277 : const int lp1 = lp + 1;
1049 12670277 : double hmatgridp[lp1][3][3];
1050 50681108 : for (int i = 0; i < 3; i++) {
1051 152043324 : for (int j = 0; j < 3; j++) {
1052 114032493 : hmatgridp[0][j][i] = 1.0;
1053 350629083 : for (int k = 1; k <= lp; k++) {
1054 236596590 : hmatgridp[k][j][i] = hmatgridp[k - 1][j][i] * dh[j][i];
1055 : }
1056 : }
1057 : }
1058 :
1059 12670277 : double trinom[3][lp1][lp1][lp1];
1060 50681108 : for (int d = 0; d < 3; d++) {
1061 154887192 : for (int kl = 0; kl <= lp; kl++) {
1062 388999854 : for (int jl = 0; jl <= lp - kl; jl++) {
1063 817348884 : for (int il = 0; il <= lp - kl - jl; il++) {
1064 545225391 : const int l = il + jl + kl;
1065 545225391 : trinom[d][kl][jl][il] = hmatgridp[il][0][d] * hmatgridp[jl][1][d] *
1066 545225391 : hmatgridp[kl][2][d] * fac(l) /
1067 545225391 : (fac(il) * fac(jl) * fac(kl));
1068 : }
1069 : }
1070 : }
1071 : }
1072 :
1073 : const int lpx = lp;
1074 51629064 : for (int klx = 0; klx <= lpx; klx++) {
1075 129666618 : for (int jlx = 0; jlx <= lpx - klx; jlx++) {
1076 272449628 : for (int ilx = 0; ilx <= lpx - klx - jlx; ilx++) {
1077 181741797 : const int lx = ilx + jlx + klx;
1078 181741797 : const double px = trinom[0][klx][jlx][ilx];
1079 181741797 : const int lpy = lp - lx;
1080 513144317 : for (int kly = 0; kly <= lpy; kly++) {
1081 897579593 : for (int jly = 0; jly <= lpy - kly; jly++) {
1082 1487909928 : for (int ily = 0; ily <= lpy - kly - jly; ily++) {
1083 921732855 : const int ly = ily + jly + kly;
1084 921732855 : const double pxy = px * trinom[1][kly][jly][ily];
1085 921732855 : const int il_xy = ilx + ily;
1086 921732855 : const int jl_xy = jlx + jly;
1087 921732855 : const int kl_xy = klx + kly;
1088 921732855 : const int lpz = lp - lx - ly;
1089 921732855 : const int cxyz_base = ly * lp1 + lx;
1090 2367142826 : for (int klz = 0; klz <= lpz; klz++) {
1091 1445409971 : const int kl_base = (kl_xy + klz) * lp1 * lp1;
1092 3644632007 : for (int jlz = 0; jlz <= lpz - klz; jlz++) {
1093 2199222036 : const int cijk_base = kl_base + (jl_xy + jlz) * lp1 + il_xy;
1094 5462639773 : for (int ilz = 0; ilz <= lpz - klz - jlz; ilz++) {
1095 3263417737 : const int lz = ilz + jlz + klz;
1096 3263417737 : const int cijk_index = cijk_base + ilz; // [kl,jl,il]
1097 3263417737 : const int cxyz_index =
1098 3263417737 : lz * lp1 * lp1 + cxyz_base; // [lz,ly,lx]
1099 3263417737 : const double p = pxy * trinom[2][klz][jlz][ilz];
1100 : #if (GRID_DO_COLLOCATE)
1101 1430625932 : cijk[cijk_index] += cxyz[cxyz_index] * p; // collocate
1102 : #else
1103 1832791805 : cxyz[cxyz_index] += cijk[cijk_index] * p; // integrate
1104 : #endif
1105 : }
1106 : }
1107 : }
1108 : }
1109 : }
1110 : }
1111 : }
1112 : }
1113 : }
1114 12670277 : }
1115 :
1116 : /*******************************************************************************
1117 : * \brief Collocates coefficients C_xyz onto the grid for general case.
1118 : * \author Ole Schuett
1119 : ******************************************************************************/
1120 : static inline void
1121 12670277 : general_cxyz_to_grid(const int border_mask, const int lp, const double zetp,
1122 : const double dh[3][3], const double dh_inv[3][3],
1123 : const double rp[3], const int npts_global[3],
1124 : const int npts_local[3], const int shift_local[3],
1125 : const int border_width[3], const double radius,
1126 : GRID_CONST_WHEN_COLLOCATE double *cxyz,
1127 12670277 : GRID_CONST_WHEN_INTEGRATE double *grid) {
1128 :
1129 12670277 : const size_t cijk_size = (lp + 1) * (lp + 1) * (lp + 1);
1130 12670277 : double cijk[cijk_size];
1131 12670277 : memset(cijk, 0, cijk_size * sizeof(double));
1132 :
1133 : #if (GRID_DO_COLLOCATE)
1134 : // collocate
1135 6936527 : general_cxyz_to_cijk(lp, dh, cxyz, cijk);
1136 6936527 : general_cijk_to_grid(border_mask, lp, zetp, dh, dh_inv, rp, npts_global,
1137 : npts_local, shift_local, border_width, radius, cijk,
1138 : grid);
1139 : #else
1140 : // integrate
1141 5733750 : general_cijk_to_grid(border_mask, lp, zetp, dh, dh_inv, rp, npts_global,
1142 : npts_local, shift_local, border_width, radius, cijk,
1143 : grid);
1144 5733750 : general_cxyz_to_cijk(lp, dh, cxyz, cijk);
1145 : #endif
1146 12670277 : }
1147 :
1148 : /*******************************************************************************
1149 : * \brief Collocates coefficients C_xyz onto the grid.
1150 : * \author Ole Schuett
1151 : ******************************************************************************/
1152 : static inline void
1153 193212983 : cxyz_to_grid(const bool orthorhombic, const int border_mask, const int lp,
1154 : const double zetp, const double dh[3][3],
1155 : const double dh_inv[3][3], const double rp[3],
1156 : const int npts_global[3], const int npts_local[3],
1157 : const int shift_local[3], const int border_width[3],
1158 : const double radius, GRID_CONST_WHEN_COLLOCATE double *cxyz,
1159 : GRID_CONST_WHEN_INTEGRATE double *grid) {
1160 :
1161 193212983 : enum grid_library_kernel k;
1162 193212983 : if (orthorhombic && border_mask == 0) {
1163 180542706 : k = (GRID_DO_COLLOCATE) ? GRID_COLLOCATE_ORTHO : GRID_INTEGRATE_ORTHO;
1164 180542706 : ortho_cxyz_to_grid(lp, zetp, dh, dh_inv, rp, npts_global, npts_local,
1165 : shift_local, radius, cxyz, grid);
1166 : } else {
1167 12670277 : k = (GRID_DO_COLLOCATE) ? GRID_COLLOCATE_GENERAL : GRID_INTEGRATE_GENERAL;
1168 12670277 : general_cxyz_to_grid(border_mask, lp, zetp, dh, dh_inv, rp, npts_global,
1169 : npts_local, shift_local, border_width, radius, cxyz,
1170 : grid);
1171 : }
1172 193212983 : grid_library_counter_add(lp, GRID_BACKEND_CPU, k, 1);
1173 193212983 : }
1174 :
1175 : /*******************************************************************************
1176 : * \brief Transforms coefficients C_ab into C_xyz.
1177 : * \author Ole Schuett
1178 : ******************************************************************************/
1179 193212983 : static inline void cab_to_cxyz(const int la_max, const int la_min,
1180 : const int lb_max, const int lb_min,
1181 : const double prefactor, const double ra[3],
1182 : const double rb[3], const double rp[3],
1183 : GRID_CONST_WHEN_COLLOCATE double *cab,
1184 193212983 : GRID_CONST_WHEN_INTEGRATE double *cxyz) {
1185 :
1186 : // Computes the polynomial expansion coefficients:
1187 : // (x-a)**lxa (x-b)**lxb -> sum_{ls} alpha(ls,lxa,lxb,1)*(x-p)**ls
1188 193212983 : const int lp = la_max + lb_max;
1189 193212983 : const int lp1 = lp + 1;
1190 193212983 : double alpha[3][lb_max + 1][la_max + 1][lp + 1];
1191 :
1192 772851932 : for (int i = 0; i < 3; i++) {
1193 579638949 : const double drpa = rp[i] - ra[i];
1194 579638949 : const double drpb = rp[i] - rb[i];
1195 579638949 : alpha[i][0][0][0] = 1.0;
1196 1089671556 : for (int lxa = 1; lxa <= la_max; lxa++) {
1197 510032607 : alpha[i][0][lxa][0] = drpa * alpha[i][0][lxa - 1][0];
1198 622556484 : for (int k = 1; k < lxa; k++) {
1199 112523877 : alpha[i][0][lxa][k] =
1200 112523877 : alpha[i][0][lxa - 1][k - 1] + drpa * alpha[i][0][lxa - 1][k];
1201 : }
1202 510032607 : alpha[i][0][lxa][lxa] = 1.0;
1203 : }
1204 1054191153 : for (int lxb = 1; lxb <= lb_max; lxb++) {
1205 1505574078 : for (int lxa = 0; lxa <= la_max; lxa++) {
1206 1031021874 : const int n = lxa + lxb;
1207 1031021874 : alpha[i][lxb][lxa][0] = drpb * alpha[i][lxb - 1][lxa][0];
1208 1964140485 : for (int k = 1; k < n; k++) {
1209 933118611 : alpha[i][lxb][lxa][k] =
1210 933118611 : alpha[i][lxb - 1][lxa][k - 1] + drpb * alpha[i][lxb - 1][lxa][k];
1211 : }
1212 1031021874 : alpha[i][lxb][lxa][n] = 1.0;
1213 : }
1214 : }
1215 : }
1216 :
1217 : // *** initialise the coefficient matrix, we transform the sum
1218 : //
1219 : // sum_{lxa,lya,lza,lxb,lyb,lzb} P_{lxa,lya,lza,lxb,lyb,lzb} *
1220 : // (x-a_x)**lxa (y-a_y)**lya (z-a_z)**lza (x-b_x)**lxb (y-a_y)**lya
1221 : // (z-a_z)**lza
1222 : //
1223 : // into
1224 : //
1225 : // sum_{lxp,lyp,lzp} P_{lxp,lyp,lzp} (x-p_x)**lxp (y-p_y)**lyp (z-p_z)**lzp
1226 : //
1227 : // where p is center of the product gaussian, and lp = la_max + lb_max
1228 : // (current implementation is l**7)
1229 : //
1230 : // Factor the 3D transformation (lxa,lya,lza,lxb,lyb,lzb) <-> (lxp,lyp,lzp)
1231 : // into three separable 1D steps along x, y, z.
1232 193212983 : const int ncoset_a = ncoset(la_max);
1233 193212983 : double txy[lp1 * lp1];
1234 193212983 : double tx[lp1];
1235 :
1236 : #if (GRID_DO_COLLOCATE)
1237 277958105 : for (int lzb = 0; lzb <= lb_max; lzb++) {
1238 533848324 : for (int lza = 0; lza <= la_max; lza++) {
1239 353913212 : const int lz_max = lza + lzb;
1240 353913212 : const int lyp_max = lp - lz_max;
1241 353913212 : memset(txy, 0, (lyp_max + 1) * lp1 * sizeof(double));
1242 921196624 : for (int lyb = 0; lyb <= lb_max - lzb; lyb++) {
1243 1495666343 : for (int lya = 0; lya <= la_max - lza; lya++) {
1244 928382931 : const int ly_max = lya + lyb;
1245 928382931 : const int lxb_min = imax(lb_min - lzb - lyb, 0);
1246 928382931 : const int lxa_min = imax(la_min - lza - lya, 0);
1247 928382931 : const int lxb_max = lb_max - lzb - lyb;
1248 928382931 : const int lxa_max = la_max - lza - lya;
1249 928382931 : const int lxp_max = lxa_max + lxb_max;
1250 928382931 : memset(tx, 0, (lxp_max + 1) * sizeof(double));
1251 2211744737 : for (int lxb = lxb_min; lxb <= lxb_max; lxb++) {
1252 1283361806 : const int jco = coset(lxb, lyb, lzb);
1253 3099929439 : for (int lxa = lxa_min; lxa <= lxa_max; lxa++) {
1254 1816567633 : const int ico = coset(lxa, lya, lza);
1255 1816567633 : const int cab_index = jco * ncoset_a + ico; // [jco, ico]
1256 1816567633 : const double cab_val = cab[cab_index] * prefactor;
1257 5133336152 : for (int lxp = 0; lxp <= lxa + lxb; lxp++) {
1258 3316768519 : tx[lxp] += cab_val * alpha[0][lxb][lxa][lxp];
1259 : }
1260 : }
1261 : }
1262 2735369382 : for (int lyp = 0; lyp <= ly_max; lyp++) {
1263 1806986451 : const double ay = alpha[1][lyb][lya][lyp];
1264 1806986451 : double *txy_row = &txy[lyp * lp1];
1265 5129706055 : for (int lxp = 0; lxp <= lxp_max; lxp++) {
1266 3322719604 : txy_row[lxp] += tx[lxp] * ay;
1267 : }
1268 : }
1269 : }
1270 : }
1271 1132921296 : for (int lzp = 0; lzp <= lz_max; lzp++) {
1272 779008084 : const double az = alpha[2][lzb][lza][lzp];
1273 779008084 : double *cxyz_z = &cxyz[lzp * lp1 * lp1];
1274 2363535985 : for (int lyp = 0; lyp <= lyp_max; lyp++) {
1275 1584527901 : const double *txy_row = &txy[lyp * lp1];
1276 1584527901 : double *cxyz_zy = &cxyz_z[lyp * lp1];
1277 4356945449 : for (int lxp = 0; lxp <= lyp_max - lyp; lxp++) {
1278 2772417548 : cxyz_zy[lxp] += txy_row[lxp] * az;
1279 : }
1280 : }
1281 : }
1282 : }
1283 : }
1284 : #else
1285 266651929 : for (int lzb = 0; lzb <= lb_max; lzb++) {
1286 524446537 : for (int lza = 0; lza <= la_max; lza++) {
1287 352984598 : const int lz_max = lza + lzb;
1288 352984598 : const int lyp_max = lp - lz_max;
1289 352984598 : memset(txy, 0, (lyp_max + 1) * lp1 * sizeof(double));
1290 1143106647 : for (int lzp = 0; lzp <= lz_max; lzp++) {
1291 790122049 : const double az = alpha[2][lzb][lza][lzp] * prefactor;
1292 790122049 : const double *cxyz_z = &cxyz[lzp * lp1 * lp1];
1293 2424659647 : for (int lyp = 0; lyp <= lyp_max; lyp++) {
1294 1634537598 : double *txy_row = &txy[lyp * lp1];
1295 1634537598 : const double *cxyz_zy = &cxyz_z[lyp * lp1];
1296 4544436126 : for (int lxp = 0; lxp <= lyp_max - lyp; lxp++) {
1297 2909898528 : txy_row[lxp] += cxyz_zy[lxp] * az;
1298 : }
1299 : }
1300 : }
1301 909164182 : for (int lyb = 0; lyb <= lb_max - lzb; lyb++) {
1302 1502119217 : for (int lya = 0; lya <= la_max - lza; lya++) {
1303 945939633 : const int ly_max = lya + lyb;
1304 945939633 : const int lxb_min = imax(lb_min - lzb - lyb, 0);
1305 945939633 : const int lxa_min = imax(la_min - lza - lya, 0);
1306 945939633 : const int lxb_max = lb_max - lzb - lyb;
1307 945939633 : const int lxa_max = la_max - lza - lya;
1308 945939633 : const int lxp_max = lxa_max + lxb_max;
1309 945939633 : memset(tx, 0, (lxp_max + 1) * sizeof(double));
1310 2821078979 : for (int lyp = 0; lyp <= ly_max; lyp++) {
1311 1875139346 : const double ay = alpha[1][lyb][lya][lyp];
1312 1875139346 : const double *txy_row = &txy[lyp * lp1];
1313 5395012489 : for (int lxp = 0; lxp <= lxp_max; lxp++) {
1314 3519873143 : tx[lxp] += txy_row[lxp] * ay;
1315 : }
1316 : }
1317 2249781799 : for (int lxb = lxb_min; lxb <= lxb_max; lxb++) {
1318 1303842166 : const int jco = coset(lxb, lyb, lzb);
1319 3249077178 : for (int lxa = lxa_min; lxa <= lxa_max; lxa++) {
1320 1945235012 : const int ico = coset(lxa, lya, lza);
1321 1945235012 : const int cab_index = jco * ncoset_a + ico; // [jco, ico]
1322 1945235012 : double sum_x = 0.0;
1323 5568973968 : for (int lxp = 0; lxp <= lxa + lxb; lxp++) {
1324 3623738956 : sum_x += tx[lxp] * alpha[0][lxb][lxa][lxp];
1325 : }
1326 1945235012 : cab[cab_index] += sum_x;
1327 : }
1328 : }
1329 : }
1330 : }
1331 : }
1332 : }
1333 : #endif
1334 193212983 : }
1335 :
1336 : /*******************************************************************************
1337 : * \brief Collocates coefficients C_ab onto the grid.
1338 : * \author Ole Schuett
1339 : ******************************************************************************/
1340 : static inline void
1341 243219233 : cab_to_grid(const bool orthorhombic, const int border_mask, const int la_max,
1342 : const int la_min, const int lb_max, const int lb_min,
1343 : const double zeta, const double zetb, const double rscale,
1344 : const double dh[3][3], const double dh_inv[3][3],
1345 : const double ra[3], const double rab[3], const int npts_global[3],
1346 : const int npts_local[3], const int shift_local[3],
1347 : const int border_width[3], const double radius,
1348 : GRID_CONST_WHEN_COLLOCATE double *cab,
1349 243219233 : GRID_CONST_WHEN_INTEGRATE double *grid) {
1350 :
1351 : // Check if radius is too small to be mapped onto grid of given resolution.
1352 243219233 : double dh_max = 0.0;
1353 972876932 : for (int i = 0; i < 3; i++) {
1354 2918630796 : for (int j = 0; j < 3; j++) {
1355 2188973097 : dh_max = fmax(dh_max, fabs(dh[i][j]));
1356 : }
1357 : }
1358 243219233 : if (2.0 * radius < dh_max) {
1359 50006250 : return;
1360 : }
1361 :
1362 193212983 : const double zetp = zeta + zetb;
1363 193212983 : const double f = zetb / zetp;
1364 193212983 : const double rab2 = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
1365 193212983 : const double prefactor = rscale * exp(-zeta * f * rab2);
1366 193212983 : double rp[3], rb[3];
1367 772851932 : for (int i = 0; i < 3; i++) {
1368 579638949 : rp[i] = ra[i] + f * rab[i];
1369 579638949 : rb[i] = ra[i] + rab[i];
1370 : }
1371 :
1372 193212983 : const int lp = la_max + lb_max;
1373 193212983 : const size_t cxyz_size = (lp + 1) * (lp + 1) * (lp + 1);
1374 193212983 : double cxyz[cxyz_size];
1375 193212983 : memset(cxyz, 0, cxyz_size * sizeof(double));
1376 :
1377 : #if (GRID_DO_COLLOCATE)
1378 : // collocate
1379 98022993 : cab_to_cxyz(la_max, la_min, lb_max, lb_min, prefactor, ra, rb, rp, cab, cxyz);
1380 98022993 : cxyz_to_grid(orthorhombic, border_mask, lp, zetp, dh, dh_inv, rp, npts_global,
1381 : npts_local, shift_local, border_width, radius, cxyz, grid);
1382 : #else
1383 : // integrate
1384 95189990 : cxyz_to_grid(orthorhombic, border_mask, lp, zetp, dh, dh_inv, rp, npts_global,
1385 : npts_local, shift_local, border_width, radius, cxyz, grid);
1386 95189990 : cab_to_cxyz(la_max, la_min, lb_max, lb_min, prefactor, ra, rb, rp, cab, cxyz);
1387 : #endif
1388 : }
1389 :
1390 : // EOF
|
