TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-11-19 20:42:50 +01:00
Code Issues Projects Releases Wiki Activity
1b02b49f51
qmckl/org/qmckl_point.org
Anthony Scemama 60ef0dc7a6 Redesign walkers
2022-08-07 15:28:00 +02:00

477 lines
13 KiB
Org Mode
Raw Blame History

#+TITLE: Point
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
This data structure contains cartesian coordinates where the
3-dimensional functions will be evaluated. For DFT codes these may be
the integration grid points. For QMC codes, these are the electron
coordinates of all the walkers.
* Headers :noexport:
#+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src
#+begin_src c :tangle (eval h_private_type)
#ifndef QMCKL_POINT_HPT
#define QMCKL_POINT_HPT
#include <stdbool.h>
#include "qmckl_blas_private_type.h"
#+end_src
#+begin_src c :tangle (eval h_private_func)
#ifndef QMCKL_POINT_HPF
#define QMCKL_POINT_HPF
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
#+end_src
#+begin_src c :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include <assert.h>
#include <math.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "chbrclf.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
int main() {
qmckl_context context;
context = qmckl_context_create();
#+end_src
#+begin_src c :tangle (eval c)
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_STDINT_H
#include <stdint.h>
#elif HAVE_INTTYPES_H
#include <inttypes.h>
#endif
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_point_private_func.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_func.h"
#+end_src
* Context
The following data stored in the context:
| Variable | Type | Description |
|--------------+----------------+-------------------------------------------|
| ~num~ | ~int64_t~ | Total number of points |
| ~date~ | ~uint64_t~ | Last modification date of the coordinates |
| ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix |
We consider that the matrix is stored 'transposed' and 'normal'
corresponds to the 3 \times ~num~ matrix.
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_point_struct {
int64_t num;
uint64_t date;
qmckl_matrix coord;
} qmckl_point_struct;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_point(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_point(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
memset(&(ctx->point), 0, sizeof(qmckl_point_struct));
return QMCKL_SUCCESS;
}
#+end_src
** Access functions
Access functions return ~QMCKL_SUCCESS~ when the data has been
successfully retrieved. They return ~QMCKL_INVALID_CONTEXT~ when
the context is not a valid context. If the function returns
successfully, the variable pointed by the pointer given in argument
contains the requested data. Otherwise, this variable is untouched.
*** Number of points
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_point_num (const qmckl_context context, int64_t* const num);
#+end_src
Returns the number of points stored in the context.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_point_num (const qmckl_context context, int64_t* const num) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (num == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_point_num",
"num is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
assert (ctx->point.num >= (int64_t) 0);
,*num = ctx->point.num;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_get_point_num(context, num) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(out) :: num
end function
end interface
#+end_src
*** Point coordinates
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_point(const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max);
#+end_src
Returns the point coordinates as sequences of (x,y,z).
The pointer is assumed to point on a memory block of size ~size_max~ \ge ~3 * point_num~.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_point(const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (coord == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_point_coord",
"coord is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t point_num = ctx->point.num;
if (point_num == 0) return QMCKL_NOT_PROVIDED;
assert (ctx->point.coord.data != NULL);
if (size_max < 3*point_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_point_coord",
"size_max too small");
}
qmckl_exit_code rc;
if (transp == 'N') {
qmckl_matrix At = qmckl_matrix_alloc( context, 3, point_num);
rc = qmckl_transpose( context, ctx->point.coord, At );
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_double_of_matrix( context, At, coord, size_max);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_matrix_free(context, &At);
} else {
rc = qmckl_double_of_matrix( context, ctx->point.coord, coord, size_max);
}
if (rc != QMCKL_SUCCESS) return rc;
return rc;
}
#+end_src
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_get_point(context, transp, coord, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character(c_char) , intent(in) , value :: transp
real (c_double ) , intent(out) :: coord(*)
integer (c_int64_t) , intent(in) :: size_max
end function
end interface
#+end_src
** Initialization functions
When the data is set in the context, if the arrays are large
enough, we overwrite the data contained in them.
To set the data relative to the points in the context, one of the
following functions need to be called. Here, ~num~ is the number of
points to set.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_point (qmckl_context context,
const char transp,
const int64_t num,
const double* coord,
const int64_t size_max);
#+end_src
Copy a sequence of ~num~ points $(x,y,z)$ into the context.
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code
qmckl_set_point (qmckl_context context,
const char transp,
const int64_t num,
const double* coord,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
if (size_max < 3*num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_point",
"Array too small");
}
if (transp != 'N' && transp != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_set_point",
"transp should be 'N' or 'T'");
}
if (coord == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_set_point",
"coord is a NULL pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_exit_code rc;
if (num != ctx->point.num) {
if (ctx->point.coord.data != NULL) {
rc = qmckl_matrix_free(context, &(ctx->point.coord));
assert (rc == QMCKL_SUCCESS);
}
ctx->point.coord = qmckl_matrix_alloc(context, num, 3);
if (ctx->point.coord.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
NULL);
}
};
ctx->point.num = num;
if (transp == 'T') {
double *a = ctx->point.coord.data;
#ifdef HAVE_OPENMP
#pragma omp for
#endif
for (int64_t i=0 ; i<3*num ; ++i) {
a[i] = coord[i];
}
} else {
#ifdef HAVE_OPENMP
#pragma omp for
#endif
for (int64_t i=0 ; i<num ; ++i) {
qmckl_mat(ctx->point.coord, i, 0) = coord[3*i ];
qmckl_mat(ctx->point.coord, i, 1) = coord[3*i+1];
qmckl_mat(ctx->point.coord, i, 2) = coord[3*i+2];
}
}
/* Increment the date of the context */
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_set_point(context, &
transp, num, coord, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character(c_char) , intent(in) , value :: transp
integer (c_int64_t) , intent(in) , value :: num
real (c_double ) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
/* Reference input data */
int64_t point_num = chbrclf_elec_num;
const double* coord = &(chbrclf_elec_coord[0][0][0]);
/* --- */
qmckl_exit_code rc;
double coord2[point_num*3];
double coord3[point_num*3];
rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_point (context, 'N', point_num, coord, (point_num*3));
assert(rc == QMCKL_SUCCESS);
int64_t n;
rc = qmckl_get_point_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == point_num);
rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*point_num ; ++i) {
assert( coord[i] == coord2[i] );
}
rc = qmckl_get_point (context, 'T', coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<point_num ; ++i) {
assert( coord[3*i+0] == coord2[i] );
assert( coord[3*i+1] == coord2[i+point_num] );
assert( coord[3*i+2] == coord2[i+point_num*2] );
}
rc = qmckl_set_point (context, 'T', point_num, coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_point (context, 'N', coord3, (point_num*3));
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*point_num ; ++i) {
assert( coord[i] == coord3[i] );
}
#+end_src
* End of files :noexport:
#+begin_src c :tangle (eval h_private_type)
#endif
#+end_src
#+begin_src c :tangle (eval h_private_func)
#endif
#+end_src
*** Test
#+begin_src c :tangle (eval c_test)
if (qmckl_context_destroy(context) != QMCKL_SUCCESS)
return QMCKL_FAILURE;
return 0;
}
#+end_src
*** Compute file names
#+begin_src emacs-lisp
; The following is required to compute the file names
(setq pwd (file-name-directory buffer-file-name))
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
(setq f (concat pwd name "_f.f90"))
(setq fh (concat pwd name "_fh.f90"))
(setq c (concat pwd name ".c"))
(setq h (concat name ".h"))
(setq h_private (concat name "_private.h"))
(setq c_test (concat pwd "test_" name ".c"))
(setq f_test (concat pwd "test_" name "_f.f90"))
; Minted
(require 'ox-latex)
(setq org-latex-listings 'minted)
(add-to-list 'org-latex-packages-alist '("" "listings"))
(add-to-list 'org-latex-packages-alist '("" "color"))
#+end_src
#+RESULTS:
| | color |
| | listings |
# -*- mode: org -*-
# vim: syntax=c
Reference in New Issue View Git Blame Copy Permalink