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