qmckl/org/qmckl_distance.org

#+TITLE: Inter-particle distances
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org

Functions for the computation of distances between particles.

* Headers                                                         :noexport:
  #+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src

  #+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#include <stdio.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
int main() {
  qmckl_context context;
  context = qmckl_context_create();


  #+end_src

* Squared distance

** ~qmckl_distance_sq~
   :PROPERTIES:
   :Name:     qmckl_distance_sq
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

 ~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
   \]

   #+NAME: qmckl_distance_sq_args
   | Variable  | Type             | In/Out | Description                                   |
   |-----------+------------------+--------+-----------------------------------------------|
   | ~context~ | ~qmckl_context~  | in     | Global state                                  |
   | ~transa~  | ~char~           | in     | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
   | ~transb~  | ~char~           | in     | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
   | ~m~       | ~int64_t~        | in     | Number of points in the first set             |
   | ~n~       | ~int64_t~        | in     | Number of points in the second set            |
   | ~A~       | ~double[][lda]~  | in     | Array containing the $m \times 3$ matrix $A$  |
   | ~lda~     | ~int64_t~        | in     | Leading dimension of array ~A~                |
   | ~B~       | ~double[][ldb]~  | in     | Array containing the $n \times 3$ matrix $B$  |
   | ~ldb~     | ~int64_t~        | in     | Leading dimension of array ~B~                |
   | ~C~       | ~double[n][ldc]~ | out    | Array containing the $m \times n$ matrix $C$  |
   | ~ldc~     | ~int64_t~        | in     | Leading dimension of array ~C~                |

   Requirements:

    - ~context~ is not ~QMCKL_NULL_CONTEXT~
    - ~m > 0~
    - ~n > 0~
    - ~lda >= 3~ if ~transa == 'N'~
    - ~lda >= m~ if ~transa == 'T'~
    - ~ldb >= 3~ if ~transb == 'N'~
    - ~ldb >= n~ if ~transb == '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

    #+CALL: generate_c_header(table=qmckl_distance_sq_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    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,
          double* const C,
          const int64_t ldc ); 
    #+end_src

    #+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(qmckl_context)  , 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 = QMCKL_SUCCESS

  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 (transb == 'T' .or. transb == '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

*** Performance

    This function is more efficient when ~A~ and ~B~ are
    transposed.

   #+CALL: generate_c_interface(table=qmckl_distance_sq_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval f) :comments org :exports none
   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           , intent(in)  , value :: transa
     character           , intent(in)  , value :: transb
     integer (c_int64_t) , intent(in)  , value :: m
     integer (c_int64_t) , intent(in)  , value :: n
     real    (c_double ) , intent(in)          :: A(lda,*)
     integer (c_int64_t) , intent(in)  , value :: lda
     real    (c_double ) , intent(in)          :: B(ldb,*)
     integer (c_int64_t) , intent(in)  , value :: ldb
     real    (c_double ) , intent(out)         :: C(ldc,n)
     integer (c_int64_t) , intent(in)  , value :: ldc

     integer(c_int32_t), 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

   #+CALL: generate_f_interface(table=qmckl_distance_sq_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   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
       import
       implicit none

       integer (c_int64_t) , intent(in)  , value :: context
       character           , intent(in)  , value :: transa
       character           , intent(in)  , value :: transb
       integer (c_int64_t) , intent(in)  , value :: m
       integer (c_int64_t) , intent(in)  , value :: n
       real    (c_double ) , intent(in)          :: A(lda,*)
       integer (c_int64_t) , intent(in)  , value :: lda
       real    (c_double ) , intent(in)          :: B(ldb,*)
       integer (c_int64_t) , intent(in)  , value :: ldb
       real    (c_double ) , intent(out)         :: C(ldc,n)
       integer (c_int64_t) , intent(in)  , value :: ldc

     end function qmckl_distance_sq
   end interface
   #+end_src

*** Test                                                           :noexport:
    #+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_distance_sq(context) bind(C)

  use qmckl
  use qmckl_verificarlo_f
  use iso_c_binding

  implicit none

  integer(qmckl_context), intent(in), value :: context
  logical(C_BOOL) :: vfc_err

  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)

  vfc_err = qmckl_probe("distance"//C_NULL_CHAR, "distance_sq_Xt_2_2"//C_NULL_CHAR, DBLE(C(2,2)))

  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)

  vfc_err = qmckl_probe("distance"//C_NULL_CHAR, "distance_sq_tX_2_2"//C_NULL_CHAR, DBLE(C(2,2)))

  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)

  vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_sq_Tt_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))

  if (test_qmckl_distance_sq == 0) return


  test_qmckl_distance_sq = QMCKL_FAILURE

  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
