mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-08-15 01:48:29 +02:00
1548 lines
47 KiB
Org Mode
1548 lines
47 KiB
Org Mode
#+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();
|
|
|
|
#ifdef VFC_CI
|
|
qmckl_init_probes();
|
|
#endif
|
|
|
|
#+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)
|
|
function qmckl_distance_sq(context, transa, transb, m, n, &
|
|
A, LDA, B, LDB, C, LDC) &
|
|
bind(C) result(info)
|
|
use, intrinsic :: iso_c_binding
|
|
|
|
use qmckl_constants
|
|
implicit none
|
|
|
|
integer (qmckl_context) , intent(in) , value :: context
|
|
character(c_char) , intent(in) , value :: transa
|
|
character(c_char) , 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(qmckl_exit_code) :: info
|
|
|
|
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. LDB < 3) then
|
|
info = QMCKL_INVALID_ARG_7
|
|
return
|
|
endif
|
|
|
|
if (iand(transab,2) == 2 .and. LDB < n) 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
|
|
#+end_src
|
|
|
|
*** Performance
|
|
|
|
This function is more efficient when ~A~ and ~B~ are
|
|
transposed.
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_distance_sq_args,fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(qmckl_exit_code) 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 (qmckl_context), intent(in) , value :: context
|
|
character(c_char ) , intent(in) , value :: transa
|
|
character(c_char ) , 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", "distance_sq_Xt_2_2", 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", "distance_sq_tX_2_2", 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", "distance_sq_Tt_2_2", C(2,2), 0.d0, 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", "distance_sq_nT_2_2", C(2,2), 0.d0, 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", "distance_sq_Tn_2_2", C(2,2), 0.d0, 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", "distance_sq_nN_2_2", C(2,2), 0.d0, 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)
|
|
function qmckl_distance(context, transa, transb, m, n, &
|
|
A, LDA, B, LDB, C, LDC) &
|
|
bind(C) result(info)
|
|
use, intrinsic :: iso_c_binding
|
|
use qmckl_constants
|
|
implicit none
|
|
integer(qmckl_context), intent(in), value :: context
|
|
character(c_char) , intent(in) , value :: transa
|
|
character(c_char) , 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 (qmckl_exit_code) :: info
|
|
|
|
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
|
|
|
|
! 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
|
|
|
|
! 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
|
|
#+end_src
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_distance_args,fname="qmckl_distance")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(qmckl_exit_code) 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 (qmckl_context), intent(in) , value :: context
|
|
character(c_char ) , intent(in) , value :: transa
|
|
character(c_char ) , 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
|
|
|
|
*** Performance
|
|
|
|
This function is more efficient when ~A~ and ~B~ are transposed.
|
|
|
|
*** 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", "distance_Xt_2_2", 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", "distance_tX_2_2", 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", "distance_Tt_2_2", C(2,2), 0.d0, 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", "distance_nT_2_2", C(2,2), 0.d0, 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", "distance_Tn_2_2", C(2,2), 0.d0, 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", "distance_nN_2_2", C(2,2), 0.d0, 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} = \frac{ 1 - e^{-\kappa r_{ij}}}{\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)
|
|
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
|
|
use qmckl_constants
|
|
implicit none
|
|
integer (qmckl_context), intent(in) , value :: context
|
|
character(c_char ) , intent(in) , value :: transa
|
|
character(c_char ) , 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(qmckl_exit_code) :: info
|
|
|
|
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
|
|
|
|
! check for LDB
|
|
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
|
|
#+end_src
|
|
|
|
*** Performance
|
|
|
|
This function is more efficient when ~A~ and ~B~ are transposed.
|
|
|
|
** C interface :noexport:
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_distance_rescaled_args,fname="qmckl_distance_rescaled")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(qmckl_exit_code) 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 (qmckl_context), intent(in) , value :: context
|
|
character(c_char ) , intent(in) , value :: transa
|
|
character(c_char ) , 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:
|
|
|
|
#+BEGIN_SRC python :results output :exports none
|
|
import numpy as np
|
|
|
|
kappa = 0.6
|
|
kappa_inv = 1./kappa
|
|
|
|
# For H2O we have the following data:
|
|
elec_num = 10
|
|
nucl_num = 2
|
|
up_num = 5
|
|
down_num = 5
|
|
nucl_coord = np.array([ [0.000000, 0.000000 ],
|
|
[0.000000, 0.000000 ],
|
|
[0.000000, 2.059801 ] ])
|
|
|
|
elec_coord = np.array( [[[-0.250655104764153 , 0.503070975550133 , -0.166554344502303],
|
|
[-0.587812193472177 , -0.128751981129274 , 0.187773606533075],
|
|
[ 1.61335569047166 , -0.615556732874863 , -1.43165470979934 ],
|
|
[-4.901239896295210E-003 , -1.120440036458986E-002 , 1.99761909330422 ],
|
|
[ 0.766647499681200 , -0.293515395797937 , 3.66454589201239 ],
|
|
[-0.127732483187947 , -0.138975497694196 , -8.669850480215846E-002],
|
|
[-0.232271834949124 , -1.059321673434182E-002 , -0.504862241464867],
|
|
[ 1.09360863531826 , -2.036103063808752E-003 , -2.702796910818986E-002],
|
|
[-0.108090166832043 , 0.189161729653261 , 2.15398313919894],
|
|
[ 0.397978144318712 , -0.254277292595981 , 2.54553335476344]]])
|
|
|
|
ee_distance_rescaled = \
|
|
np.array([ [(1.-np.exp(-kappa*np.linalg.norm(elec_coord[0,j,:]-elec_coord[0,i,:])))/kappa \
|
|
for i in range(elec_num) ]
|
|
for j in range(elec_num) ])
|
|
|
|
en_distance_rescaled = \
|
|
np.array([ [(1.-np.exp(-kappa*np.linalg.norm(elec_coord[0,j,:]-nucl_coord[:,i])))/kappa \
|
|
for j in range(elec_num) ]
|
|
for i in range(nucl_num) ])
|
|
|
|
print(ee_distance_rescaled)
|
|
print()
|
|
print(en_distance_rescaled)
|
|
#+END_SRC
|
|
|
|
#+RESULTS:
|
|
#+begin_example
|
|
[[0. 0.63475074 1.29816415 1.23147027 1.51933127 0.54402406
|
|
0.51452479 0.96538731 1.25878564 1.3722995 ]
|
|
[0.63475074 0. 1.35148664 1.13524156 1.48940503 0.4582292
|
|
0.62758076 1.06560856 1.179133 1.30763703]
|
|
[1.29816415 1.35148664 0. 1.50021375 1.59200788 1.23488312
|
|
1.20844259 1.0355537 1.52064535 1.53049239]
|
|
[1.23147027 1.13524156 1.50021375 0. 1.12016142 1.19158954
|
|
1.29762585 1.24824277 0.25292267 0.58609336]
|
|
[1.51933127 1.48940503 1.59200788 1.12016142 0. 1.50217017
|
|
1.54012828 1.48753895 1.10441805 0.84504381]
|
|
[0.54402406 0.4582292 1.23488312 1.19158954 1.50217017 0.
|
|
0.39417354 0.87009603 1.23838502 1.33419121]
|
|
[0.51452479 0.62758076 1.20844259 1.29762585 1.54012828 0.39417354
|
|
0. 0.95118809 1.33068934 1.41097406]
|
|
[0.96538731 1.06560856 1.0355537 1.24824277 1.48753895 0.87009603
|
|
0.95118809 0. 1.29422213 1.33222493]
|
|
[1.25878564 1.179133 1.52064535 0.25292267 1.10441805 1.23838502
|
|
1.33068934 1.29422213 0. 0.62196802]
|
|
[1.3722995 1.30763703 1.53049239 0.58609336 0.84504381 1.33419121
|
|
1.41097406 1.33222493 0.62196802 0. ]]
|
|
|
|
[[0.49421587 0.52486545 1.23280503 1.16396998 1.49156627 0.1952946
|
|
0.4726453 0.80211227 1.21198526 1.31411513]
|
|
[1.24641375 1.15444238 1.50565145 0.06218339 1.10153184 1.20919677
|
|
1.3111856 1.26122875 0.22122563 0.55669168]]
|
|
#+end_example
|
|
|
|
|
|
#+begin_src f90 :tangle (eval f_test)
|
|
integer(qmckl_exit_code) function test_qmckl_dist_rescaled(context) bind(C)
|
|
|
|
use qmckl
|
|
use iso_c_binding
|
|
|
|
implicit none
|
|
|
|
integer(qmckl_context), intent(in), value :: context
|
|
integer*8 :: m, n, LDA, LDB, LDC
|
|
double precision :: x
|
|
integer*8 :: i,j
|
|
|
|
double precision, parameter :: kappa = 0.6d0
|
|
double precision, parameter :: kappa_inv = 1.d0/kappa
|
|
|
|
integer*8, parameter :: elec_num = 10_8
|
|
integer*8, parameter :: nucl_num = 2_8
|
|
|
|
double precision :: nucl_coord(nucl_num,3) = reshape( (/ &
|
|
0.0d0, 0.0d0 , &
|
|
0.0d0, 0.0d0 , &
|
|
0.0d0, 2.059801d0 /), shape(nucl_coord) )
|
|
|
|
double precision :: elec_coord(3,elec_num) = reshape( (/ &
|
|
-0.250655104764153d0 , 0.503070975550133d0 , -0.166554344502303d0 , &
|
|
-0.587812193472177d0 , -0.128751981129274d0 , 0.187773606533075d0 , &
|
|
1.61335569047166d0 , -0.615556732874863d0 , -1.43165470979934d0 , &
|
|
-4.901239896295210d-003 , -1.120440036458986d-002 , 1.99761909330422d0 , &
|
|
0.766647499681200d0 , -0.293515395797937d0 , 3.66454589201239d0 , &
|
|
-0.127732483187947d0 , -0.138975497694196d0 , -8.669850480215846d-002 , &
|
|
-0.232271834949124d0 , -1.059321673434182d-002 , -0.504862241464867d0 , &
|
|
1.09360863531826d0 , -2.036103063808752d-003 , -2.702796910818986d-002 , &
|
|
-0.108090166832043d0 , 0.189161729653261d0 , 2.15398313919894d0 , &
|
|
0.397978144318712d0 , -0.254277292595981d0 , 2.54553335476344d0 /), &
|
|
shape(elec_coord))
|
|
|
|
double precision :: ref_ee(elec_num,elec_num) = reshape( (/ &
|
|
0.d0, 0.63475074d0, 1.29816415d0, 1.23147027d0, 1.51933127d0, &
|
|
0.54402406d0, 0.51452479d0, 0.96538731d0, 1.25878564d0, 1.3722995d0 , &
|
|
0.63475074d0, 0.d0, 1.35148664d0, 1.13524156d0, 1.48940503d0, &
|
|
0.4582292d0, 0.62758076d0, 1.06560856d0, 1.179133d0, 1.30763703d0 , &
|
|
1.29816415d0, 1.35148664d0, 0.d0, 1.50021375d0, 1.59200788d0, &
|
|
1.23488312d0, 1.20844259d0, 1.0355537d0, 1.52064535d0, 1.53049239d0 , &
|
|
1.23147027d0, 1.13524156d0, 1.50021375d0, 0.d0, 1.12016142d0, &
|
|
1.19158954d0, 1.29762585d0, 1.24824277d0, 0.25292267d0, 0.58609336d0 , &
|
|
1.51933127d0, 1.48940503d0, 1.59200788d0, 1.12016142d0, 0.d0, &
|
|
1.50217017d0, 1.54012828d0, 1.48753895d0, 1.10441805d0, 0.84504381d0 , &
|
|
0.54402406d0, 0.4582292d0, 1.23488312d0, 1.19158954d0, 1.50217017d0, &
|
|
0.d0, 0.39417354d0, 0.87009603d0, 1.23838502d0, 1.33419121d0 , &
|
|
0.51452479d0, 0.62758076d0, 1.20844259d0, 1.29762585d0, 1.54012828d0, &
|
|
0.39417354d0, 0.d0, 0.95118809d0, 1.33068934d0, 1.41097406d0 , &
|
|
0.96538731d0, 1.06560856d0, 1.0355537d0, 1.24824277d0, 1.48753895d0, &
|
|
0.87009603d0, 0.95118809d0, 0.d0, 1.29422213d0, 1.33222493d0 , &
|
|
1.25878564d0, 1.179133d0, 1.52064535d0, 0.25292267d0, 1.10441805d0, &
|
|
1.23838502d0, 1.33068934d0, 1.29422213d0, 0.d0, 0.62196802d0 , &
|
|
1.3722995d0, 1.30763703d0, 1.53049239d0, 0.58609336d0, 0.84504381d0, &
|
|
1.33419121d0, 1.41097406d0, 1.33222493d0, 0.62196802d0, 0.d0 /), shape(ref_ee) )
|
|
|
|
double precision :: ref_en(elec_num, nucl_num) = reshape( (/ &
|
|
0.49421587d0, 0.52486545d0, 1.23280503d0, 1.16396998d0, 1.49156627d0, &
|
|
0.1952946d0, 0.4726453d0, 0.80211227d0, 1.21198526d0, 1.31411513d0, &
|
|
1.24641375d0, 1.15444238d0, 1.50565145d0, 0.06218339d0, 1.10153184d0, &
|
|
1.20919677d0, 1.3111856d0, 1.26122875d0, 0.22122563d0, 0.55669168d0 /), shape(ref_en) )
|
|
|
|
double precision, allocatable :: distance_ee(:,:), distance_en(:,:)
|
|
|
|
allocate( distance_ee(elec_num,elec_num), distance_en(elec_num,nucl_num) )
|
|
|
|
print *, 'ee'
|
|
test_qmckl_dist_rescaled = &
|
|
qmckl_distance_rescaled(context, 'N', 'N', elec_num, elec_num, elec_coord, &
|
|
size(elec_coord,1)*1_8, elec_coord, size(elec_coord,1)*1_8, &
|
|
distance_ee, size(distance_ee,1)*1_8, kappa)
|
|
|
|
if (test_qmckl_dist_rescaled /= QMCKL_SUCCESS) return
|
|
|
|
test_qmckl_dist_rescaled = QMCKL_FAILURE
|
|
|
|
do j=1,elec_num
|
|
do i=1,elec_num
|
|
print *, i,j,real(distance_ee(i,j)), real(ref_ee(i,j))
|
|
if (dabs(distance_ee(i,j) - ref_ee(i,j)) > 1.d-7) then
|
|
return
|
|
endif
|
|
end do
|
|
end do
|
|
|
|
print *, 'en'
|
|
test_qmckl_dist_rescaled = &
|
|
qmckl_distance_rescaled(context, 'N', 'T', elec_num, nucl_num, elec_coord, &
|
|
size(elec_coord,1)*1_8, nucl_coord, size(nucl_coord,1)*1_8, &
|
|
distance_en, size(distance_en,1)*1_8, kappa)
|
|
|
|
if (test_qmckl_dist_rescaled /= QMCKL_SUCCESS) return
|
|
|
|
test_qmckl_dist_rescaled = QMCKL_FAILURE
|
|
|
|
do j=1,nucl_num
|
|
do i=1,elec_num
|
|
print *, i,j,real(distance_en(i,j)), real(ref_en(i,j))
|
|
if (dabs(distance_en(i,j) - ref_en(i,j)) > 1.d-7) then
|
|
return
|
|
endif
|
|
end do
|
|
end do
|
|
|
|
test_qmckl_dist_rescaled = QMCKL_SUCCESS
|
|
end function test_qmckl_dist_rescaled
|
|
#+end_src
|
|
|
|
#+begin_src c :comments link :tangle (eval c_test)
|
|
qmckl_exit_code test_qmckl_dist_rescaled(qmckl_context context);
|
|
assert(test_qmckl_dist_rescaled(context) == QMCKL_SUCCESS);
|
|
#+end_src
|
|
* Rescaled Distance Derivatives
|
|
|
|
** ~qmckl_distance_rescaled_gl~
|
|
:PROPERTIES:
|
|
:Name: qmckl_distance_rescaled_gl
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END: ~qmckl_distance_rescaled_gl~ 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(r_{ij}) = \left( 1 - \exp(-\kappa\, r_{ij})\right)/\kappa
|
|
\]
|
|
|
|
Here the gradient is defined as follows:
|
|
|
|
\[
|
|
\nabla_i C(r_{ij}) = \left(\frac{\partial C(r_{ij})}{\partial x_i},\frac{\partial C(r_{ij})}{\partial y_i},\frac{\partial C(r_{ij})}{\partial z_i} \right)
|
|
\]
|
|
and the Laplacian is defined as follows:
|
|
|
|
\[
|
|
\Delta_i C(r_{ij}) = \frac{\partial^2}{\partial x_i^2} + \frac{\partial^2}{\partial y_i^2} + \frac{\partial^2}{\partial z_i^2}
|
|
\]
|
|
|
|
Using the above three formulas, the expression for the gradient and Laplacian is:
|
|
|
|
\[
|
|
\frac{\partial C(r_{ij})}{\partial x_i} = \frac{|(x_i -
|
|
x_j)|}{r_{ij}} \exp (- \kappa \, r_{ij})
|
|
\]
|
|
|
|
\[
|
|
\Delta C_{ij}(r_{ij}) = \left[ \frac{2}{r_{ij}} - \kappa \right] \exp (- \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_gl_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][4]~ | 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_gl_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_gl (
|
|
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)
|
|
function qmckl_distance_rescaled_gl(context, transa, transb, m, n, &
|
|
A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
|
|
bind(C) result(info)
|
|
use qmckl_constants
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer(qmckl_exit_code) :: info
|
|
integer (qmckl_context), intent(in) , value :: context
|
|
character(c_char ) , intent(in) , value :: transa
|
|
character(c_char ) , 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(4,ldc,n)
|
|
integer (c_int64_t) , intent(in) , value :: ldc
|
|
real (c_double ) , intent(in) , value :: rescale_factor_kappa
|
|
|
|
integer*8 :: i,j
|
|
real*8 :: x, y, z, dist, dist_inv
|
|
real*8 :: rij
|
|
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
|
|
|
|
! 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
|
|
|
|
! check for LDB
|
|
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 = max(1.d-20, dsqrt(x*x + y*y + z*z))
|
|
dist_inv = 1.0d0/dist
|
|
rij = dexp(-rescale_factor_kappa * dist)
|
|
C(1,i,j) = x * dist_inv * rij
|
|
C(2,i,j) = y * dist_inv * rij
|
|
C(3,i,j) = z * dist_inv * rij
|
|
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * 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 = max(1.d-20, dsqrt(x*x + y*y + z*z))
|
|
dist_inv = 1.0d0/dist
|
|
rij = dexp(-rescale_factor_kappa * dist)
|
|
C(1,i,j) = x * dist_inv * rij
|
|
C(2,i,j) = y * dist_inv * rij
|
|
C(3,i,j) = z * dist_inv * rij
|
|
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * 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 = max(1.d-20, dsqrt(x*x + y*y + z*z))
|
|
dist_inv = 1.0d0/dist
|
|
rij = dexp(-rescale_factor_kappa * dist)
|
|
C(1,i,j) = x * dist_inv * rij
|
|
C(2,i,j) = y * dist_inv * rij
|
|
C(3,i,j) = z * dist_inv * rij
|
|
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * 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 = max(1.d-20, dsqrt(x*x + y*y + z*z))
|
|
dist_inv = 1.0d0/dist
|
|
rij = dexp(-rescale_factor_kappa * dist)
|
|
C(1,i,j) = x * dist_inv * rij
|
|
C(2,i,j) = y * dist_inv * rij
|
|
C(3,i,j) = z * dist_inv * rij
|
|
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
|
|
end do
|
|
end do
|
|
|
|
end select
|
|
|
|
end function qmckl_distance_rescaled_gl
|
|
#+end_src
|
|
|
|
This function is more efficient when ~A~ and ~B~ are transposed.
|
|
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_distance_rescaled_gl_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(qmckl_exit_code) function qmckl_distance_rescaled_gl &
|
|
(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 (qmckl_context), intent(in) , value :: context
|
|
character(c_char ) , intent(in) , value :: transa
|
|
character(c_char ) , 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(4,ldc,n)
|
|
integer (c_int64_t) , intent(in) , value :: ldc
|
|
real (c_double ) , intent(in) , value :: rescale_factor_kappa
|
|
|
|
end function qmckl_distance_rescaled_gl
|
|
end interface
|
|
#+end_src
|
|
|
|
* End of files :noexport:
|
|
|
|
#+begin_src c :comments link :tangle (eval c_test)
|
|
|
|
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
|
|
|
|
#ifdef VFC_CI
|
|
qmckl_dump_probes();
|
|
#endif
|
|
return 0;
|
|
}
|
|
|
|
#+end_src
|
|
|
|
|
|
# -*- mode: org -*-
|
|
# vim: syntax=c
|