#+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 #ifdef HAVE_CONFIG_H #include "config.h" #endif int main() { qmckl_context context; context = qmckl_context_create(); qmckl_init_probes(); #+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); qmckl_dump_probes(); return 0; } #+end_src # -*- mode: org -*- # vim: syntax=c