#ifndef VFC_CI
        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14) return
#endif
     end do
  end do

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

  vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_sq_nT_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))


  test_qmckl_distance_sq = QMCKL_FAILURE

  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
#ifndef VFC_CI
        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14) return
#endif
     end do
  end do

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

  vfc_err =  qmckl_probe_check("distance"//C_NULL_CHAR, "distance_sq_Tn_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))

  if (test_qmckl_distance_sq == 0) return

  test_qmckl_distance_sq = QMCKL_FAILURE

  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
#ifndef VFC_CI
        if ( dabs(1.d0 - C(i,j)/x) > 1.d-14) return
#endif
     end do
  end do

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

  vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_sq_nN_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))

  test_qmckl_distance_sq = QMCKL_FAILURE

  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 = QMCKL_SUCCESS

  deallocate(A,B,C)
end function test_qmckl_distance_sq
    #+end_src

    #+begin_src c :comments link :tangle (eval c_test)
    qmckl_exit_code test_qmckl_distance_sq(qmckl_context context);
    assert(test_qmckl_distance_sq(context) == QMCKL_SUCCESS);
    #+end_src
* Distance

** ~qmckl_distance~
   :PROPERTIES:
   :Name:     qmckl_distance
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

 ~qmckl_distance~ computes the matrix of the distances between all
   pairs of points in two sets, one point within each set:

   \[
   C_{ij} = \sqrt{\sum_{k=1}^3 (A_{k,i}-B_{k,j})^2}
   \]

   If the input array is normal (~'N'~), the xyz coordinates are in
   the leading dimension: ~[n][3]~ in C and ~(3,n)~ in Fortran.

   #+NAME: qmckl_distance_args
   | Variable  | Type             | In/Out | Description                                   |
   |-----------+------------------+--------+-----------------------------------------------|
   | ~context~ | ~qmckl_context~  | in     | Global state                                  |
   | ~transa~  | ~char~           | in     | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
   | ~transb~  | ~char~           | in     | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
   | ~m~       | ~int64_t~        | in     | Number of points in the first set             |
   | ~n~       | ~int64_t~        | in     | Number of points in the second set            |
   | ~A~       | ~double[][lda]~  | in     | Array containing the $m \times 3$ matrix $A$  |
   | ~lda~     | ~int64_t~        | in     | Leading dimension of array ~A~                |
   | ~B~       | ~double[][ldb]~  | in     | Array containing the $n \times 3$ matrix $B$  |
   | ~ldb~     | ~int64_t~        | in     | Leading dimension of array ~B~                |
   | ~C~       | ~double[n][ldc]~ | out    | Array containing the $m \times n$ matrix $C$  |
   | ~ldc~     | ~int64_t~        | in     | Leading dimension of array ~C~                |

*** Requirements

    - ~context~ is not ~QMCKL_NULL_CONTEXT~
    - ~m > 0~
    - ~n > 0~
    - ~lda >= 3~ if ~transa == 'N'~
    - ~lda >= m~ if ~transa == 'T'~
    - ~ldb >= 3~ if ~transb == 'N'~
    - ~ldb >= n~ if ~transb == '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

*** C header

    #+CALL: generate_c_header(table=qmckl_distance_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_distance (
          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,
          double* const C,
          const int64_t ldc ); 
    #+end_src

*** Source
    #+begin_src f90 :tangle (eval f)
integer function qmckl_distance_f(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , 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 = QMCKL_SUCCESS

  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 (transb == 'T' .or. transb == 't') then
     transab = transab + 2
  else
     transab = -100
  endif

  if (transab < 0) then
     info = QMCKL_INVALID_ARG_1
     return
  endif

  ! check for LDA
  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

  ! check for LDB
  if (iand(transab,1) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,1) == 1 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,2) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,2) == 2 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  ! check for LDC
  if (LDC < m) then
     info = QMCKL_INVALID_ARG_11
     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
        C(:,j) = dsqrt(C(:,j))
     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
        C(:,j) = dsqrt(C(:,j))
     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
        C(:,j) = dsqrt(C(:,j))
     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
        C(:,j) = dsqrt(C(:,j))
     end do

  end select

