31 KiB
Inter-particle distances
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 notQMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3
iftransa == 'N'
lda >= m
iftransa == 'T'
ldb >= 3
iftransb == 'N'
ldb >= n
iftransb == 'T'
ldc >= m
A
is allocated with at least $3 \times m \times 8$ bytesB
is allocated with at least $3 \times n \times 8$ bytesC
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 notQMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3
iftransa == 'N'
lda >= m
iftransa == 'T'
ldb >= 3
iftransb == 'N'
ldb >= n
iftransb == 'T'
ldc >= m
A
is allocated with at least $3 \times m \times 8$ bytesB
is allocated with at least $3 \times n \times 8$ bytesC
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 notQMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3
iftransa == 'N'
lda >= m
iftransa == 'T'
ldb >= 3
iftransb == 'N'
ldb >= n
iftransb == 'T'
ldc >= m
A
is allocated with at least $3 \times m \times 8$ bytesB
is allocated with at least $3 \times n \times 8$ bytesC
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.