TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-09-16 17:55:29 +02:00
Code Issues Projects Releases Wiki Activity

Added vectors and matrices in nucleus

Browse Source
This commit is contained in:
Anthony Scemama 2022-01-23 16:18:46 +01:00
parent 61e09a7870
commit 4b36005ca0
11 changed files with 1011 additions and 802 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
org/qmckl_ao.org
View File

@ -2258,10 +2258,10 @@ qmckl_exit_code rc;
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
@ -3228,7 +3228,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
ctx->nucleus.num,
ctx->ao_basis.nucleus_prim_index,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->nucleus.coord.data,
ctx->ao_basis.exponent,
ctx->ao_basis.primitive_vgl);
} else {
@ -3296,7 +3296,7 @@ print ( "[7][4][26] : %e"% lf(a,x,y))
#+begin_src c :tangle (eval c_test) :exports none
{
#define walk_num chbrclf_walk_num
#define walk_num 1 // chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define prim_num chbrclf_prim_num
@ -3629,7 +3629,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
ctx->ao_basis.shell_prim_index,
ctx->ao_basis.shell_prim_num,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->nucleus.coord.data,
ctx->ao_basis.exponent,
ctx->ao_basis.coefficient_normalized,
ctx->ao_basis.shell_vgl);
@ -3710,7 +3710,7 @@ print ( "[1][4][26] : %25.15e"% lf(a,x,y))
#+begin_src c :tangle (eval c_test) :exports none
{
#define walk_num chbrclf_walk_num
#define walk_num 1 // chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define shell_num chbrclf_shell_num
@ -4635,7 +4635,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->nucleus.coord.data,
ctx->ao_basis.nucleus_index,
ctx->ao_basis.nucleus_shell_num,
ctx->ao_basis.nucleus_range,
@ -4735,7 +4735,7 @@ print ( "[1][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) )
#+begin_src c :tangle (eval c_test) :exports none
{
#define walk_num chbrclf_walk_num
#define walk_num 1 // chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define shell_num chbrclf_shell_num
#define ao_num chbrclf_ao_num

480
org/qmckl_blas.org
View File

@ -38,8 +38,9 @@
#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_func.h"
#+end_src
@ -379,23 +380,18 @@ qmckl_tensor_of_vector(const qmckl_vector vector,
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix,
const int64_t size);
qmckl_vector_of_matrix(const qmckl_matrix matrix);
#+end_src
Reshapes a matrix into a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix,
const int64_t size)
qmckl_vector_of_matrix(const qmckl_matrix matrix)
{
/* Always true */
assert (matrix.size[0] * matrix.size[1] == size);
qmckl_vector result;
result.size = size;
result.size = matrix.size[0] * matrix.size[1];
result.data = matrix.data;
return result;
@ -438,27 +434,23 @@ qmckl_tensor_of_matrix(const qmckl_matrix matrix,
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor,
const int64_t size);
qmckl_vector_of_tensor(const qmckl_tensor tensor);
#+end_src
Reshapes a tensor into a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor,
const int64_t size)
qmckl_vector_of_tensor(const qmckl_tensor tensor)
{
/* Always true */
int64_t prod_size = (int64_t) 1;
for (int64_t i=0 ; i<tensor.order ; i++) {
int64_t prod_size = (int64_t) tensor.size[0];
for (int64_t i=1 ; i<tensor.order ; i++) {
prod_size *= tensor.size[i];
}
assert (prod_size == size);
qmckl_vector result;
result.size = size;
result.size = prod_size;
result.data = tensor.data;
return result;
@ -510,6 +502,186 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
#define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*((l) + size[3]*(m))))]
#+end_src
** Copy to/from to ~double*~
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_double_of_vector(const qmckl_context context,
const qmckl_vector vector,
double* const target,
const int64_t size_max);
#+end_src
Converts a vector to a ~double*~.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_double_of_vector(const qmckl_context context,
const qmckl_vector vector,
double* const target,
const int64_t size_max)
{
/* Always true by construction */
assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
assert (vector.size > (int64_t) 0);
assert (target != NULL);
assert (size_max > (int64_t) 0);
assert (size_max >= vector.size);
for (int64_t i=0 ; i<vector.size ; ++i) {
target[i] = vector.data[i];
}
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_double_of_matrix(const qmckl_context context,
const qmckl_matrix matrix,
double* const target,
const int64_t size_max);
#+end_src
Converts a matrix to a ~double*~.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_double_of_matrix(const qmckl_context context,
const qmckl_matrix matrix,
double* const target,
const int64_t size_max)
{
qmckl_vector vector = qmckl_vector_of_matrix(matrix);
return qmckl_double_of_vector(context, vector, target, size_max);
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_double_of_tensor(const qmckl_context context,
const qmckl_tensor tensor,
double* const target,
const int64_t size_max);
#+end_src
Converts a tensor to a ~double*~.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_double_of_tensor(const qmckl_context context,
const qmckl_tensor tensor,
double* const target,
const int64_t size_max)
{
qmckl_vector vector = qmckl_vector_of_tensor(tensor);
return qmckl_double_of_vector(context, vector, target, size_max);
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_vector_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_vector* vector);
#+end_src
Converts a ~double*~ to a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_vector_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_vector* vector_out)
{
qmckl_vector vector = *vector_out;
/* Always true by construction */
assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
if (vector.size == 0) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_double_of_vector",
"Vector not allocated");
}
if (vector.size != size_max) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_double_of_vector",
"Wrong vector size");
}
for (int64_t i=0 ; i<vector.size ; ++i) {
vector.data[i] = target[i];
}
*vector_out = vector;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_matrix_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_matrix* matrix);
#+end_src
Converts a matrix to a ~double*~.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_matrix_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_matrix* matrix)
{
qmckl_vector vector = qmckl_vector_of_matrix(*matrix);
qmckl_exit_code rc =
qmckl_vector_of_double(context, target, size_max, &vector);
,*matrix = qmckl_matrix_of_vector(vector, matrix->size[0], matrix->size[1]);
return rc;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_tensor_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_tensor* tensor);
#+end_src
Converts a matrix to a ~double*~.
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_tensor_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_tensor* tensor)
{
qmckl_vector vector = qmckl_vector_of_tensor(*tensor);
qmckl_exit_code rc =
qmckl_vector_of_double(context, target, size_max, &vector);
*tensor = qmckl_tensor_of_vector(vector, tensor->order, tensor->size);
return rc;
}
#+end_src
** Tests
#+begin_src c :comments link :tangle (eval c_test)
@ -534,7 +706,7 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
for (int64_t i=0 ; i<m ; ++i)
assert ( qmckl_mat(mat, i, j) == qmckl_vec(vec, i+j*m)) ;
qmckl_vector vec2 = qmckl_vector_of_matrix(mat, p);
qmckl_vector vec2 = qmckl_vector_of_matrix(mat);
assert (vec2.size == p);
assert (vec2.data == vec.data);
for (int64_t i=0 ; i<p ; ++i)
@ -565,14 +737,14 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
| ~m~ | ~int64_t~ | in | Number of rows of the input matrix |
| ~n~ | ~int64_t~ | in | Number of columns of the input matrix |
| ~k~ | ~int64_t~ | in | Number of columns of the input matrix |
| ~alpha~ | ~double~ | in | Number of columns of the input matrix |
| ~alpha~ | ~double~ | in | \alpha |
| ~A~ | ~double[][lda]~ | in | Array containing the matrix $A$ |
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
| ~B~ | ~double[][ldb]~ | in | Array containing the matrix $B$ |
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~beta~ | ~double~ | in | Array containing the matrix $B$ |
| ~C~ | ~double[][ldc]~ | out | Array containing the matrix $B$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~beta~ | ~double~ | in | \beta |
| ~C~ | ~double[][ldc]~ | out | Array containing the matrix $C$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~C~ |
Requirements:
@ -808,6 +980,184 @@ qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
#+end_src
** ~qmckl_matmul~
Matrix multiplication:
\[
C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
\]
# TODO: Add description about the external library dependence.
#+NAME: qmckl_matmul_args
| Variable | Type | In/Out | Description |
|-----------+-----------------+--------+-------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~TransA~ | ~char~ | in | 'T' is transposed |
| ~TransB~ | ~char~ | in | 'T' is transposed |
| ~alpha~ | ~double~ | in | \alpha |
| ~A~ | ~qmckl_matrix~ | in | Matrix $A$ |
| ~B~ | ~qmckl_matrix~ | in | Matrix $B$ |
| ~beta~ | ~double~ | in | \beta |
| ~C~ | ~qmckl_matrix~ | out | Matrix $C$ |
#+CALL: generate_c_header(table=qmckl_matmul_args,rettyp="qmckl_exit_code",fname="qmckl_matmul")
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code
qmckl_matmul (const qmckl_context context,
const char TransA,
const char TransB,
const double alpha,
const qmckl_matrix A,
const qmckl_matrix B,
const double beta,
qmckl_matrix* const C );
#+end_src
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_matmul (const qmckl_context context,
const char TransA,
const char TransB,
const double alpha,
const qmckl_matrix A,
const qmckl_matrix B,
const double beta,
qmckl_matrix* const C )
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
qmckl_exit_code rc = QMCKL_SUCCESS;
if (TransA != 'N' && TransA != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"TransA should be 'N' or 'T'");
}
if (TransB != 'N' && TransB != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_matmul",
"TransB should be 'N' or 'T'");
}
if (A.size[0] < 1) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_5,
"qmckl_matmul",
"Invalid size for A");
}
if (B.size[0] < 1) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_6,
"qmckl_matmul",
"Invalid size for B");
}
if (C == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_8,
"qmckl_matmul",
"Null pointer");
}
int t = 0;
if (TransA == 'T') t +=1;
if (TransB == 'T') t +=2;
/*
| t | TransA | TransB |
+---+--------+--------+
| 0 | N | N |
| 1 | T | N |
| 2 | N | T |
| 3 | T | T |
,*/
switch (t) {
case 0:
if (A.size[1] != B.size[0]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[0];
C->size[1] = B.size[1];
rc = qmckl_dgemm (context, 'N', 'N',
C->size[0], C->size[1], A.size[1],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
case 1:
if (A.size[0] != B.size[0]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[1];
C->size[1] = B.size[1];
rc = qmckl_dgemm (context, 'T', 'N',
C->size[0], C->size[1], A.size[0],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
case 2:
if (A.size[1] != B.size[1]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[0];
C->size[1] = B.size[0];
rc = qmckl_dgemm (context, 'N', 'T',
C->size[0], C->size[1], A.size[1],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
case 3:
if (A.size[0] != B.size[1]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[1];
C->size[1] = B.size[0];
rc = qmckl_dgemm (context, 'T', 'T',
C->size[0], C->size[1], A.size[0],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
}
return rc;
}
#+end_src
*** TODO Test :noexport:
** ~qmckl_adjugate~
Given a matrix $\mathbf{A}$, the adjugate matrix
@ -1693,6 +2043,90 @@ qmckl_exit_code test_qmckl_adjugate(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
#+end_src
** ~qmckl_transpose~
Transposes a matrix: $A^\dagger_{ji} = A_{ij}$.
| Variable | Type | In/Out | Description |
|----------+---------------+--------+-------------------|
| context | qmckl_context | in | Global state |
| A | qmckl_matrix | in | Input matrix |
| At | qmckl_matrix | out | Transposed matrix |
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code
qmckl_transpose (qmckl_context context,
const qmckl_matrix A,
qmckl_matrix At );
#+end_src
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_transpose (qmckl_context context,
const qmckl_matrix A,
qmckl_matrix At )
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (A.size[0] < 1) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_transpose",
"Invalid size for A");
}
if (At.data == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_transpose",
"Output matrix not allocated");
}
if (At.size[0] != A.size[1] || At.size[1] != A.size[0]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_transpose",
"Invalid size for At");
}
for (int64_t j=0 ; j<At.size[1] ; ++j)
for (int64_t i=0 ; i<At.size[0] ; ++i)
qmckl_mat(At, i, j) = qmckl_mat(A, j, i);
return QMCKL_SUCCESS;
}
#+end_src
*** Test
#+begin_src c :comments link :tangle (eval c_test)
{
qmckl_matrix A;
qmckl_matrix At;
A = qmckl_matrix_alloc(context, 2, 3);
At = qmckl_matrix_alloc(context, 3, 2);
for (int j=0 ; j<3 ; ++j)
for (int i=0 ; i<2 ; ++i)
qmckl_mat(A, i, j) = (double) 10*i+j;
qmckl_exit_code rc = qmckl_transpose(context, A, At);
assert(rc == QMCKL_SUCCESS);
assert(A.size[0] == At.size[1]);
assert(A.size[1] == At.size[0]);
for (int j=0 ; j<3 ; ++j)
for (int i=0 ; i<2 ; ++i)
assert (qmckl_mat(A, i, j) == qmckl_mat(At, j, i));
qmckl_matrix_free(context, A);
qmckl_matrix_free(context, At);
}
#+end_src
* End of files :noexport:

4
org/qmckl_determinant.org
View File

@ -1160,10 +1160,10 @@ assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));

115
org/qmckl_electron.org
View File

@ -26,6 +26,7 @@ up-spin and down-spin electrons, and the electron coordinates.
#+begin_src c :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include <assert.h>
#include <stdio.h>
#include <math.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
@ -78,8 +79,8 @@ int main() {
| ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor |
| ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor |
| ~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~ | ~double[3][walk_num][num]~ | New set of electron coordinates |
| ~coord_old~ | ~double[3][walk_num][num]~ | Old set of electron coordinates |
| ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates |
Computed data:
@ -397,20 +398,27 @@ qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const
to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walk_num~.
The order of the indices is:
| | Normal | Transposed |
|---------+---------------------------+---------------------------|
| C | ~[walk_num][elec_num][3]~ | ~[walk_num][3][elec_num]~ |
| Fortran | ~(3,elec_num,walk_num)~ | ~(elec_num,3,walk_num)~ |
| | Normal | Transposed |
|---------+--------------------------+--------------------------|
| C | ~[walk_num*elec_num][3]~ | ~[3][walk_num*elec_num]~ |
| Fortran | ~(3,walk_num*elec_num)~ | ~(walk_num*elec_num, 3)~ |
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_electron_coord (const qmckl_context context, const char transp, double* const coord);
qmckl_exit_code
qmckl_get_electron_coord (const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_electron_coord (const qmckl_context context, const char transp, double* const coord) {
qmckl_get_electron_coord (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;
}
@ -442,6 +450,13 @@ qmckl_get_electron_coord (const qmckl_context context, const char transp, double
int64_t elec_num = ctx->electron.num;
int64_t walk_num = ctx->electron.walk_num;
if (size_max < elec_num*walk_num*3) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_get_electron_coord",
"Array too small");
}
assert (ctx->electron.coord_new != NULL);
double* ptr1 = ctx->electron.coord_new;
@ -449,13 +464,10 @@ qmckl_get_electron_coord (const qmckl_context context, const char transp, double
if (transp == 'N') {
for (int64_t i=0 ; i<walk_num ; ++i) {
qmckl_exit_code rc;
rc = qmckl_transpose(context, elec_num, 3,
ptr1, elec_num, ptr2, 3);
if (rc != QMCKL_SUCCESS) return rc;
ptr1 += elec_num * 3;
ptr2 += elec_num * 3;
for (int64_t i=0 ; i<elec_num*walk_num ; ++i) {
coord[3*i] = ptr1[i] ;
coord[3*i+1] = ptr1[i+elec_num*walk_num] ;
coord[3*i+2] = ptr1[i+2*elec_num*walk_num] ;
}
} else {
@ -641,6 +653,7 @@ interface
integer (c_int64_t) , intent(in) , value :: beta
end function
end interface
interface
integer(c_int32_t) function qmckl_set_electron_walk_num(context, walk_num) bind(C)
use, intrinsic :: iso_c_binding
@ -713,7 +726,7 @@ qmckl_set_electron_coord(qmckl_context context,
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_electron_coord",
"destination array is too small");
"Array too small");
}
/* If num and walk_num are set, the arrays should be allocated */
@ -732,11 +745,10 @@ qmckl_set_electron_coord(qmckl_context context,
double* ptr1 = ctx->electron.coord_new;
if (transp == 'N') {
for (int64_t i=0 ; i<walk_num ; ++i) {
rc = qmckl_transpose(context, 3, elec_num,
&(coord[i*3*elec_num]), 3, ptr1, elec_num);
if (rc != QMCKL_SUCCESS) return rc;
ptr1 += elec_num * 3;
for (int64_t i=0 ; i<walk_num*elec_num ; ++i) {
ptr1[i] = coord[3*i];
ptr1[i+walk_num*elec_num] = coord[3*i+1];
ptr1[i+2*walk_num*elec_num] = coord[3*i+2];
}
} else {
@ -865,9 +877,9 @@ assert(rc == QMCKL_SUCCESS);
double elec_coord2[walk_num*3*elec_num];
rc = qmckl_get_electron_coord (context, 'N', elec_coord2);
rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*elec_num ; ++i) {
for (int64_t i=0 ; i<3*elec_num*walk_num ; ++i) {
assert( elec_coord[i] == elec_coord2[i] );
}
@ -995,7 +1007,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
| ~context~ | ~qmckl_context~ | in | Global state |
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates |
| ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron distances |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
@ -1006,10 +1018,11 @@ integer function qmckl_compute_ee_distance_f(context, elec_num, walk_num, coord,
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(in) :: coord(elec_num,walk_num,3)
double precision , intent(out) :: ee_distance(elec_num,elec_num,walk_num)
integer*8 :: k
integer*8 :: k, i, j
double precision :: x, y, z
info = QMCKL_SUCCESS
@ -1030,8 +1043,8 @@ integer function qmckl_compute_ee_distance_f(context, elec_num, walk_num, coord,
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, &
coord(1,k,1), elec_num * walk_num, &
coord(1,k,1), elec_num * walk_num, &
ee_distance(1,1,k), elec_num)
if (info /= QMCKL_SUCCESS) then
exit
@ -1274,8 +1287,8 @@ integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale
do k=1,walk_num
info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num, &
coord(1,1,k), elec_num, &
coord(1,1,k), elec_num, &
coord(1,k,1), elec_num * walk_num, &
coord(1,k,1), elec_num * walk_num, &
ee_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_ee)
if (info /= QMCKL_SUCCESS) then
exit
@ -1898,7 +1911,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->nucleus.coord.data,
ctx->electron.en_distance);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1938,7 +1951,7 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n
integer*8 , intent(in) :: elec_num
integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
double precision , intent(in) :: elec_coord(elec_num,walk_num,3)
double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(out) :: en_distance(elec_num,nucl_num,walk_num)
@ -1968,7 +1981,7 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n
do k=1,walk_num
info = qmckl_distance(context, 'T', 'T', elec_num, nucl_num, &
elec_coord(1,1,k), elec_num, &
elec_coord(1,k,1), elec_num * walk_num, &
nucl_coord, nucl_num, &
en_distance(1,1,k), elec_num)
if (info /= QMCKL_SUCCESS) then
@ -2005,7 +2018,7 @@ qmckl_exit_code qmckl_compute_en_distance (
integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) , value :: walk_num
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3)
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
real (c_double ) , intent(out) :: en_distance(elec_num,nucl_num,walk_num)
@ -2052,10 +2065,10 @@ assert(qmckl_electron_provided(context));
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge (context, charge);
rc = qmckl_set_nucleus_charge (context, charge, nucl_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
@ -2103,6 +2116,7 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled)
{
@ -2174,7 +2188,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->nucleus.coord.data,
ctx->electron.en_distance_rescaled);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -2217,7 +2231,7 @@ integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_nu
integer*8 , intent(in) :: nucl_num
double precision , intent(in) :: rescale_factor_kappa_en
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
double precision , intent(in) :: elec_coord(elec_num,walk_num,3)
double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num)
@ -2253,7 +2267,7 @@ integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_nu
do k=1,walk_num
info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, nucl_num, &
elec_coord(1,1,k), elec_num, &
elec_coord(1,k,1), elec_num*walk_num, &
nucl_coord, nucl_num, &
en_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_en)
if (info /= QMCKL_SUCCESS) then
@ -2292,7 +2306,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled (
integer (c_int64_t) , intent(in) , value :: nucl_num
real (c_double ) , intent(in) , value :: rescale_factor_kappa_en
integer (c_int64_t) , intent(in) , value :: walk_num
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3)
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
real (c_double ) , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num)
@ -2341,10 +2355,10 @@ assert(qmckl_electron_provided(context));
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge (context, charge);
rc = qmckl_set_nucleus_charge (context, charge, nucl_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
@ -2360,6 +2374,7 @@ assert (rc == QMCKL_SUCCESS);
assert(fabs(en_distance_rescaled[0][0][0] - 0.9994721712909764) < 1.e-12);
// (1,2,1)
printf("%f\n%f\n", en_distance_rescaled[0][1][0] , 0.9998448354439821);
assert(fabs(en_distance_rescaled[0][1][0] - 0.9998448354439821) < 1.e-12);
// (2,1,1)
@ -2462,7 +2477,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->nucleus.coord.data,
ctx->electron.en_distance_rescaled_deriv_e);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -2506,7 +2521,7 @@ integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num,
integer*8 , intent(in) :: nucl_num
double precision , intent(in) :: rescale_factor_kappa_en
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
double precision , intent(in) :: elec_coord(elec_num,walk_num,3)
double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(out) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,walk_num)
@ -2542,7 +2557,7 @@ integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num,
do k=1,walk_num
info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, nucl_num, &
elec_coord(1,1,k), elec_num, &
elec_coord(1,k,1), elec_num*walk_num, &
nucl_coord, nucl_num, &
en_distance_rescaled_deriv_e(1,1,k), elec_num, rescale_factor_kappa_en)
if (info /= QMCKL_SUCCESS) then
@ -2581,7 +2596,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled_deriv_e (
integer (c_int64_t) , intent(in) , value :: nucl_num
real (c_double ) , intent(in) , value :: rescale_factor_kappa_en
integer (c_int64_t) , intent(in) , value :: walk_num
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3)
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
real (c_double ) , intent(out) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,walk_num)
@ -2611,10 +2626,10 @@ assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge (context, charge);
rc = qmckl_set_nucleus_charge (context, charge, nucl_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord);
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
@ -2732,7 +2747,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->nucleus.charge,
ctx->nucleus.charge.data,
ctx->electron.en_distance,
ctx->electron.en_pot);
if (rc != QMCKL_SUCCESS) {

60
org/qmckl_jastrow.org
View File

@ -1,4 +1,5 @@
#+TITLE: Jastrow Factor
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
@ -103,26 +104,25 @@ int main() {
computed data:
| Variable | Type | In/Out | Description |
|-------------------------------+---------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------|
| ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | |
| ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | |
| ~asymp_jasb~ | ~double[2]~ | Asymptotic component | |
| ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | |
| ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | |
| ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | |
| ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | |
| ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | |
| ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | vector of non-zero coefficients |
| ~dtmp_c~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num][0:cord_num-1]~ | vector of non-zero coefficients | |
| ~een_rescaled_e~ | ~double[walk_num][walk_num][elec_num][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | |
| ~een_rescaled_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | |
| ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| Variable | Type | In/Out | Description |
|----------------------------+---------------------------------------------------------------------+-------------------------------------------------+---------------------------------|
| ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | |
| ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | |
| ~asymp_jasb~ | ~double[2]~ | Asymptotic component | |
| ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | |
| ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | |
| ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | |
| ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | |
| ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | |
| ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | vector of non-zero coefficients |
| ~dtmp_c~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num][0:cord_num-1]~ | vector of non-zero coefficients | |
| ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | |
| ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
| ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | |
| ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | |
For H2O we have the following data:
@ -1093,7 +1093,7 @@ assert(rc == QMCKL_SUCCESS);
double elec_coord2[walk_num*3*elec_num];
rc = qmckl_get_electron_coord (context, 'N', elec_coord2);
rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*elec_num ; ++i) {
assert( elec_coord[i] == elec_coord2[i] );
@ -1131,15 +1131,15 @@ assert(k == nucl_rescale_factor_kappa);
double nucl_coord2[3*nucl_num];
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, nucl_num*3);
assert(rc == QMCKL_SUCCESS);
for (int64_t k=0 ; k<3 ; ++k) {
for (int64_t i=0 ; i<nucl_num ; ++i) {
@ -1147,7 +1147,7 @@ for (int64_t k=0 ; k<3 ; ++k) {
}
}
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, nucl_num*3);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*nucl_num ; ++i) {
assert( nucl_coord[i] == nucl_coord2[i] );
@ -1155,13 +1155,13 @@ for (int64_t i=0 ; i<3*nucl_num ; ++i) {
double nucl_charge2[nucl_num];
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<nucl_num ; ++i) {
assert( nucl_charge[i] == nucl_charge2[i] );
@ -1169,7 +1169,7 @@ for (int64_t i=0 ; i<nucl_num ; ++i) {
assert(qmckl_nucleus_provided(context));
#+end_src
* Computation
The computed data is stored in the context so that it can be reused
@ -3925,7 +3925,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
ctx->jastrow.cord_num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.coord_new,
ctx->nucleus.coord,
ctx->nucleus.coord.data,
ctx->electron.en_distance,
ctx->jastrow.een_rescaled_n,
ctx->jastrow.een_rescaled_n_deriv_e);

4
org/qmckl_local_energy.org
View File

@ -576,10 +576,10 @@ assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));

