1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-09-16 17:55:29 +02:00
qmckl/src/qmckl_distance.org
2021-03-19 00:10:35 +01:00

9.8 KiB

Distances

Functions for the computation of distances between particles.

Squared distance

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 \]

context input Global state
transa input Array A is N: Normal, T: Transposed
transb input Array B is N: Normal, T: Transposed
m input Number of points in the first set
n input Number of points in the second set
A(lda,3) input Array containing the $m \times 3$ matrix $A$
lda input Leading dimension of array A
B(ldb,3) input Array containing the $n \times 3$ matrix $B$
ldb input Leading dimension of array B
C(ldc,n) output Array containing the $m \times n$ matrix $C$
ldc input Leading dimension of array C

Requirements

  • context is not 0
  • m > 0
  • n > 0
  • lda >= 3 if transa is N
  • lda >= m if transa is T
  • ldb >= 3 if transb is N
  • ldb >= n if transb is 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

Performance

This function might be more efficient when A and B are transposed.

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,
                                const double *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*8  , 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 = 0

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