Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2025 CP2K developers group <https://cp2k.org> */
4 : /* */
5 : /* SPDX-License-Identifier: BSD-3-Clause */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : #include <assert.h>
9 : #include <stdbool.h>
10 : #include <stdio.h>
11 : #include <stdlib.h>
12 : #include <string.h>
13 : #include <unistd.h>
14 :
15 : #include "../common/grid_common.h"
16 : #include "grid_dgemm_collocation_integration.h"
17 : #include "grid_dgemm_non_orthorombic_corrections.h"
18 : #include "grid_dgemm_tensor_local.h"
19 : #include "grid_dgemm_utils.h"
20 :
21 8 : struct collocation_integration_ *collocate_create_handle(void) {
22 8 : struct collocation_integration_ *handle =
23 8 : malloc(sizeof(struct collocation_integration_));
24 8 : assert(handle != NULL);
25 8 : memset(handle, 0, sizeof(struct collocation_integration_));
26 :
27 8 : handle->alpha.alloc_size_ = 8192;
28 8 : handle->coef.alloc_size_ = 1024;
29 8 : handle->pol.alloc_size_ = 1024;
30 : /* it is a cube of size 32 x 32 x 32 */
31 8 : handle->cube.alloc_size_ = 32768;
32 :
33 8 : handle->cube_alloc_size = realloc_tensor(&handle->cube);
34 8 : handle->alpha_alloc_size = realloc_tensor(&handle->alpha);
35 8 : handle->coef_alloc_size = realloc_tensor(&handle->coef);
36 8 : handle->pol_alloc_size = realloc_tensor(&handle->pol);
37 :
38 8 : handle->scratch = malloc(32768 * sizeof(double));
39 8 : assert(handle->scratch != NULL);
40 8 : handle->scratch_alloc_size = 32768;
41 8 : handle->T_alloc_size = 8192;
42 8 : handle->W_alloc_size = 2048;
43 8 : handle->blockDim[0] = 5;
44 8 : handle->blockDim[1] = 5;
45 8 : handle->blockDim[2] = 5;
46 8 : handle->device_id = (int *)malloc(sizeof(double) * 12);
47 8 : assert(handle->device_id != NULL);
48 8 : handle->number_of_devices = 1;
49 :
50 : /* to suppress when we remove the spherical cutoff */
51 8 : handle->map = (int **)malloc(3 * sizeof(int *));
52 8 : assert(handle->map != NULL);
53 8 : handle->map[0] = (int *)malloc(sizeof(int) * 512 * 3);
54 8 : assert(handle->map[0] != NULL);
55 8 : handle->map[1] = handle->map[0] + 512;
56 8 : handle->map[2] = handle->map[1] + 512;
57 8 : handle->cmax = 512 * 3;
58 8 : return handle;
59 : }
60 :
61 8 : void collocate_destroy_handle(void *gaussian_handle) {
62 8 : struct collocation_integration_ *handle =
63 : (struct collocation_integration_ *)gaussian_handle;
64 8 : if (handle->Exp.data)
65 4 : free(handle->Exp.data);
66 :
67 8 : if (handle->grid.data)
68 0 : free(handle->grid.data);
69 :
70 8 : free(handle->scratch);
71 8 : free(handle->pol.data);
72 8 : free(handle->cube.data);
73 8 : free(handle->alpha.data);
74 8 : free(handle->coef.data);
75 8 : free(handle->blocks_coordinates.data);
76 8 : handle->alpha.data = NULL;
77 8 : handle->coef.data = NULL;
78 8 : handle->blocks_coordinates.data = NULL;
79 8 : free(handle->device_id);
80 8 : free(handle->map[0]);
81 8 : free(handle->map);
82 8 : free(handle);
83 :
84 8 : handle = NULL;
85 8 : }
86 :
87 90334 : void initialize_W_and_T(collocation_integration *const handler,
88 : const tensor *cube, const tensor *coef) {
89 90334 : size_t tmp1 = compute_memory_space_tensor_3(coef->size[0] /* alpha */,
90 90334 : coef->size[1] /* gamma */,
91 90334 : cube->size[1] /* j */);
92 :
93 90334 : size_t tmp2 = compute_memory_space_tensor_3(
94 90334 : coef->size[0] /* gamma */, cube->size[1] /* j */, cube->size[2] /* i */);
95 :
96 90334 : const size_t mem_alloc_size_ =
97 90334 : imax(imax(tmp1 + tmp2, cube->alloc_size_), coef->alloc_size_);
98 :
99 90334 : handler->T_alloc_size = tmp1;
100 90334 : handler->W_alloc_size = tmp2;
101 :
102 90334 : if ((mem_alloc_size_ > handler->scratch_alloc_size) ||
103 90315 : (handler->scratch == NULL)) {
104 19 : handler->scratch_alloc_size = mem_alloc_size_;
105 :
106 19 : if (handler->scratch)
107 19 : free(handler->scratch);
108 19 : handler->scratch = malloc(sizeof(double) * handler->scratch_alloc_size);
109 19 : assert(handler->scratch != NULL);
110 : }
111 90334 : }
112 :
113 0 : void initialize_W_and_T_integrate(collocation_integration *const handler,
114 : const int num_block, const tensor *coef,
115 : const tensor *block) {
116 : /* T */
117 0 : size_t tmp1 = compute_memory_space_tensor_4(num_block, block->size[0] /* k */,
118 0 : block->size[1] /* j */,
119 0 : coef->size[1] /* alpha */);
120 :
121 : /* W */
122 0 : size_t tmp2 = compute_memory_space_tensor_4(num_block, block->size[1] /* j */,
123 : coef->size[1] /* alpha */,
124 0 : coef->size[2] /* gamma */);
125 :
126 0 : const size_t mem_alloc_size_ = tmp1 + tmp2;
127 :
128 0 : handler->T_alloc_size = tmp1;
129 0 : handler->W_alloc_size = tmp2;
130 :
131 0 : if ((mem_alloc_size_ > handler->scratch_alloc_size) ||
132 0 : (handler->scratch == NULL)) {
133 0 : handler->scratch_alloc_size = mem_alloc_size_;
134 :
135 0 : if (handler->scratch)
136 0 : free(handler->scratch);
137 0 : handler->scratch = malloc(sizeof(double) * handler->scratch_alloc_size);
138 0 : assert(handler->scratch != NULL);
139 : }
140 0 : }
141 :
142 1280 : void initialize_basis_vectors(collocation_integration *const handler,
143 : const double dh[3][3],
144 : const double dh_inv[3][3]) {
145 1280 : handler->dh[0][0] = dh[0][0];
146 1280 : handler->dh[0][1] = dh[0][1];
147 1280 : handler->dh[0][2] = dh[0][2];
148 1280 : handler->dh[1][0] = dh[1][0];
149 1280 : handler->dh[1][1] = dh[1][1];
150 1280 : handler->dh[1][2] = dh[1][2];
151 1280 : handler->dh[2][0] = dh[2][0];
152 1280 : handler->dh[2][1] = dh[2][1];
153 1280 : handler->dh[2][2] = dh[2][2];
154 :
155 1280 : handler->dh_inv[0][0] = dh_inv[0][0];
156 1280 : handler->dh_inv[0][1] = dh_inv[0][1];
157 1280 : handler->dh_inv[0][2] = dh_inv[0][2];
158 1280 : handler->dh_inv[1][0] = dh_inv[1][0];
159 1280 : handler->dh_inv[1][1] = dh_inv[1][1];
160 1280 : handler->dh_inv[1][2] = dh_inv[1][2];
161 1280 : handler->dh_inv[2][0] = dh_inv[2][0];
162 1280 : handler->dh_inv[2][1] = dh_inv[2][1];
163 1280 : handler->dh_inv[2][2] = dh_inv[2][2];
164 :
165 : /* Only used when we are in the non orthorombic case */
166 1280 : handler->dx[2] = handler->dh[0][0] * handler->dh[0][0] +
167 1280 : handler->dh[0][1] * handler->dh[0][1] +
168 1280 : handler->dh[0][2] * handler->dh[0][2];
169 1280 : handler->dx[1] = handler->dh[1][0] * handler->dh[1][0] +
170 1280 : handler->dh[1][1] * handler->dh[1][1] +
171 1280 : handler->dh[1][2] * handler->dh[1][2];
172 1280 : handler->dx[0] = handler->dh[2][0] * handler->dh[2][0] +
173 1280 : handler->dh[2][1] * handler->dh[2][1] +
174 1280 : handler->dh[2][2] * handler->dh[2][2];
175 1280 : }
|
