# -*- mode: org -*-
# vim: syntax=c
#+TITLE: 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 test   file : =test_qmckl_distance.c= 

*** Header 
    #+BEGIN_SRC C :comments link  :tangle qmckl_distance.h
#ifndef QMCKL_DISTANCE_H
#define QMCKL_DISTANCE_H
#include "qmckl_context.h"
    #+END_SRC

*** Source
    #+BEGIN_SRC f90 :comments link :tangle qmckl_distance.f90
    #+END_SRC

*** Test
    #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
#include <math.h>
#include "qmckl.h"
#include "munit.h"
MunitResult test_qmckl_distance() {
  qmckl_context context;
  int64_t  m, n, LDA, LDB, LDC;
  double *A, *B, *C ;
  int i, j;

  context = qmckl_context_create();

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

  A = (double*) qmckl_malloc (context, LDA*4*sizeof(double));
  B = (double*) qmckl_malloc (context, LDB*3*sizeof(double));
  C = (double*) qmckl_malloc (context, LDC*n*sizeof(double));

  for (j=0 ; j<3 ; j++) {
    for (i=0 ; i<m ; i++) {
      A[i+j*LDA] = -10. + (double) (i+j);
    }
  }

  for (j=0 ; j<3 ; j++) {
    for (i=0 ; i<n ; i++) {
      B[i+j*LDB] = -1. + (double) (i*j);
    }
  }

    #+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^2} = \sum_{k=1}^3 (A_{i,k}-B_{j,k})^2
   \]
   
*** Arguments

    | =context= | input  | Global state                            |
    | =m=       | input  | Number of points in the first set       |
    | =n=       | input  | Number of points in the second set      |
    | =LDA=     | input  | Leading dimension of array =A=            |
    | =A=       | input  | Array containing the $3 \times m$ matrix $A$ |
    | =LDB=     | input  | Leading dimension of array =B=            |
    | =B=       | input  | Array containing the $3 \times n$ matrix $B$ |
    | =LDC=     | input  | Leading dimension of array =C=            |
    | =C=       | output | Array containing the $m \times n$ matrix $C$ |
    | =info=    | output | exit status is zero upon success        |

*** Requirements

    - =context= is not 0
    - =m= > 0
    - =n= > 0
    - =LDA= >= m
    - =LDB= >= n
    - =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

*** Header
    #+BEGIN_SRC C :comments link :tangle qmckl_distance.h
qmckl_exit_code qmckl_distance_sq(qmckl_context context, 
				  int64_t m, int64_t n,
				  double *A, int64_t LDA,
				  double *B, int64_t LDB,
				  double *C, int64_t LDC);
    #+END_SRC

*** Source
    #+BEGIN_SRC f90 :comments link  :tangle qmckl_distance.f90
integer(c_int32_t) function qmckl_distance_sq(context, 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
  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 (c_int64_t) :: i,j
  real    (c_double)  :: x, y, z
  
  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 (LDA < m) then
     info = -4
     return
  endif
  
  if (LDB < n) then
     info = -5
     return
  endif
  
  if (LDC < m) then
     info = -6
     return
  endif
  
  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 function qmckl_distance_sq
    #+END_SRC

*** Test
  #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c

  munit_assert_int64(QMCKL_SUCCESS, ==, 
		     qmckl_distance_sq(context, m, n, A, LDA, B, LDB, C, LDC) );

  for (j=0 ; j<n ; j++) {
    for (i=0 ; i<m ; i++) {
      munit_assert_double_equal(C[i+j*LDC], 
				pow(A[i      ]-B[j      ],2) +
				pow(A[i+  LDA]-B[j+  LDB],2) +
				pow(A[i+2*LDA]-B[j+2*LDB],2) ,
				14 );
    }
  } 

  #+END_SRC
* End of files

*** Header
  #+BEGIN_SRC C :comments link :tangle qmckl_distance.h
#endif
  #+END_SRC

*** Test
  #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
  qmckl_free(A);
  qmckl_free(B);
  qmckl_free(C);
  if (qmckl_context_destroy(context) != QMCKL_SUCCESS) 
    return QMCKL_FAILURE;
  return MUNIT_OK;
} 

  #+END_SRC