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

484 lines
13 KiB
Org Mode
Raw Normal View History

2022-01-20 01:50:54 +01:00
#+TITLE: Point
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
2022-08-07 14:57:10 +02:00
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.
2022-01-20 01:50:54 +01:00
* 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>
2022-01-23 16:18:46 +01:00
#include "qmckl_blas_private_type.h"
2022-01-20 01:50:54 +01:00
#+end_src
2023-03-13 15:32:35 +01:00
2022-01-20 01:50:54 +01:00
#+begin_src c :tangle (eval h_private_func)
#ifndef QMCKL_POINT_HPF
#define QMCKL_POINT_HPF
2022-01-23 16:18:46 +01:00
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
2022-01-20 01:50:54 +01:00
#+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"
2022-01-23 16:18:46 +01:00
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
2022-01-20 01:50:54 +01:00
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"
2022-01-23 16:18:46 +01:00
#include "qmckl_blas_private_type.h"
2022-01-20 01:50:54 +01:00
#include "qmckl_point_private_func.h"
2022-01-23 16:18:46 +01:00
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_func.h"
2022-01-20 01:50:54 +01:00
#+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 |
2022-01-20 01:50:54 +01:00
2022-01-23 16:18:46 +01:00
We consider that the matrix is stored 'transposed' and 'normal'
corresponds to the 3 \times ~num~ matrix.
2022-01-20 01:50:54 +01:00
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_point_struct {
2022-01-23 16:18:46 +01:00
int64_t num;
2022-01-26 17:06:51 +01:00
uint64_t date;
2022-01-23 16:18:46 +01:00
qmckl_matrix coord;
2022-01-20 01:50:54 +01:00
} qmckl_point_struct;
#+end_src
2022-01-26 17:06:51 +01:00
#+begin_src c :comments org :tangle (eval h_private_func)
2022-01-20 01:50:54 +01:00
qmckl_exit_code qmckl_init_point(qmckl_context context);
#+end_src
2022-01-26 17:06:51 +01:00
2022-01-20 01:50:54 +01:00
#+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;
}
2022-04-05 11:03:38 +02:00
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
2022-01-20 01:50:54 +01:00
assert (ctx != NULL);
2022-01-26 17:06:51 +01:00
memset(&(ctx->point), 0, sizeof(qmckl_point_struct));
2022-01-20 01:50:54 +01:00
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.
2022-01-26 17:06:51 +01:00
2022-01-20 01:50:54 +01:00
#+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");
}
2022-04-05 11:03:38 +02:00
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
2022-01-20 01:50:54 +01:00
assert (ctx != NULL);
2022-01-26 17:06:51 +01:00
assert (ctx->point.num >= (int64_t) 0);
,*num = ctx->point.num;
2022-01-20 01:50:54 +01:00
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
2022-01-26 17:06:51 +01:00
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(out) :: num
2022-01-20 01:50:54 +01:00
end function
end interface
#+end_src
*** Point coordinates
2022-01-26 17:06:51 +01:00
2022-01-20 01:50:54 +01:00
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_point(const qmckl_context context,
2022-01-23 16:18:46 +01:00
const char transp,
2022-01-20 01:50:54 +01:00
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,
2022-01-23 16:18:46 +01:00
const char transp,
2022-01-20 01:50:54 +01:00
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");
}
2022-04-05 11:03:38 +02:00
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
2022-01-20 01:50:54 +01:00
assert (ctx != NULL);
2022-01-26 17:06:51 +01:00
int64_t point_num = ctx->point.num;
2022-01-23 16:18:46 +01:00
if (point_num == 0) return QMCKL_NOT_PROVIDED;
2022-01-20 01:50:54 +01:00
2022-01-26 17:06:51 +01:00
assert (ctx->point.coord.data != NULL);
2022-01-20 01:50:54 +01:00
if (size_max < 3*point_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_point_coord",
"size_max too small");
}
2022-01-23 16:18:46 +01:00
qmckl_exit_code rc;
if (transp == 'N') {
qmckl_matrix At = qmckl_matrix_alloc( context, 3, point_num);
2022-01-26 17:06:51 +01:00
rc = qmckl_transpose( context, ctx->point.coord, At );
2022-01-23 16:18:46 +01:00
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_double_of_matrix( context, At, coord, size_max);
if (rc != QMCKL_SUCCESS) return rc;
2022-02-16 15:14:41 +01:00
rc = qmckl_matrix_free(context, &At);
2022-01-23 16:18:46 +01:00
} else {
2022-01-26 17:06:51 +01:00
rc = qmckl_double_of_matrix( context, ctx->point.coord, coord, size_max);
2022-01-20 01:50:54 +01:00
}
2022-01-23 16:18:46 +01:00
if (rc != QMCKL_SUCCESS) return rc;
2022-01-26 17:06:51 +01:00
2022-01-23 16:18:46 +01:00
return rc;
2022-01-20 01:50:54 +01:00
}
#+end_src
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
2022-01-23 16:18:46 +01:00
integer(c_int32_t) function qmckl_get_point(context, transp, coord, size_max) bind(C)
2022-01-20 01:50:54 +01:00
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
2022-01-23 16:18:46 +01:00
character(c_char) , intent(in) , value :: transp
2022-01-20 01:50:54 +01:00
real (c_double ) , intent(out) :: coord(*)
integer (c_int64_t) , intent(in) :: size_max
end function
end interface
#+end_src
2022-01-26 17:06:51 +01:00
2022-01-23 16:18:46 +01:00
** Initialization functions
2022-01-20 01:50:54 +01:00
2022-01-23 16:18:46 +01:00
When the data is set in the context, if the arrays are large
enough, we overwrite the data contained in them.
2022-01-26 17:06:51 +01:00
2022-01-23 16:18:46 +01:00
To set the data relative to the points in the context, one of the
2022-02-27 23:31:52 +01:00
following functions need to be called. Here, ~num~ is the number of
points to set.
2022-01-26 17:06:51 +01:00
2022-01-23 16:18:46 +01:00
#+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,
2022-01-23 16:18:46 +01:00
const double* coord,
const int64_t size_max);
2022-01-23 16:18:46 +01:00
#+end_src
2022-01-26 17:06:51 +01:00
Copy a sequence of ~num~ points $(x,y,z)$ into the context.
2022-01-23 16:18:46 +01:00
#+begin_src c :comments org :tangle (eval c) :noweb yes
2022-01-20 01:50:54 +01:00
qmckl_exit_code
2022-01-23 16:18:46 +01:00
qmckl_set_point (qmckl_context context,
const char transp,
const int64_t num,
2022-01-23 16:18:46 +01:00
const double* coord,
const int64_t size_max)
2022-01-20 01:50:54 +01:00
{
2022-01-26 17:06:51 +01:00
2022-01-20 01:50:54 +01:00
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
2022-01-23 16:18:46 +01:00
return QMCKL_NULL_CONTEXT;
2022-01-20 01:50:54 +01:00
}
2022-11-18 14:08:16 +01:00
if (num <= 0) {
return qmckl_failwith( context,
2022-11-24 16:17:47 +01:00
QMCKL_INVALID_ARG_3,
2022-11-18 14:08:16 +01:00
"qmckl_set_point",
"Number of points should be >0.");
}
if (size_max < 3*num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_point",
"Array too small");
}
2022-01-23 16:18:46 +01:00
if (transp != 'N' && transp != 'T') {
2022-01-20 01:50:54 +01:00
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
2022-01-23 16:18:46 +01:00
"qmckl_set_point",
"transp should be 'N' or 'T'");
2022-01-20 01:50:54 +01:00
}
2022-01-26 17:06:51 +01:00
2022-01-23 16:18:46 +01:00
if (coord == NULL) {
2022-01-20 01:50:54 +01:00
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
2022-01-23 16:18:46 +01:00
"qmckl_set_point",
"coord is a NULL pointer");
2022-01-20 01:50:54 +01:00
}
2022-01-26 17:06:51 +01:00
2022-04-05 11:03:38 +02:00
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
2022-01-20 01:50:54 +01:00
assert (ctx != NULL);
2022-01-23 16:18:46 +01:00
qmckl_exit_code rc;
2022-08-07 14:57:10 +02:00
if (num != ctx->point.num) {
2022-01-20 01:50:54 +01:00
2022-01-26 17:06:51 +01:00
if (ctx->point.coord.data != NULL) {
2023-03-13 15:32:35 +01:00
rc = qmckl_matrix_free(context, &(ctx->point.coord));
assert (rc == QMCKL_SUCCESS);
2022-01-20 01:50:54 +01:00
}
2022-01-26 17:06:51 +01:00
ctx->point.coord = qmckl_matrix_alloc(context, num, 3);
if (ctx->point.coord.data == NULL) {
2022-01-20 01:50:54 +01:00
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
NULL);
}
};
2022-01-26 17:06:51 +01:00
ctx->point.num = num;
2022-01-20 01:50:54 +01:00
2022-01-23 16:18:46 +01:00
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];
}
2022-01-23 16:18:46 +01:00
} else {
#ifdef HAVE_OPENMP
#pragma omp for
#endif
2022-01-23 16:18:46 +01:00
for (int64_t i=0 ; i<num ; ++i) {
2022-01-26 17:06:51 +01:00
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];
2022-01-23 16:18:46 +01:00
}
2022-01-20 01:50:54 +01:00
}
2022-01-26 17:06:51 +01:00
/* Increment the date of the context */
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
2022-01-26 17:06:51 +01:00
2022-01-20 01:50:54 +01:00
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)
2022-01-20 01:50:54 +01:00
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
2022-01-23 16:18:46 +01:00
character(c_char) , intent(in) , value :: transp
2022-02-27 23:31:52 +01:00
integer (c_int64_t) , intent(in) , value :: num
2022-06-15 19:17:05 +02:00
real (c_double ) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
2022-01-20 01:50:54 +01:00
end function
end interface
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
/* Reference input data */
int64_t point_num = chbrclf_elec_num;
2022-01-23 16:18:46 +01:00
const double* coord = &(chbrclf_elec_coord[0][0][0]);
2022-01-20 01:50:54 +01:00
/* --- */
qmckl_exit_code rc;
2022-01-23 16:18:46 +01:00
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);
2022-01-20 01:50:54 +01:00
rc = qmckl_set_point (context, 'N', point_num, coord, (point_num*3));
2022-01-20 01:50:54 +01:00
assert(rc == QMCKL_SUCCESS);
int64_t n;
rc = qmckl_get_point_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == point_num);
2022-01-23 16:18:46 +01:00
rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
2022-01-20 01:50:54 +01:00
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*point_num ; ++i) {
assert( coord[i] == coord2[i] );
}
2022-01-23 16:18:46 +01:00
rc = qmckl_get_point (context, 'T', coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);
2022-01-20 01:50:54 +01:00
for (int64_t i=0 ; i<point_num ; ++i) {
2022-01-23 16:18:46 +01:00
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));
2022-01-23 16:18:46 +01:00
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] );
2022-01-20 01:50:54 +01:00
}
#+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