TREX/qmckl

mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-06-01 02:45:43 +02:00

Aurélien Delval d0eb207404

Integration of Verificarlo CI tests (#1 )

* comment

* Update distance test code

The distance test has been updated to the latest version, with a first attempt at using vfc_probes inside it

* Functional implementation of vfc_probes in the distance tests

This commit has the first functional vfc_ci tests. Verificarlo tests
should be written over the existing tests, and they can be enabled with
the following configure command:

QMCKL_DEVEL=1 ./configure --prefix=$PWD/_install --enable-maintainer-mode --enable-vfc_ci CC="verificarlo-f -Mpreprocess -D VFC_CI" FC="verificarlo-f -Mpreprocess -D VFC_CI" --host=x86_64

The --enable-vfc_ci flag will trigger the linking of the vfc_ci
library. Moreover, as of now, the "-Mpreprocess" and "-D VFC_CI" flags
have to be specified directly here. There is probably an appropriate
macro to place those flags into but I couldn't find it yet, and could
only manage to build the tests this way.

When the VFC_CI preprocessor is defined, somme additional code to
register and export the test probes can be executed (see
qmckl_distance.org).

As of now, the tests are built as normal, even though they are expected
to fail :

make all
make check

From there, the test_qmckl_distance (and potentially the others)
executable can be called at will. This will typically be done
automatically by vfc_ci, but one could manually execute the executable
by defining the following env variables :

VFC_PROBES_OUTPUT="test.csv" VFC_BACKENDS="libinterflop_ieee.so"

depending on the export file and the Verificarlo backend to be used.

The next steps will be to define more tests such as this one, and to
integrate them into a Verificarlo CI workflow (by writing a
vfc_tests_config.json file and using the automatic CI setup
command).

* Error in FOrtran interface fixed

* Added missing Fortran interfaces

* Modify distance test and install process integration

All probes are now ignored using only the preprocessor (instead
of checking for a facultative argument) in the distance test.
Moreover,preprocessing can now be enabled correctly using FCFLAGS
(the issue seemed to come from the order of the arguments passed
to gfortran/verificarlo-f with the preprocessor arg having to come
first).

* Add vfc_probes to AO tests

vfc_probes have been added to qmckl_ao.org in the same way as
qmckl_distance.org, which means that it can be enabled or disabled at
compile time using the --enable-vfc_ci option.

qmckl_distance.org has been slightly modified with a better indentation,
and configure.ac now adds the "-D VFC_CI" flag to CFLAGS when vfc_ci is
enabled.

Co-authored-by: Anthony Scemama <scemama@irsamc.ups-tlse.fr>

2021-07-07 13:42:42 +02:00

46 KiB

Raw Blame History

Inter-particle distances

Squared distance
- qmckl_distance_sq
Distance
- qmckl_distance
Rescaled Distance
- qmckl_distance_rescaled
Rescaled Distance Derivatives
- qmckl_distance_rescaled_deriv_e

Functions for the computation of distances between particles.

Squared distance

`qmckl_distance_sq`

qmckl_distance_sq computes the matrix of the squared distances between all pairs of points in two sets, one point within each set:

\[ C_{ij} = \sum_{k=1}^3 (A_{k,i}-B_{k,j})^2 \]

qmckl_context	context	in	Global state
char	transa	in	Array `A` is `'N'`: Normal, `'T'`: Transposed
char	transb	in	Array `B` is `'N'`: Normal, `'T'`: Transposed
int64_t	m	in	Number of points in the first set
int64_t	n	in	Number of points in the second set
double	A[][lda]	in	Array containing the $m \times 3$ matrix $A$
int64_t	lda	in	Leading dimension of array `A`
double	B[][ldb]	in	Array containing the $n \times 3$ matrix $B$
int64_t	ldb	in	Leading dimension of array `B`
double	C[n][ldc]	out	Array containing the $m \times n$ matrix $C$
int64_t	ldc	in	Leading dimension of array `C`

Requirements

context is not QMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3 if transa == 'N'
lda >= m if transa == 'T'
ldb >= 3 if transb == 'N'
ldb >= n if transb == 'T'
ldc >= m
A is allocated with at least $3 \times m \times 8$ bytes
B is allocated with at least $3 \times n \times 8$ bytes
C is allocated with at least $m \times n \times 8$ bytes

C header

qmckl_exit_code qmckl_distance_sq (
      const qmckl_context context,
      const char transa,
      const char transb,
      const int64_t m,
      const int64_t n,
      const double* A,
      const int64_t lda,
      const double* B,
      const int64_t ldb,
      double* const C,
      const int64_t ldc );

Source

integer function qmckl_distance_sq_f(context, transa, transb, m, n, &
 A, LDA, B, LDB, C, LDC) &
 result(info)
use qmckl
implicit none
integer(qmckl_context)  , intent(in)  :: context
character  , intent(in)  :: transa, transb
integer*8  , intent(in)  :: m, n
integer*8  , intent(in)  :: lda
real*8     , intent(in)  :: A(lda,*)
integer*8  , intent(in)  :: ldb
real*8     , intent(in)  :: B(ldb,*)
integer*8  , intent(in)  :: ldc
real*8     , intent(out) :: C(ldc,*)

integer*8 :: i,j
real*8    :: x, y, z
integer   :: transab

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
 info = QMCKL_INVALID_CONTEXT
 return
endif

if (m <= 0_8) then
 info = QMCKL_INVALID_ARG_4
 return
endif

if (n <= 0_8) then
 info = QMCKL_INVALID_ARG_5
 return
endif

if (transa == 'N' .or. transa == 'n') then
 transab = 0
else if (transa == 'T' .or. transa == 't') then
 transab = 1
else
 transab = -100
endif

if (transb == 'N' .or. transb == 'n') then
 continue
else if (transa == 'T' .or. transa == 't') then
 transab = transab + 2
else
 transab = -100
endif

if (transab < 0) then
 info = QMCKL_INVALID_ARG_1
 return
endif

if (iand(transab,1) == 0 .and. LDA < 3) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,1) == 1 .and. LDA < m) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,2) == 0 .and. LDA < 3) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,2) == 2 .and. LDA < m) then
 info = QMCKL_INVALID_ARG_7
 return