4
org/qmckl_mo.org
View File

@ -720,10 +720,10 @@ assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));

534
org/qmckl_nucleus.org
View File

@ -67,42 +67,42 @@ int main() {
The following data stored in the context:
|------------------------+----------------+-------------------------------------------|
| ~uninitialized~ | int32_t | Keeps bit set for uninitialized data |
| ~num~ | int64_t | Total number of nuclei |
| ~provided~ | bool | If true, ~nucleus~ is valid |
| ~charge~ | double[num] | Nuclear charges |
| ~coord~ | double[3][num] | Nuclear coordinates, in transposed format |
| ~coord_date~ | int64_t | Nuclear coordinates, date if modified |
| ~rescale_factor_kappa~ | double | The distance scaling factor |
|------------------------+--------------+-------------------------------------------|
| ~uninitialized~ | int32_t | Keeps bit set for uninitialized data |
| ~num~ | int64_t | Total number of nuclei |
| ~provided~ | bool | If true, ~nucleus~ is valid |
| ~charge~ | qmckl_vector | Nuclear charges |
| ~coord~ | qmckl_matrix | Nuclear coordinates, in transposed format |
| ~coord_date~ | int64_t | Nuclear coordinates, date if modified |
| ~rescale_factor_kappa~ | double | The distance scaling factor |
Computed data:
|-----------------------------+------------------+------------------------------------------------------------|
| ~nn_distance~ | double[num][num] | Nucleus-nucleus distances |
| ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed |
| ~nn_distance_rescaled~ | double[num][num] | Nucleus-nucleus rescaled distances |
| ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed |
| ~repulsion~ | double | Nuclear repulsion energy |
| ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed |
|-----------------------------+--------------+------------------------------------------------------------|
| ~nn_distance~ | qmckl_matrix | Nucleus-nucleus distances |
| ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed |
| ~nn_distance_rescaled~ | qmckl_matrix | Nucleus-nucleus rescaled distances |
| ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed |
| ~repulsion~ | double | Nuclear repulsion energy |
| ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed |
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_nucleus_struct {
int64_t num;
int64_t repulsion_date;
int64_t nn_distance_date;
int64_t nn_distance_rescaled_date;
int64_t coord_date;
double* coord;
double* charge;
double* nn_distance;
double* nn_distance_rescaled;
double repulsion;
double rescale_factor_kappa;
int32_t uninitialized;
bool provided;
int64_t num;
int64_t repulsion_date;
int64_t nn_distance_date;
int64_t nn_distance_rescaled_date;
int64_t coord_date;
qmckl_vector charge;
qmckl_matrix coord;
qmckl_matrix nn_distance;
qmckl_matrix nn_distance_rescaled;
double repulsion;
double rescale_factor_kappa;
int32_t uninitialized;
bool provided;
} qmckl_nucleus_struct;
#+end_src
@ -138,15 +138,6 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) {
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num);
qmckl_exit_code qmckl_get_nucleus_charge (const qmckl_context context, double* const charge);
qmckl_exit_code qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord);
qmckl_exit_code qmckl_get_nucleus_rescale_factor (const qmckl_context context, double* const rescale_factor_kappa);
#+end_src
#+NAME:post
#+begin_src c :exports none
if ( (ctx->nucleus.uninitialized & mask) != 0) {
@ -154,6 +145,13 @@ if ( (ctx->nucleus.uninitialized & mask) != 0) {
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_num(const qmckl_context context,
int64_t* const num);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
@ -187,10 +185,22 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) {
qmckl_get_nucleus_charge(const qmckl_context context,
double* const charge,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_charge (const qmckl_context context,
double* const charge,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
@ -215,16 +225,31 @@ qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) {
"nucleus data is not provided");
}
assert (ctx->nucleus.charge != NULL);
assert (ctx->nucleus.charge.data != NULL);
int64_t nucl_num = ctx->nucleus.num;
if (ctx->nucleus.num > size_max) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_nucleus_charge",
"Array too small");
}
memcpy(charge, ctx->nucleus.charge, nucl_num*sizeof(double));
qmckl_exit_code rc;
rc = qmckl_double_of_vector(context, ctx->nucleus.charge, charge, size_max);
return QMCKL_SUCCESS;
return rc;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_rescale_factor(const qmckl_context context,
double* const rescale_factor_kappa);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_rescale_factor (const qmckl_context context,
double* const rescale_factor_kappa)
@ -250,10 +275,24 @@ qmckl_get_nucleus_rescale_factor (const qmckl_context context,
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord) {
qmckl_get_nucleus_coord(const qmckl_context context,
const char transp,
double* const coord,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_coord (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;
@ -285,25 +324,21 @@ qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double*
"nucleus data is not provided");
}
int64_t nucl_num = ctx->nucleus.num;
assert (ctx->nucleus.coord.data != NULL);
assert (ctx->nucleus.coord != NULL);
qmckl_exit_code rc;
if (transp == 'N') {
qmckl_exit_code rc;
rc = qmckl_transpose(context, nucl_num, 3,
ctx->nucleus.coord, nucl_num,
coord, 3);
qmckl_matrix At = qmckl_matrix_alloc(context, 3, ctx->nucleus.coord.size[0]);
rc = qmckl_transpose(context, ctx->nucleus.coord, At);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_double_of_matrix(context, At, coord, size_max);
qmckl_matrix_free(context, At);
} else {
memcpy(coord, ctx->nucleus.coord, 3*nucl_num*sizeof(double));
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, size_max);
}
return QMCKL_SUCCESS;
return rc;
}
#+end_src
@ -330,51 +365,6 @@ bool qmckl_nucleus_provided(const qmckl_context context) {
** Initialization functions
To set the data relative to the nuclei in the context, the
following functions need to be called.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_nucleus_num (qmckl_context context, const int64_t num);
qmckl_exit_code qmckl_set_nucleus_charge (qmckl_context context, const double* charge);
qmckl_exit_code qmckl_set_nucleus_coord (qmckl_context context, const char transp, const double* coord);
qmckl_exit_code qmckl_set_nucleus_rescale_factor (qmckl_context context, const double rescale_factor_kappa);
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_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(in) , value :: num
end function
end interface
interface
integer(c_int32_t) function qmckl_set_nucleus_charge(context, charge) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) :: charge(*)
end function
end interface
interface
integer(c_int32_t) function qmckl_set_nucleus_coord(context, transp, coord) &
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(in) :: coord(*)
end function
end interface
#+end_src
#+NAME:pre2
#+begin_src c :exports none
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -392,12 +382,22 @@ ctx->nucleus.provided = (ctx->nucleus.uninitialized == 0);
return QMCKL_SUCCESS;
#+end_src
To set the data relative to the nuclei in the context, the
following functions need to be called.
To set the number of nuclei, use
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_num(qmckl_context context,
const int64_t num);
#+end_src
Sets the number of nuclei.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_num(qmckl_context context, const int64_t num) {
qmckl_set_nucleus_num(qmckl_context context,
const int64_t num)
{
<<pre2>>
if (num <= 0) {
@ -415,11 +415,35 @@ qmckl_set_nucleus_num(qmckl_context context, const int64_t num) {
}
#+end_src
The following function sets the nuclear charges of all the atoms.
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_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(in) , value :: num
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_charge(qmckl_context context,
const double* charge,
const int64_t size_max);
#+end_src
Sets the nuclear charges of all the atoms.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_charge(qmckl_context context, const double* charge) {
qmckl_set_nucleus_charge(qmckl_context context,
const double* charge,
const int64_t size_max)
{
<<pre2>>
if (charge == NULL) {
@ -437,34 +461,126 @@ qmckl_set_nucleus_charge(qmckl_context context, const double* charge) {
rc = qmckl_get_nucleus_num(context, &num);
if (rc != QMCKL_SUCCESS) return rc;
if (ctx->nucleus.charge != NULL) {
qmckl_free(context, ctx->nucleus.charge);
ctx->nucleus.charge= NULL;
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = num*sizeof(double);
assert (ctx->nucleus.charge == NULL);
ctx->nucleus.charge = (double*) qmckl_malloc(context, mem_info);
if (ctx->nucleus.charge == NULL) {
if (num > size_max) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
QMCKL_INVALID_ARG_3,
"qmckl_set_nucleus_charge",
NULL);
"Array too small");
}
ctx->nucleus.charge= memcpy(ctx->nucleus.charge, charge, num*sizeof(double));
assert (ctx->nucleus.charge != NULL);
ctx->nucleus.charge = qmckl_vector_alloc(context, num);
rc = qmckl_vector_of_double(context, charge, num, &(ctx->nucleus.charge));
<<post2>>
}
#+end_src
The following function sets the rescale parameter for the nuclear distances.
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_charge(context, charge, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) :: charge(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_coord(qmckl_context context,
const char transp,
const double* coord,
const int64_t size_max);
#+end_src
Sets the nuclear coordinates of all the atoms. The coordinates
are be given in atomic units.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context, const double rescale_factor_kappa) {
qmckl_set_nucleus_coord(qmckl_context context,
const char transp,
const double* coord,
const int64_t size_max)
{
<<pre2>>
qmckl_exit_code rc;
int32_t mask = 1 << 2;
const int64_t nucl_num = (int64_t) ctx->nucleus.num;
if (ctx->nucleus.coord.data != NULL) {
rc = qmckl_matrix_free(context, ctx->nucleus.coord);
if (rc != QMCKL_SUCCESS) return rc;
}
ctx->nucleus.coord = qmckl_matrix_alloc(context, nucl_num, 3);
if (ctx->nucleus.coord.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_nucleus_coord",
NULL);
}
if (size_max < 3*nucl_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_set_nucleus_coord",
"Array too small");
}
if (transp == 'N') {
qmckl_matrix At;
At = qmckl_matrix_alloc(context, 3, nucl_num);
rc = qmckl_matrix_of_double(context, coord, 3*nucl_num, &At);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_transpose(context, At, ctx->nucleus.coord);
} else {
rc = qmckl_matrix_of_double(context, coord, nucl_num*3, &(ctx->nucleus.coord));
}
if (rc != QMCKL_SUCCESS) return rc;
<<post2>>
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_coord(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(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context,
const double kappa);
#+end_src
Sets the rescale parameter for the nuclear distances.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context,
const double rescale_factor_kappa)
{
<<pre2>>
if (rescale_factor_kappa <= 0.0) {
@ -480,50 +596,17 @@ qmckl_set_nucleus_rescale_factor(qmckl_context context, const double rescale_fac
}
#+end_src
The following function sets the nuclear coordinates of all the
atoms. The coordinates should be given in atomic units.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_nucleus_coord(qmckl_context context, const char transp, const double* coord) {
<<pre2>>
int64_t nucl_num = (int64_t) 0;
qmckl_exit_code rc;
int32_t mask = 1 << 2;
rc = qmckl_get_nucleus_num(context, &nucl_num);
if (rc != QMCKL_SUCCESS) return rc;
if (ctx->nucleus.coord != NULL) {
qmckl_free(context, ctx->nucleus.coord);
ctx->nucleus.coord = NULL;
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3*nucl_num*sizeof(double);
assert(ctx->nucleus.coord == NULL);
ctx->nucleus.coord = (double*) qmckl_malloc(context, mem_info);
if (coord == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_nucleus_coord",
NULL);
}
if (transp == 'N') {
rc = qmckl_transpose(context, 3, nucl_num,
coord, 3,
ctx->nucleus.coord, nucl_num);
if (rc != QMCKL_SUCCESS) return rc;
} else {
memcpy(ctx->nucleus.coord, coord, 3*nucl_num*sizeof(double));
}
<<post2>>
}
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_rescale_factor(context, kappa) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) , value :: kappa
end function
end interface
#+end_src
** Test
@ -568,15 +651,15 @@ assert(k == nucl_rescale_factor_kappa);
double nucl_coord2[3*nucl_num];
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
for (size_t k=0 ; k<3 ; ++k) {
for (size_t i=0 ; i<nucl_num ; ++i) {
@ -584,7 +667,7 @@ for (size_t k=0 ; k<3 ; ++k) {
}
}
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
for (size_t i=0 ; i<3*nucl_num ; ++i) {
assert( nucl_coord[i] == nucl_coord2[i] );
@ -592,13 +675,13 @@ for (size_t i=0 ; i<3*nucl_num ; ++i) {
double nucl_charge2[nucl_num];
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_SUCCESS);
for (size_t i=0 ; i<nucl_num ; ++i) {
assert( nucl_charge[i] == nucl_charge2[i] );
@ -621,11 +704,17 @@ assert(qmckl_nucleus_provided(context));
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_nucleus_nn_distance(qmckl_context context, double* distance);
qmckl_exit_code
qmckl_get_nucleus_nn_distance(qmckl_context context,
double* distance,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_nucleus_nn_distance(qmckl_context context, double* distance)
qmckl_exit_code
qmckl_get_nucleus_nn_distance(qmckl_context context,
double* distance,
const int64_t size_max)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -638,22 +727,29 @@ qmckl_exit_code qmckl_get_nucleus_nn_distance(qmckl_context context, double* dis
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = ctx->nucleus.num * ctx->nucleus.num;
memcpy(distance, ctx->nucleus.nn_distance, sze * sizeof(double));
const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
if (sze > size_max) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_3,
"qmckl_get_nucleus_nn_distance",
"Array too small");
}
rc = qmckl_double_of_matrix(context, ctx->nucleus.nn_distance, distance, size_max);
return QMCKL_SUCCESS;
return rc;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_nucleus_nn_distance(context, distance) &
integer(c_int32_t) function qmckl_get_nucleus_nn_distance(context, distance, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(out) :: distance(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
@ -678,26 +774,23 @@ qmckl_exit_code qmckl_provide_nn_distance(qmckl_context context)
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
/* Allocate array */
if (ctx->nucleus.nn_distance == NULL) {
if (ctx->nucleus.nn_distance.data == NULL) {
ctx->nucleus.nn_distance =
qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->nucleus.num * ctx->nucleus.num * sizeof(double);
double* nn_distance = (double*) qmckl_malloc(context, mem_info);
if (nn_distance == NULL) {
if (ctx->nucleus.nn_distance.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_nn_distance",
NULL);
}
ctx->nucleus.nn_distance = nn_distance;
}
qmckl_exit_code rc =
qmckl_compute_nn_distance(context,
ctx->nucleus.num,
ctx->nucleus.coord,
ctx->nucleus.nn_distance);
ctx->nucleus.coord.data,
ctx->nucleus.nn_distance.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -788,7 +881,7 @@ qmckl_exit_code qmckl_compute_nn_distance (
assert(qmckl_nucleus_provided(context));
double distance[nucl_num*nucl_num];
rc = qmckl_get_nucleus_nn_distance(context, distance);
rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num);
assert(distance[0] == 0.);
assert(distance[1] == distance[nucl_num]);
assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
@ -800,11 +893,17 @@ assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled);
qmckl_exit_code
qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
double* distance_rescaled,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled)
qmckl_exit_code
qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
double* distance_rescaled,
const int64_t size_max)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -817,13 +916,35 @@ qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, do
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = ctx->nucleus.num * ctx->nucleus.num;
memcpy(distance_rescaled, ctx->nucleus.nn_distance_rescaled, sze * sizeof(double));
return QMCKL_SUCCESS;
const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
if (sze > size_max) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_3,
"qmckl_get_nucleus_nn_distance_rescaled",
"Array too small");
}
rc = qmckl_double_of_matrix(context,
ctx->nucleus.nn_distance_rescaled,
distance_rescaled,
size_max);
return rc;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_nucleus_nn_distance_rescaled(context, distance_rescaled, size_max) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(out) :: distance_rescaled(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
@ -844,27 +965,24 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context)
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
/* Allocate array */
if (ctx->nucleus.nn_distance_rescaled == NULL) {
if (ctx->nucleus.nn_distance_rescaled.data == NULL) {
ctx->nucleus.nn_distance_rescaled =
qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->nucleus.num * ctx->nucleus.num * sizeof(double);
double* nn_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
if (nn_distance_rescaled == NULL) {
if (ctx->nucleus.nn_distance_rescaled.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_nn_distance_rescaled",
NULL);
}
ctx->nucleus.nn_distance_rescaled = nn_distance_rescaled;
}
qmckl_exit_code rc =
qmckl_compute_nn_distance_rescaled(context,
ctx->nucleus.num,
ctx->nucleus.rescale_factor_kappa,
ctx->nucleus.coord,
ctx->nucleus.nn_distance_rescaled);
ctx->nucleus.coord.data,
ctx->nucleus.nn_distance_rescaled.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -959,7 +1077,7 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled (
//assert(qmckl_nucleus_provided(context));
//
//double distance[nucl_num*nucl_num];
//rc = qmckl_get_nucleus_nn_distance(context, distance);
//rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num);
//assert(distance[0] == 0.);
//assert(distance[1] == distance[nucl_num]);
//assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
@ -975,11 +1093,11 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled (
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* energy);
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const energy);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* energy)
qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const energy)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@ -1037,8 +1155,8 @@ qmckl_exit_code qmckl_provide_nucleus_repulsion(qmckl_context context)
rc = qmckl_compute_nucleus_repulsion(context,
ctx->nucleus.num,
ctx->nucleus.charge,
ctx->nucleus.nn_distance,
ctx->nucleus.charge.data,
ctx->nucleus.nn_distance.data,
&(ctx->nucleus.repulsion));
if (rc != QMCKL_SUCCESS) {
return rc;

357
org/qmckl_point.org
View File

@ -17,11 +17,14 @@ walkers.
#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
@ -33,6 +36,8 @@ walkers.
#endif
#include "chbrclf.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
int main() {
qmckl_context context;
@ -61,36 +66,31 @@ int main() {
#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_memory_private_func.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 |
| ~coord_x~ | ~double[num]~ | X coordinates |
| ~coord_y~ | ~double[num]~ | Y coordinates |
| ~coord_z~ | ~double[num]~ | Z coordinates |
| Variable | Type | Description |
|----------+----------------+------------------------|
| ~num~ | ~int64_t~ | Total number of points |
| ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 |
We consider that 'transposed' and 'normal' storage follows the convention:
| | Normal | Transposed |
|---------+------------------+------------------|
| C | ~[point_num][3]~ | ~[3][point_num]~ |
| Fortran | ~(3,point_num)~ | ~(point_num,3)~ |
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 {
double* coord_x;
double* coord_y;
double* coord_z;
int64_t num;
int64_t num;
qmckl_matrix coord;
} qmckl_point_struct;
#+end_src
@ -182,6 +182,7 @@ end interface
#+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
@ -193,6 +194,7 @@ qmckl_exit_code qmckl_get_point(const qmckl_context context,
#+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)
{
@ -213,10 +215,9 @@ qmckl_get_point(const qmckl_context context,
assert (ctx->point != NULL);
int64_t point_num = ctx->point->num;
if (point_num == 0) return QMCKL_NOT_PROVIDED;
assert (ctx->point->coord_x != NULL);
assert (ctx->point->coord_y != NULL);
assert (ctx->point->coord_z != NULL);
assert (ctx->point->coord.data != NULL);
if (size_max < 3*point_num) {
return qmckl_failwith( context,
@ -225,27 +226,33 @@ qmckl_get_point(const qmckl_context context,
"size_max too small");
}
double * ptr = coord;
for (int64_t i=0 ; i<point_num ; ++i) {
,*ptr = ctx->point->coord_x[i]; ++ptr;
,*ptr = ctx->point->coord_y[i]; ++ptr;
,*ptr = ctx->point->coord_z[i]; ++ptr;
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 QMCKL_SUCCESS;
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, coord, size_max) bind(C)
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
@ -253,165 +260,19 @@ end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_get_point_xyz (const qmckl_context context,
double* const coord_x,
double* const coord_y,
double* const coord_z,
const int64_t size_max);
#+end_src
Returns the point coordinates in three different arrays, one for
each component x,y,z.
The pointers are assumed to point on a memory block of size
~size_max~ \ge ~point_num~.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_point_xyz (const qmckl_context context,
double* const coord_x,
double* const coord_y,
double* const coord_z,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (coord_x == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_point_coord_xyz",
"coord_x is a null pointer");
}
if (coord_y == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_point_coord_xyz",
"coord_y is a null pointer");
}
if (coord_z == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_get_point_coord_xyz",
"coord_z is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
assert (ctx->point != NULL);
int64_t point_num = ctx->point->num;
assert (ctx->point->coord_x != NULL);
assert (ctx->point->coord_y != NULL);
assert (ctx->point->coord_z != NULL);
if (size_max < point_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_5,
"qmckl_get_point_coord_xyz",
"size_max too small");
}
memcpy(coord_x, ctx->point->coord_x, point_num*sizeof(double));
memcpy(coord_y, ctx->point->coord_y, point_num*sizeof(double));
memcpy(coord_z, ctx->point->coord_z, point_num*sizeof(double));
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_xyz(context, &
coord_x, coord_y, coord_z, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(out) :: coord_x(*)
real (c_double ) , intent(out) :: coord_y(*)
real (c_double ) , intent(out) :: coord_z(*)
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.
#+NAME: check_alloc
#+begin_src c :exports none
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
assert (ctx->point != NULL);
if (ctx->point->num < num) {
if (ctx->point->coord_x != NULL) {
qmckl_free(context, ctx->point->coord_x);
ctx->point->coord_x = NULL;
}
if (ctx->point->coord_y != NULL) {
qmckl_free(context, ctx->point->coord_y);
ctx->point->coord_y = NULL;
}
if (ctx->point->coord_z != NULL) {
qmckl_free(context, ctx->point->coord_z);
ctx->point->coord_z = NULL;
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = num*sizeof(double);
ctx->point->coord_x = (double*) qmckl_malloc(context, mem_info);
if (ctx->point->coord_x == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
NULL);
}
ctx->point->coord_y = (double*) qmckl_malloc(context, mem_info);
if (ctx->point->coord_y == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
NULL);
}
ctx->point->coord_z = (double*) qmckl_malloc(context, mem_info);
if (ctx->point->coord_z == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_point",
NULL);
}
};
ctx->point->num = num;
#+end_src
To set the data relative to the points in the context, one of the
following functions need to be called.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_point (qmckl_context context,
const char transp,
const double* coord,
const int64_t num);
const int64_t size_max);
#+end_src
Copy a sequence of (x,y,z) into the context.
@ -419,16 +280,61 @@ qmckl_exit_code qmckl_set_point (qmckl_context context,
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code
qmckl_set_point (qmckl_context context,
const char transp,
const double* coord,
const int64_t num)
{
<<check_alloc>>
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
for (int64_t i=0 ; i<num ; ++i) {
ctx->point->coord_x[i] = coord[3*i ];
ctx->point->coord_y[i] = coord[3*i+1];
ctx->point->coord_z[i] = coord[3*i+2];
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* const) context;
assert (ctx != NULL);
assert (ctx->point != NULL);
qmckl_exit_code rc;
if (ctx->point->num < num) {
if (ctx->point->coord.data != NULL) {
rc = qmckl_matrix_free(context, ctx->point->coord);
return rc;
}
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') {
memcpy(ctx->point->coord.data, coord, 3*num*sizeof(double));
} else {
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];
}
}
return QMCKL_SUCCESS;
@ -439,76 +345,38 @@ qmckl_set_point (qmckl_context context,
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_set_point(context, &
coord_x, coord_y, coord_z, size_max) bind(C)
transp, coord, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(in) :: coord_x(*)
real (c_double ) , intent(in) :: coord_y(*)
real (c_double ) , intent(in) :: coord_z(*)
character(c_char) , intent(in) , value :: transp
real (c_double ) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_point_xyz (qmckl_context context,
const double* coord_x,
const double* coord_y,
const double* coord_z,
const int64_t num);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code
qmckl_set_point_xyz (qmckl_context context,
const double* coord_x,
const double* coord_y,
const double* coord_z,
const int64_t num)
{
<<check_alloc>>
memcpy(ctx->point->coord_x, coord_x, num*sizeof(double));
memcpy(ctx->point->coord_y, coord_y, num*sizeof(double));
memcpy(ctx->point->coord_z, coord_z, num*sizeof(double));
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_xyz(context, &
coord_x, coord_y, coord_z, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double ) , intent(in) :: coord_x(*)
real (c_double ) , intent(in) :: coord_y(*)
real (c_double ) , intent(in) :: coord_z(*)
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;
double* coord = &(chbrclf_elec_coord[0][0][0]);
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_set_point (context, coord, point_num);
rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_point (context, 'N', coord, point_num);
assert(rc == QMCKL_SUCCESS);
int64_t n;
@ -516,25 +384,30 @@ rc = qmckl_get_point_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == point_num);
double coord2[point_num*3];
double coord_x[point_num];
double coord_y[point_num];
double coord_z[point_num];
rc = qmckl_get_point_xyz (context, coord_x, coord_y, coord_z, point_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_point (context, coord2, (point_num*3));
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] == coord_x[i] );
assert( coord[3*i+1] == coord_y[i] );
assert( coord[3*i+2] == coord_z[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', coord2, point_num);
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

10
org/qmckl_trexio.org
View File

@ -208,8 +208,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)
trexio_string_of_error(rcio));
}
//rc = qmckl_set_nucleus_charge(context, nucl_charge, nucleus_num);
rc = qmckl_set_nucleus_charge(context, nucl_charge);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucleus_num);
qmckl_free(context, nucl_charge);
nucl_charge = NULL;
@ -249,8 +248,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)
trexio_string_of_error(rcio));
}
//rc = qmckl_set_nucleus_coord(context, 'N', nucl_coord, 3*nucleus_num);
rc = qmckl_set_nucleus_coord(context, 'N', nucl_coord);
rc = qmckl_set_nucleus_coord(context, 'N', nucl_coord, 3*nucleus_num);
qmckl_free(context, nucl_coord);
nucl_coord = NULL;
@ -1183,7 +1181,7 @@ assert (nucl_num == chbrclf_nucl_num);
printf("Nuclear charges\n");
double * charge = (double*) malloc (nucl_num * sizeof(double));
rc = qmckl_get_nucleus_charge(context, charge);
rc = qmckl_get_nucleus_charge(context, charge, nucl_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i<nucl_num ; i++) {
assert (charge[i] == chbrclf_charge[i]);
@ -1193,7 +1191,7 @@ charge = NULL;
printf("Nuclear coordinates\n");
double * coord = (double*) malloc (nucl_num * 3 * sizeof(double));
rc = qmckl_get_nucleus_coord(context, 'T', coord);
rc = qmckl_get_nucleus_coord(context, 'T', coord, 3*nucl_num);
assert (rc == QMCKL_SUCCESS);
int k=0;
for (int j=0 ; j<3 ; ++j) {

229
org/qmckl_utils.org
View File

@ -1,229 +0,0 @@
#+TITLE: Utility functions
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
* 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"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
int main() {
qmckl_context context;
context = qmckl_context_create();
#+end_src
* Matrix operations
** ~qmckl_transpose~
Transposes a matrix: $B_{ji} = A_{ij}$
#+NAME: qmckl_transpose_args
| qmckl_context | context | in | Global state |
| int64_t | m | in | Number of rows of the input matrix |
| int64_t | n | in | Number of columns of the input matrix |
| double | A[][lda] | in | Array containing the $m \times n$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[][ldb] | out | Array containing the $n \times m$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
*** Requirements
- ~context~ is not ~QMCKL_NULL_CONTEXT~
- ~m > 0~
- ~n > 0~
- ~lda >= m~
- ~ldb >= n~
- ~A~ is allocated with at least $m \times n \times 8$ bytes
- ~B~ is allocated with at least $n \times m \times 8$ bytes
*** C header
#+CALL: generate_c_header(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_transpose (
const qmckl_context context,
const int64_t m,
const int64_t n,
const double* A,
const int64_t lda,
double* const B,
const int64_t ldb );
#+end_src
*** Source
#+begin_src f90 :tangle (eval f)
integer function qmckl_transpose_f(context, m, n, A, LDA, B, LDB) &
result(info)
use qmckl
implicit none
integer(qmckl_context) , intent(in) :: context
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(out) :: B(ldb,*)
integer*8 :: i,j
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (m <= 0_8) then
info = QMCKL_INVALID_ARG_2
return
endif
if (n <= 0_8) then
info = QMCKL_INVALID_ARG_3
return
endif
if (LDA < m) then
info = QMCKL_INVALID_ARG_5
return
endif
if (LDB < n) then
info = QMCKL_INVALID_ARG_7
return
endif
do j=1,m
do i=1,n
B(i,j) = A(j,i)
end do
end do
end function qmckl_transpose_f
#+end_src
*** C interface :noexport:
#+CALL: generate_c_interface(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_transpose &
(context, m, n, A, lda, B, ldb) &
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 :: 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(out) :: B(ldb,*)
integer (c_int64_t) , intent(in) , value :: ldb
integer(c_int32_t), external :: qmckl_transpose_f
info = qmckl_transpose_f &
(context, m, n, A, lda, B, ldb)
end function qmckl_transpose
#+end_src
#+CALL: generate_f_interface(table=qmckl_transpose_args,rettyp="qmckl_exit_code",fname="qmckl_transpose")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_transpose &
(context, m, n, A, lda, B, ldb) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
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(out) :: B(ldb,*)
integer (c_int64_t) , intent(in) , value :: ldb
end function qmckl_transpose
end interface
#+end_src
*** Test :noexport:
#+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_transpose(context) bind(C)
use qmckl
implicit none
integer(qmckl_context), intent(in), value :: context
double precision, allocatable :: A(:,:), B(:,:)
integer*8 :: m, n, LDA, LDB
integer*8 :: i,j
double precision :: x
m = 5
n = 6
LDA = m+3
LDB = n+1
allocate( A(LDA,n), B(LDB,m) )
A = 0.d0
B = 0.d0
do j=1,n
do i=1,m
A(i,j) = -10.d0 + dble(i+j)
end do
end do
test_qmckl_transpose = qmckl_transpose(context, m, n, A, LDA, B, LDB)
if (test_qmckl_transpose /= QMCKL_SUCCESS) return
test_qmckl_transpose = QMCKL_FAILURE
x = 0.d0
do j=1,n
do i=1,m
x = x + (A(i,j)-B(j,i))**2
end do
end do
if (dabs(x) <= 1.d-15) then
test_qmckl_transpose = QMCKL_SUCCESS
endif
deallocate(A,B)
end function test_qmckl_transpose
#+end_src
#+begin_src c :comments link :tangle (eval c_test)
qmckl_exit_code test_qmckl_transpose(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_transpose(context));
#+end_src
* End of files :noexport:
#+begin_src c :comments link :tangle (eval c_test)
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0;
}
#+end_src
# -*- mode: org -*-
# vim: syntax=c