mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-08-17 19:01:43 +02:00
10 KiB
10 KiB
Computation of distances
Function for the computation of distances between particles.
3 files are produced:
- a header file :
qmckl_distance.h
- a source file :
qmckl_distance.f90
- a C test file :
test_qmckl_distance.c
- a Fortran test file :
test_qmckl_distance_f.f90
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 \]
Arguments
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 0m
> 0n
> 0lda
>= 3 iftransa
isN
lda
>= m iftransa
isT
ldb
>= 3 iftransb
isN
ldb
>= n iftransb
isT
ldc
>= m iftransa
is =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
Performance
This function might be more efficient when A
and B
are
transposed.
Header
qmckl_exit_code qmckl_distance_sq(qmckl_context context,
char transa, char transb,
int64_t m, int64_t n,
double *A, int64_t lda,
double *B, int64_t ldb,
double *C, int64_t ldc);
Source
integer function qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC) result(info)
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 == 0_8) then
info = -1
return
endif
if (m <= 0_8) then
info = -2
return
endif
if (n <= 0_8) then
info = -3
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 = -4
return
endif
if (iand(transab,1) == 0 .and. LDA < 3) then
info = -5
return
endif
if (iand(transab,1) == 1 .and. LDA < m) then
info = -6
return
endif
if (iand(transab,2) == 0 .and. LDA < 3) then
info = -6
return
endif
if (iand(transab,2) == 2 .and. LDA < m) then
info = -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