endif


select case (transab)

case(0)

 do j=1,n
    do i=1,m
       x = A(1,i) - B(1,j)
       y = A(2,i) - B(2,j)
       z = A(3,i) - B(3,j)
       C(i,j) = x*x + y*y + z*z
    end do
 end do

case(1)

 do j=1,n
    do i=1,m
       x = A(i,1) - B(1,j)
       y = A(i,2) - B(2,j)
       z = A(i,3) - B(3,j)
       C(i,j) = x*x + y*y + z*z
    end do
 end do

case(2)

 do j=1,n
    do i=1,m
       x = A(1,i) - B(j,1)
       y = A(2,i) - B(j,2)
       z = A(3,i) - B(j,3)
       C(i,j) = x*x + y*y + z*z
    end do
 end do

case(3)

 do j=1,n
    do i=1,m
       x = A(i,1) - B(j,1)
       y = A(i,2) - B(j,2)
       z = A(i,3) - B(j,3)
       C(i,j) = x*x + y*y + z*z
    end do
 end do

end select

end function qmckl_distance_sq_f

Performance

This function is more efficient when A and B are transposed.

Distance

`qmckl_distance`

qmckl_distance computes the matrix of the distances between all pairs of points in two sets, one point within each set:

\[ C_{ij} = \sqrt{\sum_{k=1}^3 (A_{k,i}-B_{k,j})^2} \]

If the input array is normal ('N'), the xyz coordinates are in the leading dimension: [n][3] in C and (3,n) in Fortran.

qmckl_context	context	in	Global state
char	transa	in	Array `A` is `'N'`: Normal, `'T'`: Transposed
char	transb	in	Array `B` is `'N'`: Normal, `'T'`: Transposed
int64_t	m	in	Number of points in the first set
int64_t	n	in	Number of points in the second set
double	A[][lda]	in	Array containing the $m \times 3$ matrix $A$
int64_t	lda	in	Leading dimension of array `A`
double	B[][ldb]	in	Array containing the $n \times 3$ matrix $B$
int64_t	ldb	in	Leading dimension of array `B`
double	C[n][ldc]	out	Array containing the $m \times n$ matrix $C$
int64_t	ldc	in	Leading dimension of array `C`

Requirements

