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 <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 = NULL;
      23           8 :   handle = (struct collocation_integration_ *)malloc(
      24             :       sizeof(struct collocation_integration_));
      25           8 :   assert(handle != NULL);
      26           8 :   memset(handle, 0, sizeof(struct collocation_integration_));
      27             : 
      28           8 :   handle->alpha.alloc_size_ = 8192;
      29           8 :   handle->coef.alloc_size_ = 1024;
      30           8 :   handle->pol.alloc_size_ = 1024;
      31             :   /* it is a cube of size 32 x 32 x 32 */
      32           8 :   handle->cube.alloc_size_ = 32768;
      33             : 
      34           8 :   handle->cube_alloc_size = realloc_tensor(&handle->cube);
      35           8 :   handle->alpha_alloc_size = realloc_tensor(&handle->alpha);
      36           8 :   handle->coef_alloc_size = realloc_tensor(&handle->coef);
      37           8 :   handle->pol_alloc_size = realloc_tensor(&handle->pol);
      38             : 
      39           8 :   handle->scratch = malloc(32768 * sizeof(double));
      40           8 :   assert(handle->scratch != NULL);
      41           8 :   handle->scratch_alloc_size = 32768;
      42           8 :   handle->T_alloc_size = 8192;
      43           8 :   handle->W_alloc_size = 2048;
      44           8 :   handle->blockDim[0] = 5;
      45           8 :   handle->blockDim[1] = 5;
      46           8 :   handle->blockDim[2] = 5;
      47           8 :   handle->device_id = (int *)malloc(sizeof(double) * 12);
      48           8 :   assert(handle->device_id != NULL);
      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 :   assert(handle->map != NULL);
      54           8 :   handle->map[0] = (int *)malloc(sizeof(int) * 512 * 3);
      55           8 :   assert(handle->map[0] != NULL);
      56           8 :   handle->map[1] = handle->map[0] + 512;
      57           8 :   handle->map[2] = handle->map[1] + 512;
      58           8 :   handle->cmax = 512 * 3;
      59           8 :   return handle;
      60             : }
      61             : 
      62           8 : void collocate_destroy_handle(void *gaussian_handle) {
      63           8 :   struct collocation_integration_ *handle =
      64             :       (struct collocation_integration_ *)gaussian_handle;
      65           8 :   if (handle->Exp.data)
      66           4 :     free(handle->Exp.data);
      67             : 
      68           8 :   if (handle->grid.data)
      69           0 :     free(handle->grid.data);
      70             : 
      71           8 :   free(handle->scratch);
      72           8 :   free(handle->pol.data);
      73           8 :   free(handle->cube.data);
      74           8 :   free(handle->alpha.data);
      75           8 :   free(handle->coef.data);
      76           8 :   free(handle->blocks_coordinates.data);
      77           8 :   handle->alpha.data = NULL;
      78           8 :   handle->coef.data = NULL;
      79           8 :   handle->blocks_coordinates.data = NULL;
      80           8 :   free(handle->device_id);
      81           8 :   free(handle->map[0]);
      82           8 :   free(handle->map);
      83           8 :   free(handle);
      84             : 
      85           8 :   handle = NULL;
      86           8 : }
      87             : 
      88       88000 : void initialize_W_and_T(collocation_integration *const handler,
      89             :                         const tensor *cube, const tensor *coef) {
      90       88000 :   size_t tmp1 = compute_memory_space_tensor_3(coef->size[0] /* alpha */,
      91       88000 :                                               coef->size[1] /* gamma */,
      92       88000 :                                               cube->size[1] /* j */);
      93             : 
      94       88000 :   size_t tmp2 = compute_memory_space_tensor_3(
      95       88000 :       coef->size[0] /* gamma */, cube->size[1] /* j */, cube->size[2] /* i */);
      96             : 
      97       88000 :   const size_t mem_alloc_size_ =
      98       88000 :       imax(imax(tmp1 + tmp2, cube->alloc_size_), coef->alloc_size_);
      99             : 
     100       88000 :   handler->T_alloc_size = tmp1;
     101       88000 :   handler->W_alloc_size = tmp2;
     102             : 
     103       88000 :   if ((mem_alloc_size_ > handler->scratch_alloc_size) ||
     104       87981 :       (handler->scratch == NULL)) {
     105          19 :     handler->scratch_alloc_size = mem_alloc_size_;
     106             : 
     107          19 :     if (handler->scratch)
     108          19 :       free(handler->scratch);
     109          19 :     handler->scratch = malloc(sizeof(double) * handler->scratch_alloc_size);
     110          19 :     assert(handler->scratch != NULL);
     111             :   }
     112       88000 : }
     113             : 
     114           0 : void initialize_W_and_T_integrate(collocation_integration *const handler,
     115             :                                   const int num_block, const tensor *coef,
     116             :                                   const tensor *block) {
     117             :   /* T */
     118           0 :   size_t tmp1 = compute_memory_space_tensor_4(num_block, block->size[0] /* k */,
     119           0 :                                               block->size[1] /* j */,
     120           0 :                                               coef->size[1] /* alpha */);
     121             : 
     122             :   /* W */
     123           0 :   size_t tmp2 = compute_memory_space_tensor_4(num_block, block->size[1] /* j */,
     124             :                                               coef->size[1] /* alpha */,
     125           0 :                                               coef->size[2] /* gamma */);
     126             : 
     127           0 :   const size_t mem_alloc_size_ = tmp1 + tmp2;
     128             : 
     129           0 :   handler->T_alloc_size = tmp1;
     130           0 :   handler->W_alloc_size = tmp2;
     131             : 
     132           0 :   if ((mem_alloc_size_ > handler->scratch_alloc_size) ||
     133           0 :       (handler->scratch == NULL)) {
     134           0 :     handler->scratch_alloc_size = mem_alloc_size_;
     135             : 
     136           0 :     if (handler->scratch)
     137           0 :       free(handler->scratch);
     138           0 :     handler->scratch = malloc(sizeof(double) * handler->scratch_alloc_size);
     139           0 :     assert(handler->scratch != NULL);
     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 : }