Computation of distances

1. Squared distance
- 1.1. qmckl_distance_sq

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

1 Squared distance

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

1.1.1 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`

1.1.2 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 if transa is =
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.3 Performance

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

1.1.4 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);

1.1.5 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