context is not QMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3 if transa == 'N'
lda >= m if transa == 'T'
ldb >= 3 if transb == 'N'
ldb >= n if transb == 'T'
ldc >= m
A is allocated with at least $3 \times m \times 8$ bytes
B is allocated with at least $3 \times n \times 8$ bytes
C is allocated with at least $m \times n \times 8$ bytes

C header

qmckl_exit_code qmckl_distance (
      const qmckl_context context,
      const char transa,
      const char transb,
      const int64_t m,
      const int64_t n,
      const double* A,
      const int64_t lda,
      const double* B,
      const int64_t ldb,
      double* const C,
      const int64_t ldc );

Source

integer function qmckl_distance_f(context, transa, transb, m, n, &
 A, LDA, B, LDB, C, LDC) &
 result(info)
use qmckl
implicit none
integer(qmckl_context)  , intent(in)  :: context
character  , intent(in)  :: transa, transb
integer*8  , intent(in)  :: m, n
integer*8  , intent(in)  :: lda
real*8     , intent(in)  :: A(lda,*)
integer*8  , intent(in)  :: ldb
real*8     , intent(in)  :: B(ldb,*)
integer*8  , intent(in)  :: ldc
real*8     , intent(out) :: C(ldc,*)

integer*8 :: i,j
real*8    :: x, y, z
integer   :: transab

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
 info = QMCKL_INVALID_CONTEXT
 return
endif

if (m <= 0_8) then
 info = QMCKL_INVALID_ARG_4
 return
endif

if (n <= 0_8) then
 info = QMCKL_INVALID_ARG_5
 return
endif

if (transa == 'N' .or. transa == 'n') then
 transab = 0
else if (transa == 'T' .or. transa == 't') then
 transab = 1
else
 transab = -100
endif

if (transb == 'N' .or. transb == 'n') then
 continue
else if (transa == 'T' .or. transa == 't') then
 transab = transab + 2
else
 transab = -100
endif

if (transab < 0) then
 info = QMCKL_INVALID_ARG_1
 return
endif

! check for LDA
if (iand(transab,1) == 0 .and. LDA < 3) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,1) == 1 .and. LDA < m) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,2) == 0 .and. LDA < 3) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,2) == 2 .and. LDA < m) then
 info = QMCKL_INVALID_ARG_7
 return
endif

! check for LDB
if (iand(transab,1) == 0 .and. LDB < 3) then
 info = QMCKL_INVALID_ARG_9
 return
endif

if (iand(transab,1) == 1 .and. LDB < n) then
 info = QMCKL_INVALID_ARG_9
 return
endif

if (iand(transab,2) == 0 .and. LDB < 3) then
 info = QMCKL_INVALID_ARG_9
 return
endif

if (iand(transab,2) == 2 .and. LDB < n) then
 info = QMCKL_INVALID_ARG_9
 return
endif

! check for LDC
if (LDC < m) then
 info = QMCKL_INVALID_ARG_11
 return
endif


select case (transab)

case(0)

 do j=1,n
    do i=1,m
       x = A(1,i) - B(1,j)
       y = A(2,i) - B(2,j)
       z = A(3,i) - B(3,j)
       C(i,j) = x*x + y*y + z*z
    end do
    C(:,j) = dsqrt(C(:,j))
 end do

case(1)

 do j=1,n
    do i=1,m
       x = A(i,1) - B(1,j)
       y = A(i,2) - B(2,j)
       z = A(i,3) - B(3,j)
       C(i,j) = x*x + y*y + z*z
    end do
    C(:,j) = dsqrt(C(:,j))
 end do

case(2)

 do j=1,n
    do i=1,m
       x = A(1,i) - B(j,1)
       y = A(2,i) - B(j,2)
       z = A(3,i) - B(j,3)
       C(i,j) = x*x + y*y + z*z
    end do
    C(:,j) = dsqrt(C(:,j))
 end do

case(3)

 do j=1,n
    do i=1,m
       x = A(i,1) - B(j,1)
       y = A(i,2) - B(j,2)
       z = A(i,3) - B(j,3)
       C(i,j) = x*x + y*y + z*z
    end do
    C(:,j) = dsqrt(C(:,j))
 end do

end select

end function qmckl_distance_f

Performance

This function is more efficient when A and B are transposed.

Rescaled Distance

`qmckl_distance_rescaled`

qmckl_distance_rescaled computes the matrix of the rescaled distances between all pairs of points in two sets, one point within each set:

