# -*- mode: org -*-
# vim: syntax=c
#+TITLE: Computation of distances

#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/readtheorg.css"/>
#+HTML_HEAD: <script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.3/jquery.min.js"></script>
#+HTML_HEAD: <script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.4/js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/lib/js/jquery.stickytableheaders.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/readtheorg/js/readtheorg.js"></script>

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                                                         :noexport:
    #+BEGIN_SRC C :comments link  :tangle qmckl_distance.h
#ifndef QMCKL_DISTANCE_H
#define QMCKL_DISTANCE_H
#include "qmckl_context.h"
    #+END_SRC

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

*** Test                                                           :noexport:
    #+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;
  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_{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           |
    | =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= >= 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 function qmckl_distance_sq_f(context, m, n, A, LDA, B, LDB, C, LDC) result(info)
  implicit none
  integer*8  , intent(in)  :: context
  integer*8  , intent(in)  :: m, n
  integer*8  , intent(in)  :: lda
  real*8     , intent(in)  :: A(lda,3)
  integer*8  , intent(in)  :: ldb
  real*8     , intent(in)  :: B(ldb,3)
  integer*8  , intent(in)  :: ldc
  real*8     , intent(out) :: C(ldc,n)

  integer*8 :: i,j
  real*8    :: 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_f
    #+END_SRC

*** C interface                                                   :noexport:
    #+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, external :: qmckl_distance_sq_f
  info = qmckl_distance_sq_f(context, m, n, A, LDA, B, LDB, C, LDC)
end function qmckl_distance_sq
    #+END_SRC

*** Test                                                           :noexport:
  #+BEGIN_SRC f90 :comments link :tangle test_qmckl_distance_f.f90
  integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C)
    use iso_c_binding
    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 

    integer, external             :: qmckl_distance_sq_f

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

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

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

    test_qmckl_distance_sq = qmckl_distance_sq_f(context, 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-12 ) 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 test_qmckl_distance.c
  int test_qmckl_distance_sq(qmckl_context context);
  munit_assert_int(0, ==, test_qmckl_distance_sq(context));
  #+END_SRC
* End of files

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

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

  #+END_SRC