Skip to content

Commit

Permalink
grid: Simplify prepare_pab and ready it for Cuda
Browse files Browse the repository at this point in the history
  • Loading branch information
@oschuett
oschuett committed Nov 15, 2020
1 parent 76181cb commit d8fd13a
Show file tree
Hide file tree
Showing 8 changed files with 563 additions and 1,067 deletions.
120 changes: 85 additions & 35 deletions src/grid/common/grid_common.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,73 +7,123 @@
#ifndef GRID_COMMON_H
#define GRID_COMMON_H

#if defined(__CUDACC__)
#define GRID_HOST_DEVICE __host__ __device__
#else
#define GRID_HOST_DEVICE
#endif

/*******************************************************************************
* \brief Number of Cartesian orbitals up to given angular momentum quantum.
* \brief Factorial function, e.g. fac(5) = 5! = 120.
* \author Ole Schuett
******************************************************************************/
static const int ncoset[] = {1, // l=0
4, // l=1
10, // l=2 ...
20, 35, 56, 84, 120, 165, 220, 286,
364, 455, 560, 680, 816, 969, 1140, 1330};
GRID_HOST_DEVICE static inline double fac(const int i) {
static const double table[] = {
0.10000000000000000000E+01, 0.10000000000000000000E+01,
0.20000000000000000000E+01, 0.60000000000000000000E+01,
0.24000000000000000000E+02, 0.12000000000000000000E+03,
0.72000000000000000000E+03, 0.50400000000000000000E+04,
0.40320000000000000000E+05, 0.36288000000000000000E+06,
0.36288000000000000000E+07, 0.39916800000000000000E+08,
0.47900160000000000000E+09, 0.62270208000000000000E+10,
0.87178291200000000000E+11, 0.13076743680000000000E+13,
0.20922789888000000000E+14, 0.35568742809600000000E+15,
0.64023737057280000000E+16, 0.12164510040883200000E+18,
0.24329020081766400000E+19, 0.51090942171709440000E+20,
0.11240007277776076800E+22, 0.25852016738884976640E+23,
0.62044840173323943936E+24, 0.15511210043330985984E+26,
0.40329146112660563558E+27, 0.10888869450418352161E+29,
0.30488834461171386050E+30, 0.88417619937397019545E+31,
0.26525285981219105864E+33};
return table[i];
}

/*******************************************************************************
* \brief Tabulation of the Factorial function, e.g. 5! = fac[5] = 120.
* \brief Number of Cartesian orbitals up to given angular momentum quantum.
* \author Ole Schuett
******************************************************************************/
static const double fac[] = {
0.10000000000000000000E+01, 0.10000000000000000000E+01,
0.20000000000000000000E+01, 0.60000000000000000000E+01,
0.24000000000000000000E+02, 0.12000000000000000000E+03,
0.72000000000000000000E+03, 0.50400000000000000000E+04,
0.40320000000000000000E+05, 0.36288000000000000000E+06,
0.36288000000000000000E+07, 0.39916800000000000000E+08,
0.47900160000000000000E+09, 0.62270208000000000000E+10,
0.87178291200000000000E+11, 0.13076743680000000000E+13,
0.20922789888000000000E+14, 0.35568742809600000000E+15,
0.64023737057280000000E+16, 0.12164510040883200000E+18,
0.24329020081766400000E+19, 0.51090942171709440000E+20,
0.11240007277776076800E+22, 0.25852016738884976640E+23,
0.62044840173323943936E+24, 0.15511210043330985984E+26,
0.40329146112660563558E+27, 0.10888869450418352161E+29,
0.30488834461171386050E+30, 0.88417619937397019545E+31,
0.26525285981219105864E+33};
GRID_HOST_DEVICE static inline int ncoset(const int l) {
static const int table[] = {1, // l=0
4, // l=1
10, // l=2 ...
20, 35, 56, 84, 120, 165, 220, 286,
364, 455, 560, 680, 816, 969, 1140, 1330};
return table[l];
}

/*******************************************************************************
* \brief Maps three angular momentum components to a single zero based index.
* \author Ole Schuett
******************************************************************************/
static inline int coset(int lx, int ly, int lz) {
GRID_HOST_DEVICE static inline int coset(int lx, int ly, int lz) {
const int l = lx + ly + lz;
if (l == 0) {
return 0;
} else {
return ncoset[l - 1] + ((l - lx) * (l - lx + 1)) / 2 + lz;
return ncoset(l - 1) + ((l - lx) * (l - lx + 1)) / 2 + lz;
}
}

/*******************************************************************************
* \brief Minimum and maximum for integers (missing from the C standard)
* \brief Returns the smaller of two given integer (missing from the C standard)
* \author Ole Schuett
******************************************************************************/
GRID_HOST_DEVICE static inline int imin(int x, int y) {
return (x < y ? x : y);
}

/*******************************************************************************
* \brief Returns the larger of two given integer (missing from the C standard)
* \author Ole Schuett
******************************************************************************/
static inline int imin(int x, int y) { return (x < y ? x : y); }
static inline int imax(int x, int y) { return (x > y ? x : y); }
GRID_HOST_DEVICE static inline int imax(int x, int y) {
return (x > y ? x : y);
}

/*******************************************************************************
* \brief Equivalent of Fortran's MODULO, which always return a positive number.
* https://gcc.gnu.org/onlinedocs/gfortran/MODULO.html
* \author Ole Schuett
******************************************************************************/
static inline int modulo(int a, int m) { return ((a % m + m) % m); }
GRID_HOST_DEVICE static inline int modulo(int a, int m) {
return ((a % m + m) % m);
}

/*******************************************************************************
* \brief Prototype for BLAS dgemm.
* \brief Orbital angular momentum.
* \author Ole Schuett
******************************************************************************/
void dgemm_(const char *transa, const char *transb, const int *m, const int *n,
const int *k, const double *alpha, const double *a, const int *lda,
const double *b, const int *ldb, const double *beta, double *c,
const int *ldc);
typedef struct {
int l[3];
} orbital;

/*******************************************************************************
* \brief Increase i'th component of given orbital angular momentum.
* \author Ole Schuett
******************************************************************************/
GRID_HOST_DEVICE static inline orbital up(const int i, const orbital a) {
orbital b = a;
b.l[i] += 1;
return b;
}

/*******************************************************************************
* \brief Decrease i'th component of given orbital angular momentum.
* \author Ole Schuett
******************************************************************************/
GRID_HOST_DEVICE static inline orbital down(const int i, const orbital a) {
orbital b = a;
b.l[i] = imax(0, a.l[i] - 1);
return b;
}

/*******************************************************************************
* \brief Return coset index of given orbital angular momentum.
* \author Ole Schuett
******************************************************************************/
GRID_HOST_DEVICE static inline int idx(const orbital a) {
return coset(a.l[0], a.l[1], a.l[2]);
}

#endif // GRID_COMMON_H

Expand Down

0 comments on commit d8fd13a

Please sign in to comment.