1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-22 20:36:01 +01:00

Added ee-distance computation

This commit is contained in:
Anthony Scemama 2021-04-26 01:45:25 +02:00
parent b56a13be66
commit f0f6ac7d85
5 changed files with 612 additions and 18 deletions

View File

@ -92,7 +92,7 @@ endif
# The Makefile.generated is the one that will be distributed with the library.
.PHONY: clean shared static doc all check install uninstall
.PHONY: clean shared static doc all check install uninstall syntax
.SECONDARY: # Needed to keep the produced C and Fortran files
$(shared_lib) $(static_lib) install uninstall: $(qmckl_h) $(qmckl_f) Makefile.generated
@ -108,6 +108,9 @@ all: shared static doc check
check: $(static_lib)
$(MAKE) -f Makefile.generated check
syntax:
cppcheck --addon=cert qmckl_*.c
doc: $(ORG_SOURCE_FILES)
$(QMCKL_ROOT)/tools/build_doc.sh

View File

@ -207,7 +207,7 @@ qmckl_context qmckl_context_create() {
ctx->numprec.range = QMCKL_DEFAULT_RANGE;
ctx->ao_basis.uninitialized = (1 << 10) - 1;
ctx->electron.uninitialized = (1 << 4) - 1;
ctx->electron.uninitialized = (1 << 2) - 1;
/* Allocate qmckl_memory_struct */
{

View File

@ -83,7 +83,9 @@ MunitResult test_<<filename()>>() {
*** Source
#+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)
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
@ -100,7 +102,7 @@ integer function qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, L
real*8 :: x, y, z
integer :: transab
info = 0
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
@ -390,6 +392,384 @@ end function test_qmckl_distance_sq
int test_qmckl_distance_sq(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_distance_sq(context));
#+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}
\]
#+NAME: qmckl_distance_args
| qmckl_context | context | in | Global state |
| char | transa | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| char | transb | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| int64_t | m | in | Number of points in the first set |
| int64_t | n | in | Number of points in the second set |
| double | A[][lda] | in | Array containing the $m \times 3$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[][ldb] | in | Array containing the $n \times 3$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
| double | C[n][ldc] | out | Array containing the $m \times n$ matrix $C$ |
| int64_t | ldc | 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 (transa == 'T' .or. transa == '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
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 might be 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
implicit none
integer(qmckl_context), intent(in), value :: context
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)
if (test_qmckl_dist == 0) return
test_qmckl_dist = &
qmckl_distance(context, 't', 'X', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_dist == 0) return
test_qmckl_dist = &
qmckl_distance(context, 'T', 't', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_dist /= 0) return
test_qmckl_dist = -1
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)
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
end do
end do
test_qmckl_dist = &
qmckl_distance(context, 'n', 'T', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_dist /= 0) return
test_qmckl_dist = -1
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)
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
end do
end do
test_qmckl_dist = &
qmckl_distance(context, 'T', 'n', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_dist /= 0) return
test_qmckl_dist = -1
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)
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
end do
end do
test_qmckl_dist = &
qmckl_distance(context, 'n', 'N', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_dist /= 0) return
test_qmckl_dist = -1
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)
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
end do
end do
test_qmckl_dist = 0
deallocate(A,B,C)
end function test_qmckl_dist
#+end_src
#+begin_src c :comments link :tangle (eval c_test)
int test_qmckl_dist(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_dist(context));
#+end_src
* End of files :noexport:
#+begin_src c :comments link :tangle (eval c_test)

View File