\[ C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa \]

If the input array is normal ('N'), the xyz coordinates are in the leading dimension: [n][3] in C and (3,n) in Fortran.

qmckl_context	context	in	Global state
char	transa	in	Array `A` is `'N'`: Normal, `'T'`: Transposed
char	transb	in	Array `B` is `'N'`: Normal, `'T'`: Transposed
int64_t	m	in	Number of points in the first set
int64_t	n	in	Number of points in the second set
double	A[][lda]	in	Array containing the $m \times 3$ matrix $A$
int64_t	lda	in	Leading dimension of array `A`
double	B[][ldb]	in	Array containing the $n \times 3$ matrix $B$
int64_t	ldb	in	Leading dimension of array `B`
double	C[n][ldc]	out	Array containing the $m \times n$ matrix $C$
int64_t	ldc	in	Leading dimension of array `C`
double	rescale_factor_kappa	in	Factor for calculating rescaled distances

Requirements

context is not QMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3 if transa == 'N'
lda >= m if transa == 'T'
ldb >= 3 if transb == 'N'
ldb >= n if transb == 'T'
ldc >= m
A is allocated with at least $3 \times m \times 8$ bytes
B is allocated with at least $3 \times n \times 8$ bytes
C is allocated with at least $m \times n \times 8$ bytes

C header

qmckl_exit_code qmckl_distance_rescaled (
      const qmckl_context context,
      const char transa,
      const char transb,
      const int64_t m,
      const int64_t n,
      const double* A,
      const int64_t lda,
      const double* B,
      const int64_t ldb,
      double* const C,
      const int64_t ldc,
      const double  rescale_factor_kappa);

Source

integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, &
 A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
 result(info)
use qmckl
implicit none
integer(qmckl_context)  , intent(in)  :: context
character  , intent(in)  :: transa, transb
integer*8  , intent(in)  :: m, n
integer*8  , intent(in)  :: lda
real*8     , intent(in)  :: A(lda,*)
integer*8  , intent(in)  :: ldb
real*8     , intent(in)  :: B(ldb,*)
integer*8  , intent(in)  :: ldc
real*8     , intent(out) :: C(ldc,*)
real*8     , intent(in)  :: rescale_factor_kappa

integer*8 :: i,j
real*8    :: x, y, z, dist, rescale_factor_kappa_inv
integer   :: transab

rescale_factor_kappa_inv = 1.0d0/rescale_factor_kappa;

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
 info = QMCKL_INVALID_CONTEXT
 return
endif

if (m <= 0_8) then
 info = QMCKL_INVALID_ARG_4
 return
endif

if (n <= 0_8) then
 info = QMCKL_INVALID_ARG_5
 return
endif

if (transa == 'N' .or. transa == 'n') then
 transab = 0
else if (transa == 'T' .or. transa == 't') then
 transab = 1
else
 transab = -100
endif

if (transb == 'N' .or. transb == 'n') then
 continue
else if (transa == 'T' .or. transa == 't') then
 transab = transab + 2
else
 transab = -100
endif

! check for LDA
if (transab < 0) then
 info = QMCKL_INVALID_ARG_1
 return
endif

if (iand(transab,1) == 0 .and. LDA < 3) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,1) == 1 .and. LDA < m) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,2) == 0 .and. LDA < 3) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,2) == 2 .and. LDA < m) then
 info = QMCKL_INVALID_ARG_7
 return
endif

! check for LDB
if (iand(transab,1) == 0 .and. LDB < 3) then
 info = QMCKL_INVALID_ARG_9
 return
endif

if (iand(transab,1) == 1 .and. LDB < n) then
 info = QMCKL_INVALID_ARG_9
 return
endif

if (iand(transab,2) == 0 .and. LDB < 3) then
 info = QMCKL_INVALID_ARG_9
 return
endif

if (iand(transab,2) == 2 .and. LDB < n) then
 info = QMCKL_INVALID_ARG_9
 return
endif

! check for LDC
if (LDC < m) then
 info = QMCKL_INVALID_ARG_11
 return
endif


select case (transab)

case(0)

 do j=1,n
    do i=1,m
       x = A(1,i) - B(1,j)
       y = A(2,i) - B(2,j)
       z = A(3,i) - B(3,j)
       dist = dsqrt(x*x + y*y + z*z)
       C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
    end do
 end do