end function qmckl_distance_f
    #+end_src

*** Performance

    This function is more efficient when ~A~ and ~B~ are transposed.

** C interface                                                     :noexport:

   #+CALL: generate_c_interface(table=qmckl_distance_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval f) :comments org :exports none
   integer(c_int32_t) function qmckl_distance &
       (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           , intent(in)  , value :: transa
     character           , intent(in)  , value :: transb
     integer (c_int64_t) , intent(in)  , value :: m
     integer (c_int64_t) , intent(in)  , value :: n
     real    (c_double ) , intent(in)          :: A(lda,*)
     integer (c_int64_t) , intent(in)  , value :: lda
     real    (c_double ) , intent(in)          :: B(ldb,*)
     integer (c_int64_t) , intent(in)  , value :: ldb
     real    (c_double ) , intent(out)         :: C(ldc,n)
     integer (c_int64_t) , intent(in)  , value :: ldc

     integer(c_int32_t), external :: qmckl_distance_f
     info = qmckl_distance_f &
            (context, transa, transb, m, n, A, lda, B, ldb, C, ldc)

   end function qmckl_distance
   #+end_src

   #+CALL: generate_f_interface(table=qmckl_distance_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
     integer(c_int32_t) function qmckl_distance &
         (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) &
         bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

       integer (c_int64_t) , intent(in)  , value :: context
       character           , intent(in)  , value :: transa
       character           , intent(in)  , value :: transb
       integer (c_int64_t) , intent(in)  , value :: m
       integer (c_int64_t) , intent(in)  , value :: n
       real    (c_double ) , intent(in)          :: A(lda,*)
       integer (c_int64_t) , intent(in)  , value :: lda
       real    (c_double ) , intent(in)          :: B(ldb,*)
       integer (c_int64_t) , intent(in)  , value :: ldb
       real    (c_double ) , intent(out)         :: C(ldc,n)
       integer (c_int64_t) , intent(in)  , value :: ldc

     end function qmckl_distance
   end interface
   #+end_src

*** Test                                                           :noexport:
    #+begin_src f90 :tangle (eval f_test)

integer(qmckl_exit_code) function test_qmckl_dist(context) bind(C)

  use qmckl
  use qmckl_verificarlo_f
  use iso_c_binding

  implicit none

  integer(qmckl_context), intent(in), value :: context
  logical(C_BOOL) :: vfc_err

  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_dist = &
       qmckl_distance(context, 'X', 't', m, n, A, LDA, B, LDB, C, LDC)

  vfc_err = qmckl_probe("distance"//C_NULL_CHAR, "distance_Xt_2_2"//C_NULL_CHAR, DBLE(C(2,2)))

  if (test_qmckl_dist == 0) return

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

  vfc_err = qmckl_probe("distance"//C_NULL_CHAR, "distance_tX_2_2"//C_NULL_CHAR, DBLE(C(2,2)))

  if (test_qmckl_dist == 0) return

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

  vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_Tt_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))

  if (test_qmckl_dist == 0) return


  test_qmckl_dist = QMCKL_FAILURE

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

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

  vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_nT_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))

  if (test_qmckl_dist == 0) return


  test_qmckl_dist = QMCKL_FAILURE

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

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

  vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_Tn_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))

  if (test_qmckl_dist == 0) return


  test_qmckl_dist = QMCKL_FAILURE

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

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

  vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_nN_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))

  if (test_qmckl_dist == 0) return


  test_qmckl_dist = QMCKL_FAILURE

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

  test_qmckl_dist = QMCKL_SUCCESS

  deallocate(A,B,C)
end function test_qmckl_dist
    #+end_src

    #+begin_src c :comments link :tangle (eval c_test)
    qmckl_exit_code test_qmckl_dist(qmckl_context context);
    assert(test_qmckl_dist(context) == QMCKL_SUCCESS);
    #+end_src

* Rescaled Distance