@ -32,6 +32,9 @@ MunitResult test_<<filename()>>() {
#include <string.h>
#include <stdbool.h>
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include "qmckl_error_type.h"
#include "qmckl_context_type.h"
@ -42,33 +45,38 @@ MunitResult test_<<filename()>>() {
#include "qmckl_memory_private_func.h"
#include "qmckl_memory_func.h"
#include "qmckl_context_func.h"
#include "qmckl_electron_private_func.h"
#+end_src
* Context
The following data stored in the context:
| ~date~ | uint64_t | Last modification date of the coordinates |
| ~uninitialized~ | int32_t | Keeps bit set for uninitialized data |
| ~num~ | int64_t | Total number of electrons |
| ~up_num~ | int64_t | Number of up-spin electrons |
| ~down_num~ | int64_t | Number of down-spin electrons |
| ~walk_num~ | int64_t | Number of walkers |
| ~provided~ | bool | If true, ~electron~ is valid |
| ~coord_new~ | double[walk_num][3][num] | New set of electron coordinates |
| ~coord_old~ | double[walk_num][3][num] | Old set of electron coordinates |
| ~uninitialized~ | int32_t | Keeps bit set for uninitialized data |
| ~num~ | int64_t | Total number of electrons |
| ~up_num~ | int64_t | Number of up-spin electrons |
| ~down_num~ | int64_t | Number of down-spin electrons |
| ~walk_num~ | int64_t | Number of walkers |
| ~provided~ | bool | If true, ~electron~ is valid |
| ~coord_new~ | double[walk_num][3][num] | New set of electron coordinates |
| ~coord_old~ | double[walk_num][3][num] | Old set of electron coordinates |
| ~coord_new_date~ | uint64_t | Last modification date of the coordinates |
| ~ee_distance~e | double[walk_num][num][num] | Electron-electron distances |
| ~ee_distance_date~ | uint64_t | Last modification date of the electron-electron distances |
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_electron_struct {
int64_t date;
int64_t num;
int64_t up_num;
int64_t down_num;
int64_t walk_num;
int64_t coord_new_date;
int64_t ee_distance_date;
double* coord_new;
double* coord_old;
double* ee_distance;
int32_t uninitialized;
bool provided;
} qmckl_electron_struct;
@ -252,6 +260,7 @@ if (ctx->electron.provided) {
NULL);
}
ctx->electron.coord_old = coord_old;
}
return QMCKL_SUCCESS;
@ -280,7 +289,7 @@ qmckl_exit_code qmckl_set_electron_num(qmckl_context context,
"down_num <= 0");
}
int32_t mask = (1 << 3) -1;
int32_t mask = 1;
ctx->electron.up_num = up_num;
ctx->electron.down_num = down_num;
@ -302,7 +311,7 @@ qmckl_exit_code qmckl_set_electron_walk_num(qmckl_context context, const int64_t
"walk_num <= 0");
}
int32_t mask = 1 << 3;
int32_t mask = 2;
ctx->electron.walk_num = walk_num;
<<post2>>
@ -350,7 +359,7 @@ qmckl_exit_code qmckl_set_electron_coord(qmckl_context context, const double* c
ctx->electron.coord_new = swap;
memcpy(ctx->electron.coord_new, coord, walk_num * num * 3 * sizeof(double));
ctx->electron.date = ctx->date;
ctx->electron.coord_new_date = ctx->date;
return QMCKL_SUCCESS;
@ -398,6 +407,207 @@ munit_assert_int64(rc, ==, QMCKL_SUCCESS);
#+end_src
* Computation
The computed data is stored in the context so that it can be reused
by different kernels. To ensure that the data is valid, for each
computed data the date of the context is stored when it is computed.
To know if some data needs to be recomputed, we check if the date of
the dependencies are more recent than the date of the data to
compute. If it is the case, then the data is recomputed and the
current date is stored.
** Electron-electron distances
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_electron_ee_distance(qmckl_context context, double* distance);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_electron_ee_distance(qmckl_context context, double* distance)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
}
qmckl_exit_code rc = qmckl_provide_ee_distance(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
memcpy(distance, ctx->electron.ee_distance, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->electron.ee_distance_date) {
fprintf(stderr, "%10ld: provide ee_distance", ctx->date);
/* Allocate array */
if (ctx->electron.ee_distance == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.num * ctx->electron.num *
ctx->electron.walk_num * sizeof(double);
double* ee_distance = (double*) qmckl_malloc(context, mem_info);
if (ee_distance == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_ee_distance",
NULL);
}
ctx->electron.ee_distance = ee_distance;
}
qmckl_exit_code rc =
qmckl_compute_ee_distance(context,
ctx->electron.num,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->electron.ee_distance);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->electron.ee_distance_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute
:PROPERTIES:
:Name: qmckl_compute_ee_distance
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_ee_distance_args
| qmckl_context | context | in | Global state |
| int64_t | elec_num | in | Number of electrons |
| int64_t | walk_num | in | Number of walkers |
| double | coord[walk_num][3][elec_num] | in | Electron coordinates |
| double | ee_distance[walk_num][elec_num][elec_num] | out | Electron-electron distances |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ee_distance_f(context, elec_num, walk_num, coord, ee_distance) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: elec_num
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: coord(elec_num,3,walk_num)
double precision , intent(out) :: ee_distance(elec_num,elec_num,walk_num)
integer*8 :: k
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP SHARED(elec_num, walk_num, coord, ee_distance)
!$OMP PRIVATE(k)
do k=1,walk_num
info = qmckl_distance(context, 'T', 'T', elec_num, elec_num, &
coord(1,1,k), elec_num, &
coord(1,1,k), elec_num, &
ee_distance(1,1,k), elec_num)
end do
!$OMP END PARALLEL DO
end function qmckl_compute_ee_distance_f
#+end_src
#+begin_src c :tangle (eval h_private_func) :comments org :exports none
qmckl_exit_code qmckl_compute_ee_distance (
const qmckl_context context,
const int64_t elec_num,
const int64_t walk_num,
const double* coord,
double* const ee_distance );
#+end_src
#+CALL: generate_c_interface(table=qmckl_ee_distance_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_ee_distance &
(context, elec_num, walk_num, coord, ee_distance) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) , value :: walk_num
real (c_double ) , intent(in) :: coord(elec_num,3,walk_num)
real (c_double ) , intent(out) :: ee_distance(elec_num,elec_num,walk_num)
integer(c_int32_t), external :: qmckl_compute_ee_distance_f
info = qmckl_compute_ee_distance_f &
(context, elec_num, walk_num, coord, ee_distance)
end function qmckl_compute_ee_distance
#+end_src
*** Test
#+begin_src c :tangle (eval c_test)
/* Reference input data */
munit_assert(qmckl_electron_provided(context));
double distance[walk_num*num*num];
rc = qmckl_get_electron_ee_distance(context, distance);
rc = qmckl_get_electron_ee_distance(context, distance);
munit_assert_double(distance[0], ==, 0.);
munit_assert_double(distance[1], ==, distance[num]);
munit_assert_double_equal(distance[1], 8.6114953086801, 12);
#+end_src
* End of files :noexport:
#+begin_src c :tangle (eval h_private_type)

View File

@ -39,6 +39,7 @@ MunitResult test_<<filename()>>() {
#define QMCKL_MEMORY_HPT
#include <stdint.h>
#include <malloc.h>
#+end_src
* Memory data structure for the context