case(1)

 do j=1,n
    do i=1,m
       x = A(i,1) - B(1,j)
       y = A(i,2) - B(2,j)
       z = A(i,3) - B(3,j)
       dist = dsqrt(x*x + y*y + z*z)
       C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
    end do
 end do

case(2)

 do j=1,n
    do i=1,m
       x = A(1,i) - B(j,1)
       y = A(2,i) - B(j,2)
       z = A(3,i) - B(j,3)
       dist = dsqrt(x*x + y*y + z*z)
       C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
    end do
 end do

case(3)

 do j=1,n
    do i=1,m
       x = A(i,1) - B(j,1)
       y = A(i,2) - B(j,2)
       z = A(i,3) - B(j,3)
       dist = dsqrt(x*x + y*y + z*z)
       C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
    end do
 end do

end select

end function qmckl_distance_rescaled_f

Performance

This function is more efficient when A and B are transposed.

Rescaled Distance Derivatives

`qmckl_distance_rescaled_deriv_e`

qmckl_distance_rescaled_deriv_e computes the matrix of the gradient and laplacian of the rescaled distance with respect to the electron coordinates. The derivative is a rank 3 tensor. The first dimension has a dimension of 4 of which the first three coordinates contains the gradient vector and the last index is the laplacian.

\[ C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa \]

Here the gradient is defined as follows:

\[ \nabla (C_{ij}(\mathbf{r}_{ee})) = \left(\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta x},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta y},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta z} \right) \] and the laplacian is defined as follows:

\[ \triangle (C_{ij}(r_{ee})) = \frac{\delta^2}{\delta x^2} + \frac{\delta^2}{\delta y^2} + \frac{\delta^2}{\delta z^2} \]

Using the above three formulae, the expression for the gradient and laplacian is as follows:

\[ \frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta x} = \frac{|(x_i - x_j)|}{r_{ij}} (1 - \kappa R_{ij}) \]

\[ \frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta y} = \frac{|(y_i - y_j)|}{r_{ij}} (1 - \kappa R_{ij}) \]

\[ \frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta z} = \frac{|(z_i - z_j)|}{r_{ij}} (1 - \kappa R_{ij}) \]

\[ \Delta(C_{ij}(r_{ee}) = \left[ \frac{2}{r_{ij}} - \kappa \right] (1-\kappa R_{ij}) \]

If the input array is normal ('N'), the xyz coordinates are in the leading dimension: [n][3] in C and (3,n) in Fortran.

qmckl_context	context	in	Global state
char	transa	in	Array `A` is `'N'`: Normal, `'T'`: Transposed
char	transb	in	Array `B` is `'N'`: Normal, `'T'`: Transposed
int64_t	m	in	Number of points in the first set
int64_t	n	in	Number of points in the second set
double	A[][lda]	in	Array containing the $m \times 3$ matrix $A$
int64_t	lda	in	Leading dimension of array `A`
double	B[][ldb]	in	Array containing the $n \times 3$ matrix $B$
int64_t	ldb	in	Leading dimension of array `B`
double	C[4][n][ldc]	out	Array containing the $4 \times m \times n$ matrix $C$
int64_t	ldc	in	Leading dimension of array `C`
double	rescale_factor_kappa	in	Factor for calculating rescaled distances derivatives

Requirements

context is not QMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3 if transa == 'N'
lda >= m if transa == 'T'
ldb >= 3 if transb == 'N'
ldb >= n if transb == 'T'
ldc >= m
A is allocated with at least $3 \times m \times 8$ bytes
B is allocated with at least $3 \times n \times 8$ bytes
C is allocated with at least $4 \times m \times n \times 8$ bytes

C header

qmckl_exit_code qmckl_distance_rescaled_deriv_e (
      const qmckl_context context,
      const char transa,
      const char transb,
      const int64_t m,
      const int64_t n,
      const double* A,
      const int64_t lda,
      const double* B,
      const int64_t ldb,
      double* const C,
      const int64_t ldc,
      const double  rescale_factor_kappa);

Source

integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n, &
 A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
 result(info)
use qmckl
implicit none
integer(qmckl_context)  , intent(in)  :: context
character  , intent(in)  :: transa, transb
integer*8  , intent(in)  :: m, n
integer*8  , intent(in)  :: lda
real*8     , intent(in)  :: A(lda,*)
integer*8  , intent(in)  :: ldb
real*8     , intent(in)  :: B(ldb,*)
integer*8  , intent(in)  :: ldc
real*8     , intent(out) :: C(4,ldc,*)
real*8     , intent(in)  :: rescale_factor_kappa

integer*8 :: i,j
real*8    :: x, y, z, dist, dist_inv
real*8    :: rescale_factor_kappa_inv, rij
integer   :: transab

rescale_factor_kappa_inv = 1.0d0/rescale_factor_kappa;

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
 info = QMCKL_INVALID_CONTEXT
 return
endif

if (m <= 0_8) then
 info = QMCKL_INVALID_ARG_4
 return
endif

if (n <= 0_8) then
 info = QMCKL_INVALID_ARG_5
 return
endif

if (transa == 'N' .or. transa == 'n') then
 transab = 0
else if (transa == 'T' .or. transa == 't') then
 transab = 1
else
 transab = -100
endif

if (transb == 'N' .or. transb == 'n') then
 continue
else if (transa == 'T' .or. transa == 't') then
 transab = transab + 2
else
 transab = -100
endif

! check for LDA
if (transab < 0) then
 info = QMCKL_INVALID_ARG_1
 return
endif

if (iand(transab,1) == 0 .and. LDA < 3) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,1) == 1 .and. LDA < m) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,2) == 0 .and. LDA < 3) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (iand(transab,2) == 2 .and. LDA < m) then
 info = QMCKL_INVALID_ARG_7
 return