** ~qmckl_distance_rescaled~
   :PROPERTIES:
   :Name:     qmckl_distance_rescaled
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

 ~qmckl_distance_rescaled~ computes the matrix of the rescaled distances between all
   pairs of points in two sets, one point within each set:

   \[
   C_{ij} = \left( 1 - \exp \left(-\kappa C_{ij} \right) \right)/\kappa
   \]

   If the input array is normal (~'N'~), the xyz coordinates are in
   the leading dimension: ~[n][3]~ in C and ~(3,n)~ in Fortran.

   #+NAME: qmckl_distance_rescaled_args
   | Variable               | Type             | In/Out | Description                                   |
   |------------------------+------------------+--------+-----------------------------------------------|
   | ~context~              | ~qmckl_context~  | in     | Global state                                  |
   | ~transa~               | ~char~           | in     | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
   | ~transb~               | ~char~           | in     | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
   | ~m~                    | ~int64_t~        | in     | Number of points in the first set             |
   | ~n~                    | ~int64_t~        | in     | Number of points in the second set            |
   | ~A~                    | ~double[][lda]~  | in     | Array containing the $m \times 3$ matrix $A$  |
   | ~lda~                  | ~int64_t~        | in     | Leading dimension of array ~A~                |
   | ~B~                    | ~double[][ldb]~  | in     | Array containing the $n \times 3$ matrix $B$  |
   | ~ldb~                  | ~int64_t~        | in     | Leading dimension of array ~B~                |
   | ~C~                    | ~double[n][ldc]~ | out    | Array containing the $m \times n$ matrix $C$  |
   | ~ldc~                  | ~int64_t~        | in     | Leading dimension of array ~C~                |
   | ~rescale_factor_kappa~ | ~double~         | in     | Factor for calculating rescaled distances     |

*** Requirements

    - ~context~ is not ~QMCKL_NULL_CONTEXT~
    - ~m > 0~
    - ~n > 0~
    - ~lda >= 3~ if ~transa == 'N'~
    - ~lda >= m~ if ~transa == 'T'~
    - ~ldb >= 3~ if ~transb == 'N'~
    - ~ldb >= n~ if ~transb == '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

*** C header

    #+CALL: generate_c_header(table=qmckl_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_distance_rescaled (
          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,
          double* const C,
          const int64_t ldc,
          const double rescale_factor_kappa ); 
    #+end_src

*** Source
    #+begin_src f90 :tangle (eval f)
integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , 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,*)
  real*8     , intent(in)  :: rescale_factor_kappa

  integer*8 :: i,j
  real*8    :: x, y, z, dist, rescale_factor_kappa_inv
  integer   :: transab

  rescale_factor_kappa_inv = 1.0d0/rescale_factor_kappa;

  info = QMCKL_SUCCESS

  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 (transb == 'T' .or. transb == 't') then
     transab = transab + 2
  else
     transab = -100
  endif

  ! check for LDA
  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

  ! check for LDB
  if (iand(transab,1) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,1) == 1 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,2) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,2) == 2 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  ! check for LDC
  if (LDC < m) then
     info = QMCKL_INVALID_ARG_11
     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)
           dist = dsqrt(x*x + y*y + z*z)
           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
        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)
           dist = dsqrt(x*x + y*y + z*z)
           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
        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)
           dist = dsqrt(x*x + y*y + z*z)
           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
        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)
           dist = dsqrt(x*x + y*y + z*z)
           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
        end do
     end do

  end select

end function qmckl_distance_rescaled_f
    #+end_src

*** Performance

    This function is more efficient when ~A~ and ~B~ are transposed.

** C interface                                                     :noexport:

   #+CALL: generate_c_interface(table=qmckl_distance_rescaled_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval f) :comments org :exports none
   integer(c_int32_t) function qmckl_distance_rescaled &
       (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) &
       bind(C) result(info)

     use, intrinsic :: iso_c_binding
     implicit none

     integer (c_int64_t) , intent(in)  , value :: context
     character           , intent(in)  , value :: transa
     character           , intent(in)  , value :: transb
     integer (c_int64_t) , intent(in)  , value :: m
     integer (c_int64_t) , intent(in)  , value :: n
     real    (c_double ) , intent(in)          :: A(lda,*)
     integer (c_int64_t) , intent(in)  , value :: lda
     real    (c_double ) , intent(in)          :: B(ldb,*)
     integer (c_int64_t) , intent(in)  , value :: ldb
     real    (c_double ) , intent(out)         :: C(ldc,n)
     integer (c_int64_t) , intent(in)  , value :: ldc
     real    (c_double ) , intent(in)  , value :: rescale_factor_kappa

     integer(c_int32_t), external :: qmckl_distance_rescaled_f
     info = qmckl_distance_rescaled_f &
            (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa)

   end function qmckl_distance_rescaled
   #+end_src

   #+CALL: generate_f_interface(table=qmckl_distance_rescaled_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
     integer(c_int32_t) function qmckl_distance_rescaled &
         (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) &
         bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

       integer (c_int64_t) , intent(in)  , value :: context
       character           , intent(in)  , value :: transa
       character           , intent(in)  , value :: transb
       integer (c_int64_t) , intent(in)  , value :: m
       integer (c_int64_t) , intent(in)  , value :: n
       real    (c_double ) , intent(in)          :: A(lda,*)
       integer (c_int64_t) , intent(in)  , value :: lda
       real    (c_double ) , intent(in)          :: B(ldb,*)
       integer (c_int64_t) , intent(in)  , value :: ldb
       real    (c_double ) , intent(out)         :: C(ldc,n)
       integer (c_int64_t) , intent(in)  , value :: ldc
       real    (c_double ) , intent(in)  , value :: rescale_factor_kappa

     end function qmckl_distance_rescaled
   end interface
   #+end_src

