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 <math.h>
10 : #include <stddef.h>
11 : #include <stdio.h>
12 : #include <stdlib.h>
13 : #include <string.h>
14 :
15 : #include "grid_task_list.h"
16 : #include "grid_task_list_internal.h"
17 :
18 : /*******************************************************************************
19 : * \brief Allocates a task list which can be passed to grid_collocate_task_list.
20 : * See grid_task_list.h for details.
21 : * \author Ole Schuett
22 : ******************************************************************************/
23 17168 : void grid_create_task_list(
24 : const bool orthorhombic, const int ntasks, const int nlevels,
25 : const int natoms, const int nkinds, const int nblocks,
26 : const int block_offsets[nblocks], const double atom_positions[natoms][3],
27 : const int atom_kinds[natoms], const grid_basis_set *basis_sets[nkinds],
28 : const int level_list[ntasks], const int iatom_list[ntasks],
29 : const int jatom_list[ntasks], const int iset_list[ntasks],
30 : const int jset_list[ntasks], const int ipgf_list[ntasks],
31 : const int jpgf_list[ntasks], const int border_mask_list[ntasks],
32 : const int block_num_list[ntasks], const double radius_list[ntasks],
33 : const double rab_list[ntasks][3], const int npts_global[nlevels][3],
34 : const int npts_local[nlevels][3], const int shift_local[nlevels][3],
35 : const int border_width[nlevels][3], const double dh[nlevels][3][3],
36 : const double dh_inv[nlevels][3][3], grid_task_list **task_list_out) {
37 :
38 17168 : const grid_library_config config = grid_library_get_config();
39 :
40 17168 : grid_task_list_internal *task_list = NULL;
41 :
42 17168 : if (*task_list_out == NULL) {
43 10334 : task_list = malloc(sizeof(grid_task_list_internal));
44 : // not a proper handling of errors. assert is a debug tool
45 10334 : assert(task_list != NULL);
46 10334 : memset(task_list, 0, sizeof(grid_task_list_internal));
47 10334 : *task_list_out = task_list;
48 :
49 : // Resolve AUTO to a concrete backend.
50 10334 : if (config.backend == GRID_BACKEND_AUTO) {
51 : #if (defined(__OFFLOAD_CUDA) || defined(__OFFLOAD_HIP)) && \
52 : !defined(__NO_OFFLOAD_GRID)
53 : task_list->backend = GRID_BACKEND_GPU;
54 : #else
55 10314 : task_list->backend = GRID_BACKEND_CPU;
56 : #endif
57 : } else {
58 20 : task_list->backend = config.backend;
59 : }
60 : } else {
61 : // Reuse existing task list.
62 : task_list = (grid_task_list_internal *)*task_list_out;
63 : }
64 :
65 17168 : if ((nblocks == 0) || (ntasks == 0) || (nlevels == 0)) {
66 322 : task_list->empty = true;
67 322 : return;
68 : } else {
69 16846 : task_list->empty = false;
70 : }
71 :
72 16846 : size_t size = (size_t)nlevels * 3 * sizeof(int);
73 :
74 16846 : if (task_list->nlevels < nlevels) {
75 10063 : free(task_list->npts_local);
76 10063 : task_list->npts_local = malloc(size);
77 : }
78 :
79 : // Store npts_local for bounds checking and validation.
80 16846 : task_list->nlevels = nlevels;
81 16846 : assert(task_list->npts_local != NULL);
82 16846 : memcpy(task_list->npts_local, npts_local, size);
83 :
84 : // Always create reference backend because it might be needed for validation.
85 16846 : grid_ref_create_task_list(
86 : orthorhombic, ntasks, nlevels, natoms, nkinds, nblocks, block_offsets,
87 : atom_positions, atom_kinds, basis_sets, level_list, iatom_list,
88 : jatom_list, iset_list, jset_list, ipgf_list, jpgf_list, border_mask_list,
89 : block_num_list, radius_list, rab_list, npts_global, npts_local,
90 : shift_local, border_width, dh, dh_inv, &task_list->ref);
91 :
92 : // Create other backend, if selected.
93 16846 : switch (task_list->backend) {
94 : case GRID_BACKEND_REF:
95 : break; // was already created above
96 16824 : case GRID_BACKEND_CPU:
97 16824 : grid_cpu_create_task_list(
98 : orthorhombic, ntasks, nlevels, natoms, nkinds, nblocks, block_offsets,
99 : atom_positions, atom_kinds, basis_sets, level_list, iatom_list,
100 : jatom_list, iset_list, jset_list, ipgf_list, jpgf_list,
101 : border_mask_list, block_num_list, radius_list, rab_list, npts_global,
102 : npts_local, shift_local, border_width, dh, dh_inv, &task_list->cpu);
103 16824 : break;
104 19 : case GRID_BACKEND_DGEMM:
105 19 : grid_dgemm_create_task_list(
106 : orthorhombic, ntasks, nlevels, natoms, nkinds, nblocks, block_offsets,
107 : atom_positions, atom_kinds, basis_sets, level_list, iatom_list,
108 : jatom_list, iset_list, jset_list, ipgf_list, jpgf_list,
109 : border_mask_list, block_num_list, radius_list, rab_list, npts_global,
110 : npts_local, shift_local, border_width, dh, dh_inv, &task_list->dgemm);
111 19 : break;
112 :
113 0 : case GRID_BACKEND_GPU:
114 : #if (defined(__OFFLOAD_CUDA) || defined(__OFFLOAD_HIP)) && \
115 : !defined(__NO_OFFLOAD_GRID)
116 : grid_gpu_create_task_list(
117 : orthorhombic, ntasks, nlevels, natoms, nkinds, nblocks, block_offsets,
118 : &atom_positions[0][0], atom_kinds, basis_sets, level_list, iatom_list,
119 : jatom_list, iset_list, jset_list, ipgf_list, jpgf_list,
120 : border_mask_list, block_num_list, radius_list, &rab_list[0][0],
121 : &npts_global[0][0], &npts_local[0][0], &shift_local[0][0],
122 : &border_width[0][0], &dh[0][0][0], &dh_inv[0][0][0], &task_list->gpu);
123 : #else
124 0 : fprintf(stderr, "Error: The GPU grid backend is not available. "
125 : "Please re-compile with -D__OFFLOAD");
126 0 : abort();
127 : #endif
128 0 : break;
129 :
130 0 : default:
131 0 : printf("Error: Unknown grid backend: %i.\n", config.backend);
132 0 : abort();
133 16846 : break;
134 : }
135 : }
136 :
137 : /*******************************************************************************
138 : * \brief Deallocates given task list, basis_sets have to be freed separately.
139 : * \author Ole Schuett
140 : ******************************************************************************/
141 10334 : void grid_free_task_list(grid_task_list *ptr) {
142 10334 : if (ptr == NULL)
143 : return;
144 :
145 10334 : grid_task_list_internal *task_list = (grid_task_list_internal *)ptr;
146 :
147 10334 : if (task_list->ref != NULL) {
148 10063 : grid_ref_free_task_list(task_list->ref);
149 10063 : task_list->ref = NULL;
150 : }
151 10334 : if (task_list->cpu != NULL) {
152 10053 : grid_cpu_free_task_list(task_list->cpu);
153 10053 : task_list->cpu = NULL;
154 : }
155 10334 : if (task_list->dgemm != NULL) {
156 7 : grid_dgemm_free_task_list(task_list->dgemm);
157 7 : task_list->dgemm = NULL;
158 : }
159 : #if (defined(__OFFLOAD_CUDA) || defined(__OFFLOAD_HIP)) && \
160 : !defined(__NO_OFFLOAD_GRID)
161 : if (task_list->gpu != NULL) {
162 : grid_gpu_free_task_list(task_list->gpu);
163 : task_list->gpu = NULL;
164 : }
165 : #endif
166 :
167 10334 : free(task_list->npts_local);
168 10334 : free(task_list);
169 : }
170 :
171 : /*******************************************************************************
172 : * \brief Collocate all tasks of in given list onto given grids.
173 : * See grid_task_list.h for details.
174 : * \author Ole Schuett
175 : ******************************************************************************/
176 260376 : void grid_collocate_task_list(const grid_task_list *ptr,
177 : const enum grid_func func, const int nlevels,
178 : const int npts_local[nlevels][3],
179 : const offload_buffer *pab_blocks,
180 : offload_buffer *grids[nlevels]) {
181 260376 : if (ptr == NULL)
182 : return;
183 :
184 260376 : grid_task_list_internal *task_list = (grid_task_list_internal *)ptr;
185 :
186 260376 : if (task_list->empty) {
187 22006 : for (int level = 0; level < nlevels; level++)
188 17554 : memset(grids[level]->host_buffer, 0, grids[level]->size);
189 : return;
190 : }
191 :
192 : // Bounds check.
193 255924 : assert(task_list->nlevels == nlevels);
194 1268378 : for (int ilevel = 0; ilevel < nlevels; ilevel++) {
195 1012454 : assert(task_list->npts_local[ilevel][0] == npts_local[ilevel][0]);
196 1012454 : assert(task_list->npts_local[ilevel][1] == npts_local[ilevel][1]);
197 1012454 : assert(task_list->npts_local[ilevel][2] == npts_local[ilevel][2]);
198 : }
199 :
200 255924 : switch (task_list->backend) {
201 21 : case GRID_BACKEND_REF:
202 21 : grid_ref_collocate_task_list(task_list->ref, func, nlevels, pab_blocks,
203 : grids);
204 21 : break;
205 255750 : case GRID_BACKEND_CPU:
206 255750 : grid_cpu_collocate_task_list(task_list->cpu, func, nlevels, pab_blocks,
207 : grids);
208 255750 : break;
209 153 : case GRID_BACKEND_DGEMM:
210 153 : grid_dgemm_collocate_task_list(task_list->dgemm, func, nlevels, pab_blocks,
211 : grids);
212 153 : break;
213 : #if (defined(__OFFLOAD_CUDA) || defined(__OFFLOAD_HIP)) && \
214 : !defined(__NO_OFFLOAD_GRID)
215 : case GRID_BACKEND_GPU:
216 : grid_gpu_collocate_task_list(task_list->gpu, func, nlevels, pab_blocks,
217 : grids);
218 : break;
219 : #endif
220 0 : default:
221 0 : printf("Error: Unknown grid backend: %i.\n", task_list->backend);
222 0 : abort();
223 255924 : break;
224 : }
225 :
226 : // Perform validation if enabled.
227 255924 : if (grid_library_get_config().validate) {
228 : // Allocate space for reference results.
229 21 : offload_buffer *grids_ref[nlevels];
230 105 : for (int level = 0; level < nlevels; level++) {
231 84 : const int npts_local_total =
232 84 : npts_local[level][0] * npts_local[level][1] * npts_local[level][2];
233 84 : grids_ref[level] = NULL;
234 84 : offload_create_buffer(npts_local_total, &grids_ref[level]);
235 : }
236 :
237 : // Call reference implementation.
238 21 : grid_ref_collocate_task_list(task_list->ref, func, nlevels, pab_blocks,
239 : grids_ref);
240 :
241 : // Compare results.
242 21 : const double tolerance = 1e-12;
243 21 : double max_rel_diff = 0.0;
244 105 : for (int level = 0; level < nlevels; level++) {
245 1491 : for (int i = 0; i < npts_local[level][0]; i++) {
246 43692 : for (int j = 0; j < npts_local[level][1]; j++) {
247 1914354 : for (int k = 0; k < npts_local[level][2]; k++) {
248 1872069 : const int idx = k * npts_local[level][1] * npts_local[level][0] +
249 1872069 : j * npts_local[level][0] + i;
250 1872069 : const double ref_value = grids_ref[level]->host_buffer[idx];
251 1872069 : const double test_value = grids[level]->host_buffer[idx];
252 1872069 : const double diff = fabs(test_value - ref_value);
253 1872069 : const double rel_diff = diff / fmax(1.0, fabs(ref_value));
254 1872069 : max_rel_diff = fmax(max_rel_diff, rel_diff);
255 1872069 : if (rel_diff > tolerance) {
256 0 : fprintf(stderr, "Error: Validation failure in grid collocate\n");
257 0 : fprintf(stderr, " diff: %le\n", diff);
258 0 : fprintf(stderr, " rel_diff: %le\n", rel_diff);
259 0 : fprintf(stderr, " value: %le\n", ref_value);
260 0 : fprintf(stderr, " level: %i\n", level);
261 0 : fprintf(stderr, " ijk: %i %i %i\n", i, j, k);
262 0 : abort();
263 : }
264 : }
265 : }
266 : }
267 84 : offload_free_buffer(grids_ref[level]);
268 84 : printf("Validated grid collocate, max rel. diff: %le\n", max_rel_diff);
269 : }
270 : }
271 : }
272 :
273 : /*******************************************************************************
274 : * \brief Integrate all tasks of in given list from given grids.
275 : * See grid_task_list.h for details.
276 : * \author Ole Schuett
277 : ******************************************************************************/
278 237320 : void grid_integrate_task_list(const grid_task_list *ptr, const bool compute_tau,
279 : const int natoms, const int nlevels,
280 : const int npts_local[nlevels][3],
281 : const offload_buffer *pab_blocks,
282 : const offload_buffer *grids[nlevels],
283 : offload_buffer *hab_blocks,
284 : double forces[natoms][3], double virial[3][3]) {
285 :
286 237320 : if (ptr == NULL)
287 : return;
288 :
289 237320 : grid_task_list_internal *task_list = (grid_task_list_internal *)ptr;
290 :
291 237320 : if (task_list->empty) {
292 3353 : memset(hab_blocks->host_buffer, 0, hab_blocks->size);
293 3353 : if (virial) {
294 29 : virial[0][0] = 0.0;
295 29 : virial[0][1] = 0.0;
296 29 : virial[0][2] = 0.0;
297 29 : virial[1][0] = 0.0;
298 29 : virial[1][1] = 0.0;
299 29 : virial[1][2] = 0.0;
300 29 : virial[2][0] = 0.0;
301 29 : virial[2][1] = 0.0;
302 29 : virial[2][2] = 0.0;
303 : }
304 3353 : if (forces) {
305 342 : for (int atom = 0; atom < natoms; atom++) {
306 213 : forces[atom][0] = 0.0;
307 213 : forces[atom][1] = 0.0;
308 213 : forces[atom][2] = 0.0;
309 : }
310 : }
311 3353 : return;
312 : }
313 :
314 : // Bounds check.
315 233967 : assert(task_list->nlevels == nlevels);
316 1160655 : for (int ilevel = 0; ilevel < nlevels; ilevel++) {
317 926688 : assert(task_list->npts_local[ilevel][0] == npts_local[ilevel][0]);
318 926688 : assert(task_list->npts_local[ilevel][1] == npts_local[ilevel][1]);
319 926688 : assert(task_list->npts_local[ilevel][2] == npts_local[ilevel][2]);
320 : }
321 :
322 233967 : assert(forces == NULL || pab_blocks != NULL);
323 233967 : assert(virial == NULL || pab_blocks != NULL);
324 :
325 233967 : switch (task_list->backend) {
326 : #if (defined(__OFFLOAD_CUDA) || defined(__OFFLOAD_HIP)) && \
327 : !defined(__NO_OFFLOAD_GRID)
328 : case GRID_BACKEND_GPU:
329 : grid_gpu_integrate_task_list(task_list->gpu, compute_tau, nlevels,
330 : pab_blocks, grids, hab_blocks, &forces[0][0],
331 : &virial[0][0]);
332 : break;
333 : #endif
334 149 : case GRID_BACKEND_DGEMM:
335 149 : grid_dgemm_integrate_task_list(task_list->dgemm, compute_tau, natoms,
336 : nlevels, pab_blocks, grids, hab_blocks,
337 : forces, virial);
338 149 : break;
339 233801 : case GRID_BACKEND_CPU:
340 233801 : grid_cpu_integrate_task_list(task_list->cpu, compute_tau, natoms, nlevels,
341 : pab_blocks, grids, hab_blocks, forces, virial);
342 233801 : break;
343 17 : case GRID_BACKEND_REF:
344 17 : grid_ref_integrate_task_list(task_list->ref, compute_tau, natoms, nlevels,
345 : pab_blocks, grids, hab_blocks, forces, virial);
346 17 : break;
347 0 : default:
348 0 : printf("Error: Unknown grid backend: %i.\n", task_list->backend);
349 0 : abort();
350 233967 : break;
351 : }
352 :
353 : // Perform validation if enabled.
354 233967 : if (grid_library_get_config().validate) {
355 : // Allocate space for reference results.
356 17 : const int hab_length = hab_blocks->size / sizeof(double);
357 17 : offload_buffer *hab_blocks_ref = NULL;
358 17 : offload_create_buffer(hab_length, &hab_blocks_ref);
359 17 : double forces_ref[natoms][3], virial_ref[3][3];
360 :
361 : // Call reference implementation.
362 34 : grid_ref_integrate_task_list(task_list->ref, compute_tau, natoms, nlevels,
363 : pab_blocks, grids, hab_blocks_ref,
364 : (forces != NULL) ? forces_ref : NULL,
365 : (virial != NULL) ? virial_ref : NULL);
366 :
367 : // Compare hab.
368 17 : const double hab_tolerance = 1e-12;
369 17 : double hab_max_rel_diff = 0.0;
370 2178 : for (int i = 0; i < hab_length; i++) {
371 2161 : const double ref_value = hab_blocks_ref->host_buffer[i];
372 2161 : const double test_value = hab_blocks->host_buffer[i];
373 2161 : const double diff = fabs(test_value - ref_value);
374 2161 : const double rel_diff = diff / fmax(1.0, fabs(ref_value));
375 2161 : hab_max_rel_diff = fmax(hab_max_rel_diff, rel_diff);
376 2161 : if (rel_diff > hab_tolerance) {
377 0 : fprintf(stderr, "Error: Validation failure in grid integrate\n");
378 0 : fprintf(stderr, " hab diff: %le\n", diff);
379 0 : fprintf(stderr, " hab rel_diff: %le\n", rel_diff);
380 0 : fprintf(stderr, " hab value: %le\n", ref_value);
381 0 : fprintf(stderr, " hab i: %i\n", i);
382 0 : abort();
383 : }
384 : }
385 :
386 : // Compare forces.
387 17 : const double forces_tolerance = 1e-8; // account for higher numeric noise
388 17 : double forces_max_rel_diff = 0.0;
389 17 : if (forces != NULL) {
390 12 : for (int iatom = 0; iatom < natoms; iatom++) {
391 36 : for (int idir = 0; idir < 3; idir++) {
392 27 : const double ref_value = forces_ref[iatom][idir];
393 27 : const double test_value = forces[iatom][idir];
394 27 : const double diff = fabs(test_value - ref_value);
395 27 : const double rel_diff = diff / fmax(1.0, fabs(ref_value));
396 27 : forces_max_rel_diff = fmax(forces_max_rel_diff, rel_diff);
397 27 : if (rel_diff > forces_tolerance) {
398 0 : fprintf(stderr, "Error: Validation failure in grid integrate\n");
399 0 : fprintf(stderr, " forces diff: %le\n", diff);
400 0 : fprintf(stderr, " forces rel_diff: %le\n", rel_diff);
401 0 : fprintf(stderr, " forces value: %le\n", ref_value);
402 0 : fprintf(stderr, " forces atom: %i\n", iatom);
403 0 : fprintf(stderr, " forces dir: %i\n", idir);
404 0 : abort();
405 : }
406 : }
407 : }
408 : }
409 :
410 : // Compare virial.
411 17 : const double virial_tolerance = 1e-8; // account for higher numeric noise
412 17 : double virial_max_rel_diff = 0.0;
413 17 : if (virial != NULL) {
414 0 : for (int i = 0; i < 3; i++) {
415 0 : for (int j = 0; j < 3; j++) {
416 0 : const double ref_value = virial_ref[i][j];
417 0 : const double test_value = virial[i][j];
418 0 : const double diff = fabs(test_value - ref_value);
419 0 : const double rel_diff = diff / fmax(1.0, fabs(ref_value));
420 0 : virial_max_rel_diff = fmax(virial_max_rel_diff, rel_diff);
421 0 : if (rel_diff > virial_tolerance) {
422 0 : fprintf(stderr, "Error: Validation failure in grid integrate\n");
423 0 : fprintf(stderr, " virial diff: %le\n", diff);
424 0 : fprintf(stderr, " virial rel_diff: %le\n", rel_diff);
425 0 : fprintf(stderr, " virial value: %le\n", ref_value);
426 0 : fprintf(stderr, " virial ij: %i %i\n", i, j);
427 0 : abort();
428 : }
429 : }
430 : }
431 : }
432 :
433 17 : printf("Validated grid_integrate, max rel. diff: %le %le %le\n",
434 : hab_max_rel_diff, forces_max_rel_diff, virial_max_rel_diff);
435 17 : offload_free_buffer(hab_blocks_ref);
436 : }
437 : }
438 :
439 : // EOF
|
