Inter-particle distances
+Table of Contents
++Functions for the computation of distances between particles. +
+ +1 Squared distance
+1.1 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 + \] +
+ +qmcklcontext | +context | +in | +Global state | +
char | +transa | +in | +Array A is 'N' : Normal, 'T' : Transposed |
+
char | +transb | +in | +Array B is 'N' : Normal, 'T' : Transposed |
+
int64t | +m | +in | +Number of points in the first set | +
int64t | +n | +in | +Number of points in the second set | +
double | +A[3][lda] | +in | +Array containing the \(m \times 3\) matrix \(A\) | +
int64t | +lda | +in | +Leading dimension of array A |
+
double | +B[3][ldb] | +in | +Array containing the \(n \times 3\) matrix \(B\) | +
int64t | +ldb | +in | +Leading dimension of array B |
+
double | +C[n][ldc] | +out | +Array containing the \(m \times n\) matrix \(C\) | +
int64t | +ldc | +in | +Leading dimension of array C |
+
1.1.1 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\) 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
+
1.1.2 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* C, + const int64_t ldc ); ++
1.1.3 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 ++
1.1.4 Performance
+
+This function might be more efficient when A
and B
are
+transposed.
+