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: MIT                                              */
       6              : /*----------------------------------------------------------------------------*/
       7              : 
       8              : /*
       9              :  *  libgrpp - a library for the evaluation of integrals over
      10              :  *            generalized relativistic pseudopotentials.
      11              :  *
      12              :  *  Copyright (C) 2021-2023 Alexander Oleynichenko
      13              :  */
      14              : #include <math.h>
      15              : #include <stdlib.h>
      16              : 
      17              : #include "grpp_utils.h"
      18              : 
      19      9088126 : inline int int_max2(int x, int y) { return (x > y) ? x : y; }
      20              : 
      21       194176 : inline int int_max3(int x, int y, int z) { return int_max2(int_max2(x, y), z); }
      22              : 
      23      1202574 : double *alloc_zeros_1d(int n) { return (double *)calloc(n, sizeof(double)); }
      24              : 
      25       469812 : double **alloc_zeros_2d(int n, int m) {
      26       469812 :   double **array = (double **)calloc(n, sizeof(double *));
      27              : 
      28      2088882 :   for (int i = 0; i < n; i++) {
      29      1619070 :     array[i] = (double *)calloc(m, sizeof(double));
      30              :   }
      31              : 
      32       469812 :   return array;
      33              : }
      34              : 
      35       469812 : void free_2d(double **array, int n) {
      36      2088882 :   for (int i = 0; i < n; i++) {
      37      1619070 :     free(array[i]);
      38              :   }
      39       469812 :   free(array);
      40       469812 : }
      41              : 
      42              : /*
      43              :  * constant times a vector plus a vector:
      44              :  * y = a * x + y
      45              :  */
      46            0 : void libgrpp_daxpy(int n, double a, double *x, double *y) {
      47            0 :   for (int i = 0; i < n; i++) {
      48            0 :     y[i] += a * x[i];
      49              :   }
      50            0 : }
      51              : 
      52              : /*
      53              :  * naive matrix multiplication
      54              :  */
      55            0 : void libgrpp_multiply_matrices(int M, int N, int K, double *A, double *B,
      56              :                                double *C) {
      57            0 :   for (int i = 0; i < M; i++) {
      58            0 :     for (int j = 0; j < N; j++) {
      59              :       double sum = 0.0;
      60            0 :       for (int k = 0; k < K; k++) {
      61            0 :         sum += A[i * K + k] * B[k * N + j];
      62              :       }
      63            0 :       C[i * N + j] += sum;
      64              :     }
      65              :   }
      66            0 : }
      67              : 
      68            0 : double distance_squared(double *A, double *B) {
      69            0 :   double dx = A[0] - B[0];
      70            0 :   double dy = A[1] - B[1];
      71            0 :   double dz = A[2] - B[2];
      72            0 :   return dx * dx + dy * dy + dz * dz;
      73              : }
      74              : 
      75            0 : double distance(double *A, double *B) { return sqrt(distance_squared(A, B)); }
      76              : 
      77              : /**
      78              :  * Checks if two 3d points coincide with each other.
      79              :  */
      80            0 : int points_are_equal(double *a, double *b) {
      81            0 :   double const thresh = 1e-12;
      82              : 
      83            0 :   if (fabs(a[0] - b[0]) < thresh && fabs(a[1] - b[1]) < thresh &&
      84            0 :       fabs(a[2] - b[2]) < thresh) {
      85            0 :     return 1;
      86              :   }
      87              : 
      88              :   return 0;
      89              : }