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 : }