#+TITLE: Distances
#+SETUPFILE: ../docs/theme.setup

Functions for the computation of distances between particles.

* Headers                                                         :noexport:

  #+NAME: filename
  #+begin_src elisp tangle: no 
(file-name-nondirectory (substring buffer-file-name 0 -4))
  #+end_src

  #+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "munit.h"
MunitResult test_<<filename()>>() {
  qmckl_context context;
  context = qmckl_context_create();

  #+end_src

* 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.

  #+begin_src c :comments org :tangle (eval h)
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);
    #+end_src

*** Source
    #+begin_src f90 :tangle (eval f)
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
    #+end_src
    
*** C interface                                                    :noexport:
    #+begin_src f90 :tangle (eval f)
integer(c_int32_t) function qmckl_distance_sq(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC) &
     bind(C) result(info)
  use, intrinsic :: iso_c_binding
  implicit none
  integer (c_int64_t) , intent(in) , value :: context
  character (c_char)  , intent(in) , value :: transa, transb
  integer (c_int64_t) , intent(in) , value :: m, n
  integer (c_int64_t) , intent(in) , value :: lda
  real    (c_double)  , intent(in)         :: A(lda,3)
  integer (c_int64_t) , intent(in) , value :: ldb
  real    (c_double)  , intent(in)         :: B(ldb,3)
  integer (c_int64_t) , intent(in) , value :: ldc
  real    (c_double)  , intent(out)        :: C(ldc,n)

  integer, external :: qmckl_distance_sq_f
  info = qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC)
end function qmckl_distance_sq
    #+end_src

    #+begin_src f90 :tangle (eval fh)
  interface
     integer(c_int32_t) function qmckl_distance_sq(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC) &
          bind(C) 
       use, intrinsic :: iso_c_binding
       implicit none
       integer (c_int64_t) , intent(in) , value :: context
       character (c_char)  , intent(in) , value :: transa, transb
       integer (c_int64_t) , intent(in) , value :: m, n
       integer (c_int64_t) , intent(in) , value :: lda
       integer (c_int64_t) , intent(in) , value :: ldb
       integer (c_int64_t) , intent(in) , value :: ldc
       real    (c_double)  , intent(in)         :: A(lda,3)
       real    (c_double)  , intent(in)         :: B(ldb,3)
       real    (c_double)  , intent(out)        :: C(ldc,n)
     end function qmckl_distance_sq
  end interface
    #+end_src

*** Test                                                           :noexport:
    #+begin_src f90 :tangle (eval f_test)
integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C)
  use qmckl
  implicit none
  integer(c_int64_t), intent(in), value :: context

  double precision, allocatable :: A(:,:), B(:,:), C(:,:)
  integer*8                     :: m, n, LDA, LDB, LDC
  double precision              :: x
  integer*8                     :: i,j 

  m = 5
  n = 6
  LDA = m
  LDB = n
  LDC = 5

  allocate( A(LDA,m), B(LDB,n), C(LDC,n) )

  do j=1,m
     do i=1,m
        A(i,j) = -10.d0 + dble(i+j)
     end do
  end do
  do j=1,n
     do i=1,n
        B(i,j) = -1.d0 + dble(i*j)
     end do
  end do

  test_qmckl_distance_sq = &
       qmckl_distance_sq(context, 'X', 't', m, n, A, LDA, B, LDB, C, LDC)

  if (test_qmckl_distance_sq == 0) return 

  test_qmckl_distance_sq = &
       qmckl_distance_sq(context, 't', 'X', m, n, A, LDA, B, LDB, C, LDC)

  if (test_qmckl_distance_sq == 0) return 

  test_qmckl_distance_sq = &
       qmckl_distance_sq(context, 'T', 't', m, n, A, LDA, B, LDB, C, LDC)

  if (test_qmckl_distance_sq /= 0) return 

  test_qmckl_distance_sq = -1

  do j=1,n
     do i=1,m
        x =  (A(i,1)-B(j,1))**2 + &
             (A(i,2)-B(j,2))**2 + &
             (A(i,3)-B(j,3))**2
        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
     end do
  end do
  
  test_qmckl_distance_sq = &
       qmckl_distance_sq(context, 'n', 'T', m, n, A, LDA, B, LDB, C, LDC)

  if (test_qmckl_distance_sq /= 0) return 

  test_qmckl_distance_sq = -1

  do j=1,n
     do i=1,m
        x =  (A(1,i)-B(j,1))**2 + &
             (A(2,i)-B(j,2))**2 + &
             (A(3,i)-B(j,3))**2
        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
     end do
  end do

  test_qmckl_distance_sq = &
       qmckl_distance_sq(context, 'T', 'n', m, n, A, LDA, B, LDB, C, LDC)

  if (test_qmckl_distance_sq /= 0) return 

  test_qmckl_distance_sq = -1

  do j=1,n
     do i=1,m
        x =  (A(i,1)-B(1,j))**2 + &
             (A(i,2)-B(2,j))**2 + &
             (A(i,3)-B(3,j))**2
        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
     end do
  end do

  test_qmckl_distance_sq = &
       qmckl_distance_sq(context, 'n', 'N', m, n, A, LDA, B, LDB, C, LDC)

  if (test_qmckl_distance_sq /= 0) return 

  test_qmckl_distance_sq = -1

  do j=1,n
     do i=1,m
        x =  (A(1,i)-B(1,j))**2 + &
             (A(2,i)-B(2,j))**2 + &
             (A(3,i)-B(3,j))**2
        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
     end do
  end do
  
  test_qmckl_distance_sq = 0
  
  deallocate(A,B,C)
end function test_qmckl_distance_sq
    #+end_src
    
    #+begin_src c :comments link :tangle (eval c_test)
int test_qmckl_distance_sq(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_distance_sq(context));
    #+end_src
* End of files                                                     :noexport:

   #+begin_src c :comments link :tangle (eval c_test)
  if (qmckl_context_destroy(context) != QMCKL_SUCCESS)
    return QMCKL_FAILURE;
  return MUNIT_OK;
}

   #+end_src
  

# -*- mode: org -*-
# vim: syntax=c