endif

! check for LDB
if (iand(transab,1) == 0 .and. LDB < 3) then
 info = QMCKL_INVALID_ARG_9
 return
endif

if (iand(transab,1) == 1 .and. LDB < n) then
 info = QMCKL_INVALID_ARG_9
 return
endif

if (iand(transab,2) == 0 .and. LDB < 3) then
 info = QMCKL_INVALID_ARG_9
 return
endif

if (iand(transab,2) == 2 .and. LDB < n) then
 info = QMCKL_INVALID_ARG_9
 return
endif

! check for LDC
if (LDC < m) then
 info = QMCKL_INVALID_ARG_11
 return
endif


select case (transab)

case(0)

 do j=1,n
    do i=1,m
       x = A(1,i) - B(1,j)
       y = A(2,i) - B(2,j)
       z = A(3,i) - B(3,j)
       dist = dsqrt(x*x + y*y + z*z)
       dist_inv = 1.0d0/dist
       rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
       C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
    end do
 end do

case(1)

 do j=1,n
    do i=1,m
       x = A(i,1) - B(1,j)
       y = A(i,2) - B(2,j)
       z = A(i,3) - B(3,j)
       dist = dsqrt(x*x + y*y + z*z)
       dist_inv = 1.0d0/dist
       rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
       C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
    end do
 end do

case(2)

 do j=1,n
    do i=1,m
       x = A(1,i) - B(j,1)
       y = A(2,i) - B(j,2)
       z = A(3,i) - B(j,3)
       dist = dsqrt(x*x + y*y + z*z)
       dist_inv = 1.0d0/dist
       rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
       C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
    end do
 end do

case(3)

 do j=1,n
    do i=1,m
       x = A(i,1) - B(j,1)
       y = A(i,2) - B(j,2)
       z = A(i,3) - B(j,3)
       dist = dsqrt(x*x + y*y + z*z)
       dist_inv = 1.0d0/dist
       rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
       C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
       C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
    end do
 end do

end select

end function qmckl_distance_rescaled_deriv_e_f

Performance

This function is more efficient when A and B are transposed.

46 KiB Raw Blame History

Inter-particle distances

Squared distance

qmckl_distance_sq

Requirements

C header

Source

Performance

Distance

qmckl_distance

Requirements

C header

Source

Performance

Rescaled Distance

qmckl_distance_rescaled

Requirements

C header

Source

Performance

Rescaled Distance Derivatives

qmckl_distance_rescaled_deriv_e

Requirements

C header

Source

Performance

46 KiB

Raw Blame History

`qmckl_distance_sq`

`qmckl_distance`

`qmckl_distance_rescaled`

`qmckl_distance_rescaled_deriv_e`