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     7671796 : inline int int_max2(int x, int y) { return (x > y) ? x : y; }
      20             : 
      21      149632 : inline int int_max3(int x, int y, int z) { return int_max2(int_max2(x, y), z); }
      22             : 
      23     1008243 : double *alloc_zeros_1d(int n) { return (double *)calloc(n, sizeof(double)); }
      24             : 
      25      413121 : double **alloc_zeros_2d(int n, int m) {
      26      413121 :   double **array = (double **)calloc(n, sizeof(double *));
      27             : 
      28     1882110 :   for (int i = 0; i < n; i++) {
      29     1468989 :     array[i] = (double *)calloc(m, sizeof(double));
      30             :   }
      31             : 
      32      413121 :   return array;
      33             : }
      34             : 
      35      413121 : void free_2d(double **array, int n) {
      36     1882110 :   for (int i = 0; i < n; i++) {
      37     1468989 :     free(array[i]);
      38             :   }
      39      413121 :   free(array);
      40      413121 : }
      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             : }