*** Test                                                           :noexport:
* Rescaled Distance Derivatives

** ~qmckl_distance_rescaled_deriv_e~
   :PROPERTIES:
   :Name:     qmckl_distance_rescaled_deriv_e
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

 ~qmckl_distance_rescaled_deriv_e~ computes the matrix of the gradient and laplacian of the
   rescaled distance with respect to the electron coordinates. The derivative is a rank 3 tensor.
   The first dimension has a dimension of 4 of which the first three coordinates
   contains the gradient vector and the last index is the laplacian.


   \[
   C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa
   \]

   Here the gradient is defined as follows:

   \[
   \nabla (C_{ij}(\mathbf{r}_{ee})) = \left(\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta x},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta y},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta z} \right)
   \]
   and the laplacian is defined as follows:

   \[
   \triangle (C_{ij}(r_{ee})) = \frac{\delta^2}{\delta x^2} + \frac{\delta^2}{\delta y^2} + \frac{\delta^2}{\delta z^2}
   \]

   Using the above three formulae, the expression for the gradient and laplacian is
   as follows:

   \[
   \frac{\delta  C_{ij}(\mathbf{r}_{ee})}{\delta x} = \frac{|(x_i - x_j)|}{r_{ij}} (1 - \kappa R_{ij})
   \]

   \[
   \frac{\delta  C_{ij}(\mathbf{r}_{ee})}{\delta y} = \frac{|(y_i - y_j)|}{r_{ij}} (1 - \kappa R_{ij})
   \]

   \[
   \frac{\delta  C_{ij}(\mathbf{r}_{ee})}{\delta z} = \frac{|(z_i - z_j)|}{r_{ij}} (1 - \kappa R_{ij})
   \]

   \[
   \Delta(C_{ij}(r_{ee}) = \left[ \frac{2}{r_{ij}} - \kappa  \right] (1-\kappa R_{ij})
   \]

   If the input array is normal (~'N'~), the xyz coordinates are in
   the leading dimension: ~[n][3]~ in C and ~(3,n)~ in Fortran.

   #+NAME: qmckl_distance_rescaled_deriv_e_args
   | Variable               | Type                | In/Out | Description                                           |
   |------------------------+---------------------+--------+-------------------------------------------------------|
   | ~context~              | ~qmckl_context~     | in     | Global state                                          |
   | ~transa~               | ~char~              | in     | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed         |
   | ~transb~               | ~char~              | in     | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed         |
   | ~m~                    | ~int64_t~           | in     | Number of points in the first set                     |
   | ~n~                    | ~int64_t~           | in     | Number of points in the second set                    |
   | ~A~                    | ~double[][lda]~     | in     | Array containing the $m \times 3$ matrix $A$          |
   | ~lda~                  | ~int64_t~           | in     | Leading dimension of array ~A~                        |
   | ~B~                    | ~double[][ldb]~     | in     | Array containing the $n \times 3$ matrix $B$          |
   | ~ldb~                  | ~int64_t~           | in     | Leading dimension of array ~B~                        |
   | ~C~                    | ~double[4][n][ldc]~ | out    | Array containing the $4 \times m \times n$ matrix $C$ |
   | ~ldc~                  | ~int64_t~           | in     | Leading dimension of array ~C~                        |
   | ~rescale_factor_kappa~ | ~double~            | in     | Factor for calculating rescaled distances derivatives |

   Requirements:

    - ~context~ is not ~QMCKL_NULL_CONTEXT~
    - ~m > 0~
    - ~n > 0~
    - ~lda >= 3~ if ~transa == 'N'~
    - ~lda >= m~ if ~transa == 'T'~
    - ~ldb >= 3~ if ~transb == 'N'~
    - ~ldb >= n~ if ~transb == '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 $4 \times m \times n \times 8$ bytes
      
    #+CALL: generate_c_header(table=qmckl_distance_rescaled_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_distance_rescaled_deriv_e (
          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,
          double* const C,
          const int64_t ldc,
          const double rescale_factor_kappa ); 
    #+end_src

    #+begin_src f90 :tangle (eval f)
integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , 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(4,ldc,*)
  real*8     , intent(in)  :: rescale_factor_kappa

  integer*8 :: i,j
  real*8    :: x, y, z, dist, dist_inv
  real*8    :: rescale_factor_kappa_inv, rij
  integer   :: transab

  rescale_factor_kappa_inv = 1.0d0/rescale_factor_kappa;

  info = QMCKL_SUCCESS

  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 (transb == 'T' .or. transb == 't') then
     transab = transab + 2
  else
     transab = -100
  endif

  ! check for LDA
  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

  ! check for LDB
  if (iand(transab,1) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,1) == 1 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,2) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,2) == 2 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  ! check for LDC
  if (LDC < m) then
     info = QMCKL_INVALID_ARG_11
     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)
           dist = dsqrt(x*x + y*y + z*z)
           dist_inv = 1.0d0/dist
           rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
           C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
        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)
           dist = dsqrt(x*x + y*y + z*z)
           dist_inv = 1.0d0/dist
           rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
           C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
        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)
           dist = dsqrt(x*x + y*y + z*z)
           dist_inv = 1.0d0/dist
           rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
           C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
        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)
           dist = dsqrt(x*x + y*y + z*z)
           dist_inv = 1.0d0/dist
           rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
           C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
        end do
     end do

  end select

end function qmckl_distance_rescaled_deriv_e_f
    #+end_src

    This function is more efficient when ~A~ and ~B~ are transposed.
    
    #+CALL: generate_c_interface(table=qmckl_distance_rescaled_deriv_e_args,fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval f) :comments org :exports none
   integer(c_int32_t) function qmckl_distance_rescaled_deriv_e &
       (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) &
       bind(C) result(info)

     use, intrinsic :: iso_c_binding
     implicit none

     integer (c_int64_t) , intent(in)  , value :: context
     character           , intent(in)  , value :: transa
     character           , intent(in)  , value :: transb
     integer (c_int64_t) , intent(in)  , value :: m
     integer (c_int64_t) , intent(in)  , value :: n
     real    (c_double ) , intent(in)          :: A(lda,*)
     integer (c_int64_t) , intent(in)  , value :: lda
     real    (c_double ) , intent(in)          :: B(ldb,*)
     integer (c_int64_t) , intent(in)  , value :: ldb
     real    (c_double ) , intent(out)         :: C(ldc,n,4)
     integer (c_int64_t) , intent(in)  , value :: ldc
     real    (c_double ) , intent(in)  , value :: rescale_factor_kappa

     integer(c_int32_t), external :: qmckl_distance_rescaled_deriv_e_f
     info = qmckl_distance_rescaled_deriv_e_f &
            (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa)

   end function qmckl_distance_rescaled_deriv_e
   #+end_src

   #+CALL: generate_f_interface(table=qmckl_distance_rescaled_deriv_e_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
     integer(c_int32_t) function qmckl_distance_rescaled_deriv_e &
         (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) &
         bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

       integer (c_int64_t) , intent(in)  , value :: context
       character           , intent(in)  , value :: transa
       character           , intent(in)  , value :: transb
       integer (c_int64_t) , intent(in)  , value :: m
       integer (c_int64_t) , intent(in)  , value :: n
       real    (c_double ) , intent(in)          :: A(lda,*)
       integer (c_int64_t) , intent(in)  , value :: lda
       real    (c_double ) , intent(in)          :: B(ldb,*)
       integer (c_int64_t) , intent(in)  , value :: ldb
       real    (c_double ) , intent(out)         :: C(ldc,n,4)
       integer (c_int64_t) , intent(in)  , value :: ldc
       real    (c_double ) , intent(in)  , value :: rescale_factor_kappa

     end function qmckl_distance_rescaled_deriv_e
   end interface
   #+end_src

* End of files                                                     :noexport:

   #+begin_src c :comments link :tangle (eval c_test)

  assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);

  return 0;
}

   #+end_src


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