TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-22 20:36:01 +01:00
Code Issues Projects Releases Wiki Activity

Merge branch 'jastrow_c' into qmckl_compute_dim_cord_vect_f

Browse Source
This commit is contained in:
Gianfranco Abrusci 2022-03-09 11:23:43 +01:00
commit 712f595ed6
9 changed files with 687 additions and 233 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
.github/workflows/test-build.yml vendored
View File

@ -47,6 +47,7 @@ jobs:
- name: Run test - name: Run test
run: make -j 4 check run: make -j 4 check
working-directory: _build
- name: Archive test log file - name: Archive test log file
if: failure() if: failure()
@ -57,13 +58,7 @@ jobs:
- name: Dist test - name: Dist test
run: make distcheck run: make distcheck
working-directory: _build
- name: Archive test log file
if: failure()
uses: actions/upload-artifact@v2
with:
name: dist-report-ubuntu
path: test-suite.log
x86_macos: x86_macos:
@ -91,37 +86,26 @@ jobs:
git clone https://github.com/TREX-CoE/trexio.git git clone https://github.com/TREX-CoE/trexio.git
cd trexio cd trexio
./autogen.sh ./autogen.sh
mkdir _build
cd _build
../configure --enable-silent-rules
make -j 4
- name: Run test
run: make -j 4 check
- name: Archive test log file
if: failure()
uses: actions/upload-artifact@v2
with:
name: test-report-ubuntu
path: test-suite.log
- name: Dist test
run: make distcheck
- name: Archive test log file
if: failure()
uses: actions/upload-artifact@v2
with:
./configure --prefix=${PWD}/_install --enable-silent-rules ./configure --prefix=${PWD}/_install --enable-silent-rules
make -j 4 make -j 4
make install make install
- name: Test TREXIO
run: make -j 4 check
working-directory: trexio
- name: Archive TREXIO test log file
if: failure()
uses: actions/upload-artifact@v2
with:
name: test-report-trexio-macos
path: trexio/test-suite.log
- name: Build QMCkl - name: Build QMCkl
run: | run: |
export PKG_CONFIG_PATH=${PWD}/trexio/_install/lib/pkgconfig:$PKG_CONFIG_PATH export PKG_CONFIG_PATH=${PWD}/trexio/_install/lib/pkgconfig:$PKG_CONFIG_PATH
./autogen.sh ./autogen.sh
./configure --enable-silent-rules --enable-debug ./configure CC=gcc-10 FC=gfortran-10 --enable-silent-rules --enable-debug
make -j 4 make -j 4
- name: Run test - name: Run test

18
README.md
View File

@ -73,6 +73,24 @@ sudo make install
sudo make installcheck sudo make installcheck
``` ```
## Linking to your program
The `make install` command takes care of installing the QMCkl shared library on the user machine.
Once installed, add `-lqmckl` to the list of compiler options.
In some cases (e.g. when using custom `prefix` during configuration), the QMCkl library might end up installed in a directory, which is absent in the default `$LIBRARY_PATH`.
In order to link the program against QMCkl, the search paths can be modified as follows:
`export LIBRARY_PATH=$LIBRARY_PATH:<path_to_qmckl>/lib`
(same holds for `$LD_LIBRARY_PATH`). The `<path_to_qmckl>` has to be replaced with the prefix used during the installation.
If your project relies on the CMake build system, feel free to use the
[FindQMCKL.cmake](https://github.com/TREX-CoE/qmckl/blob/master/cmake/FindQMCKL.cmake)
module to find and link the QMCkl library automatically.
## Verificarlo CI ## Verificarlo CI
Since Verificarlo should not be a dependency of QMCkl, all Verificarlo Since Verificarlo should not be a dependency of QMCkl, all Verificarlo

77
cmake/FindQMCKL.cmake Normal file
View File

@ -0,0 +1,77 @@
#===========================================
# Try to find the QMCkl library;
# If found, it will define the following variables (note the plural form):
# QMCKL_FOUND - System has libqmckl;
# QMCKL_INCLUDE_DIRS - The QMCKL include directories;
# QMCKL_LIBRARIES - The libraries needed to use QMCKL;
# If QMCKL has been installed in a non-standard location, one can set an
# environment variable $QMCKL_DIR in the current shell:
# $ export QMCKL_DIR=<custom_path>
# to indicate the prefix used during the QMCKL installation
# (typically `./configure prefix=<custom_path> ..` or `cmake -DCMAKE_INSTALL_DIR=<custom_path> ..`)
# This file should be located WITHIN your project source tree.
# (e.g. in cmake/FindQMCKL.cmake)
# How to use it in your project CMakeLists.txt:
# This is needed to locate FindQMCKL.cmake file, modify it according to your source tree.
# list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake/")
# find_package(QMCKL)
# if (QMCKL_FOUND)
# include_directories(${QMCKL_INCLUDE_DIRS})
# target_link_libraries(your_target ${QMCKL_LIBRARIES})
# endif()
#===========================================
# This file is distirbuted under the BSD 3-Clause License.
# Copyright (c) 2021, TREX Center of Excellence
#===========================================
message("<FindQMCKL.cmake>")
set(QMCKL_SEARCH_PATHS
~/Library/Frameworks
/Library/Frameworks
/usr/local
/usr
/sw # Fink
/opt/local # DarwinPorts
/opt/csw # Blastwave
/opt
)
find_path(QMCKL_INCLUDE_DIR
NAMES qmckl.h
HINTS $ENV{QMCKL_DIR}
PATH_SUFFIXES include/qmckl include
PATHS ${QMCKL_SEARCH_PATHS}
)
# No need to specify platform-specific prefix (e.g. libqmckl on Unix) or
# suffix (e.g. .so on Unix or .dylib on MacOS) in NAMES. CMake takes care of that.
find_library(QMCKL_LIBRARY
NAMES qmckl
HINTS $ENV{QMCKL_DIR}
PATH_SUFFIXES lib64 lib
PATHS ${QMCKL_SEARCH_PATHS}
)
message("<FindQMCKL.cmake>")
# Handle the QUIETLY and REQUIRED arguments and set QMCKL_FOUND to TRUE if
# all listed variables are TRUE.
INCLUDE(FindPackageHandleStandardArgs)
FIND_PACKAGE_HANDLE_STANDARD_ARGS(QMCKL DEFAULT_MSG QMCKL_LIBRARY QMCKL_INCLUDE_DIR )
MARK_AS_ADVANCED(QMCKL_INCLUDE_DIR QMCKL_LIBRARY)
# Mot setting _INCLUDE_DIR and _LIBRARIES is considered a bug,
# see https://gitlab.kitware.com/cmake/community/-/wikis/doc/tutorials/How-To-Find-Libraries
set(QMCKL_LIBRARIES ${QMCKL_LIBRARY})
set(QMCKL_INCLUDE_DIRS ${QMCKL_INCLUDE_DIR})

4
configure.ac
View File

@ -49,12 +49,12 @@ AS_IF([test -d ${srcdir}/.git],
AC_ARG_WITH(ifort, [AS_HELP_STRING([--with-ifort],[Use Intel Fortran compiler])], with_ifort=$withval, with_ifort=no) AC_ARG_WITH(ifort, [AS_HELP_STRING([--with-ifort],[Use Intel Fortran compiler])], with_ifort=$withval, with_ifort=no)
AS_IF([test "$with_ifort" == "yes"], [ AS_IF([test "$with_ifort" == "yes"], [
FC=ifort FC=ifort
FCFLAGS="-xHost -ip -O2 -ftz -finline -g -mkl=sequential" ]) FCFLAGS="-xHost -ip -Ofast -ftz -finline -g -mkl=sequential" ])
AC_ARG_WITH(icc, [AS_HELP_STRING([--with-icc],[Use Intel C compiler])], with_icc=$withval, with_icc=no) AC_ARG_WITH(icc, [AS_HELP_STRING([--with-icc],[Use Intel C compiler])], with_icc=$withval, with_icc=no)
AS_IF([test "$with_icc" == "yes"], [ AS_IF([test "$with_icc" == "yes"], [
CC=icc CC=icc
CFLAGS="-xHost -ip -O2 -ftz -finline -g -mkl=sequential" ]) CFLAGS="-xHost -ip -Ofast -ftz -finline -g -mkl=sequential" ])
AS_IF([test "$with_icc"."$with_ifort" == "yes.yes"], [ AS_IF([test "$with_icc"."$with_ifort" == "yes.yes"], [
ax_blas_ok="yes" ax_blas_ok="yes"

250
org/qmckl_ao.org
View File

@ -1,3 +1,4 @@
B
#+TITLE: Atomic Orbitals #+TITLE: Atomic Orbitals
#+SETUPFILE: ../tools/theme.setup #+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org #+INCLUDE: ../tools/lib.org
@ -268,7 +269,7 @@ end interface
*** Data structure :noexport: *** Data structure :noexport:
#+begin_src c :comments org :tangle (eval h_private_type) #+begin_src c :comments org :tangle (eval h_private_type) :noweb yes
typedef struct qmckl_ao_basis_struct { typedef struct qmckl_ao_basis_struct {
int64_t shell_num; int64_t shell_num;
int64_t prim_num; int64_t prim_num;
@ -299,6 +300,9 @@ typedef struct qmckl_ao_basis_struct {
bool provided; bool provided;
bool ao_cartesian; bool ao_cartesian;
char type; char type;
#ifdef HAVE_HPC
<<HPC_struct>>
#endif
} qmckl_ao_basis_struct; } qmckl_ao_basis_struct;
#+end_src #+end_src
@ -2485,7 +2489,6 @@ free(ao_factor_test);
| ~ao_vgl~ | ~double[point_num][5][ao_num]~ | Value, gradients, Laplacian of the primitives at current positions | | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | Value, gradients, Laplacian of the primitives at current positions |
| ~ao_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions | | ~ao_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions |
*** After initialization *** After initialization
When the basis set is completely entered, extra data structures may be When the basis set is completely entered, extra data structures may be
@ -2625,9 +2628,112 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
} }
} }
} }
/* TODO : sort the basis set here */
rc = QMCKL_SUCCESS;
#ifdef HAVE_HPC
rc = qmckl_finalize_basis_hpc(context);
#endif
return rc;
}
#+end_src
*** HPC-specific data structures
For faster access, we provide extra arrays for the shell information as:
- $C_{psa}$ = =coef_per_nucleus (iprim, ishell, inucl)=
- $\gamma_{pa}$ =expo_per_nucleus (iprim, inucl)=
such that the computation of the radial parts can be vectorized.
Exponents are sorted in increasing order to exit loops quickly,
and only unique exponents are kept. This also allows to pack the
exponents to enable vectorization of exponentials.
The computation of the radial part is made as
\[
R_{sa} = \sum_p C_{psa} \gamma_{pa}
\]
which is a matrix-vector product.
#+NAME:HPC_struct
#+begin_src c :comments org :exports none
/* HPC specific data structures */
qmckl_tensor coef_per_nucleus;
qmckl_matrix expo_per_nucleus;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
#ifdef HAVE_HPC
qmckl_exit_code qmckl_finalize_basis_hpc (qmckl_context context);
#endif
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
#ifdef HAVE_HPC
qmckl_exit_code qmckl_finalize_basis_hpc (qmckl_context context)
{
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
int64_t shell_max = 0;
int64_t prim_max = 0;
int64_t nucl_num = ctx->nucleus.num;
for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
shell_max = ctx->ao_basis.nucleus_shell_num[inucl] > shell_max ?
ctx->ao_basis.nucleus_shell_num[inucl] : shell_max;
int64_t prim_num = 0;
const int64_t ishell_start = ctx->ao_basis.nucleus_index[inucl];
const int64_t ishell_end = ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl];
for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
prim_num += ctx->ao_basis.shell_prim_num[ishell];
}
prim_max = prim_num > prim_max ?
prim_num : prim_max;
}
int64_t size[3] = { prim_max, shell_max, nucl_num };
ctx->ao_basis.coef_per_nucleus = qmckl_tensor_alloc( context, 3, size );
ctx->ao_basis.coef_per_nucleus = qmckl_tensor_set(ctx->ao_basis.coef_per_nucleus, 0.);
ctx->ao_basis.expo_per_nucleus = qmckl_matrix_alloc( context, prim_max, nucl_num );
ctx->ao_basis.expo_per_nucleus = qmckl_matrix_set(ctx->ao_basis.expo_per_nucleus, 0.);
double expo[prim_max];
double coef[shell_max][prim_max];
for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
memset(expo, 0, sizeof(expo));
memset(coef, 0, sizeof(expo));
int64_t idx = 0;
const int64_t ishell_start = ctx->ao_basis.nucleus_index[inucl];
const int64_t ishell_end = ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl];
for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
const int64_t iprim_start = ctx->ao_basis.shell_prim_index[ishell];
const int64_t iprim_end = ctx->ao_basis.shell_prim_index[ishell] + ctx->ao_basis.shell_prim_num[ishell];
for (int64_t iprim = iprim_start ; iprim < iprim_end ; ++iprim) {
expo[idx] = ctx->ao_basis.coefficient_normalized[iprim];
coef[ishell - ishell_start][idx] = ctx->ao_basis.coefficient_normalized[iprim];
idx += 1;
}
}
for (int64_t i=0 ; i<prim_max ; ++i) {
qmckl_mat( ctx->ao_basis.expo_per_nucleus, i, inucl ) = expo[i];
}
for (int64_t j=0 ; j<shell_max ; ++j) {
for (int64_t i=0 ; i<prim_max ; ++i) {
qmckl_ten3( ctx->ao_basis.coef_per_nucleus, i, j, inucl ) = coef[j][i];
}
}
}
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
} }
#endif
#+end_src #+end_src
*** Access functions *** Access functions
@ -3040,11 +3146,11 @@ integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
integer(c_int64_t), intent(in), value :: context integer(c_int64_t), intent(in), value :: context
integer*8 :: n, ldv, j, i integer*8 :: n, ldv, j, i
double precision :: X(3), R(3), Y(3), r2 double precision :: X(3), R(3), Y(3), r2, z
double precision, allocatable :: VGL(:,:), A(:) double precision, allocatable :: VGL(:,:), A(:)
double precision :: epsilon double precision :: epsilon
epsilon = qmckl_get_numprec_epsilon(context) epsilon = 3.d0 * qmckl_get_numprec_epsilon(context)
X = (/ 1.1 , 2.2 , 3.3 /) X = (/ 1.1 , 2.2 , 3.3 /)
R = (/ 0.1 , 1.2 , -2.3 /) R = (/ 0.1 , 1.2 , -2.3 /)
@ -3068,29 +3174,43 @@ integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
do i=1,n do i=1,n
test_qmckl_ao_gaussian_vgl = -11 test_qmckl_ao_gaussian_vgl = -11
if (dabs(1.d0 - VGL(i,1) / (& z = dabs(1.d0 - VGL(i,1) / (dexp(-A(i) * r2)) )
dexp(-A(i) * r2) & if ( z > epsilon ) then
)) > epsilon ) return print *, z, epsilon
return
end if
test_qmckl_ao_gaussian_vgl = -12 test_qmckl_ao_gaussian_vgl = -12
if (dabs(1.d0 - VGL(i,2) / (& z = dabs(1.d0 - VGL(i,2) / (&
-2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) & -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) ))
)) > epsilon ) return if ( z > epsilon ) then
print *, z, epsilon
return
end if
test_qmckl_ao_gaussian_vgl = -13 test_qmckl_ao_gaussian_vgl = -13
if (dabs(1.d0 - VGL(i,3) / (& z = dabs(1.d0 - VGL(i,3) / (&
-2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) & -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) ))
)) > epsilon ) return if ( z > epsilon ) then
print *, z, epsilon
return
end if
test_qmckl_ao_gaussian_vgl = -14 test_qmckl_ao_gaussian_vgl = -14
if (dabs(1.d0 - VGL(i,4) / (& z = dabs(1.d0 - VGL(i,4) / (&
-2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) & -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) ))
)) > epsilon ) return if ( z > epsilon ) then
print *, z, epsilon
return
end if
test_qmckl_ao_gaussian_vgl = -15 test_qmckl_ao_gaussian_vgl = -15
if (dabs(1.d0 - VGL(i,5) / (& z = dabs(1.d0 - VGL(i,5) / (&
A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) & A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) ))
)) > epsilon ) return if ( z > epsilon ) then
print *, z, epsilon
return
end if
end do end do
test_qmckl_ao_gaussian_vgl = 0 test_qmckl_ao_gaussian_vgl = 0
@ -5062,8 +5182,8 @@ for (int32_t ldl=3 ; ldl<=5 ; ++ldl) {
for (int32_t k=0 ; k<5 ; ++k) { for (int32_t k=0 ; k<5 ; ++k) {
for (int32_t l=0 ; l<n ; ++l) { for (int32_t l=0 ; l<n ; ++l) {
printf("VGL[%d][%d] = %e %e\n", k, l, VGL0[k][l], VGL1[k][l]); printf("VGL[%d][%d] = %e %e %e\n", k, l, VGL0[k][l], VGL1[k][l], VGL0[k][l]-VGL1[k][l]);
assert( VGL0[k][l] == VGL1[k][l] ); assert( fabs(1.-(VGL0[k][l]+1.e-100)/(VGL1[k][l]+1.e-100)) < 1.e-15 );
} }
} }
} }
@ -5256,6 +5376,7 @@ end function qmckl_compute_ao_vgl_doc_f
| ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs | | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs |
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
#ifdef HAVE_HPC
qmckl_exit_code qmckl_exit_code
qmckl_compute_ao_vgl_hpc_gaussian ( qmckl_compute_ao_vgl_hpc_gaussian (
const qmckl_context context, const qmckl_context context,
@ -5423,40 +5544,61 @@ qmckl_compute_ao_vgl_hpc_gaussian (
const double s4 = s6_*z; const double s4 = s6_*z;
const double s5 = s5_; const double s5 = s5_;
const int64_t k = ao_index[ishell];
double* __restrict__ const ao_vgl_1 = ao_vgl + ipoint*5*ao_num + k;
const int32_t l = shell_ang_mom[ishell]; const int32_t l = shell_ang_mom[ishell];
const int32_t n = lstart[l+1]-lstart[l]; const int32_t n = lstart[l+1]-lstart[l];
const int64_t k = ao_index[ishell];
double* const ao_vgl_1 = &(ao_vgl[ipoint*5*ao_num+k]); double* __restrict__ const ao_vgl_2 = ao_vgl_1 + ao_num;
double* const ao_vgl_2 = &(ao_vgl_1[ ao_num]); double* __restrict__ const ao_vgl_3 = ao_vgl_1 + (ao_num<<1);
double* const ao_vgl_3 = &(ao_vgl_1[2*ao_num]); double* __restrict__ const ao_vgl_4 = ao_vgl_1 + (ao_num<<1) + ao_num;
double* const ao_vgl_4 = &(ao_vgl_1[3*ao_num]); double* __restrict__ const ao_vgl_5 = ao_vgl_1 + (ao_num<<2);
double* const ao_vgl_5 = &(ao_vgl_1[4*ao_num]);
double* poly_vgl_1; double* __restrict__ poly_vgl_1;
double* poly_vgl_2; double* __restrict__ poly_vgl_2;
double* poly_vgl_3; double* __restrict__ poly_vgl_3;
double* poly_vgl_4; double* __restrict__ poly_vgl_4;
double* poly_vgl_5; double* __restrict__ poly_vgl_5;
if (nidx > 0) { if (nidx > 0) {
const double* f = &(ao_factor[k]); const double* __restrict__ f = ao_factor + k;
const int64_t idx = lstart[l]; const int64_t idx = lstart[l];
switch (nucleus_max_ang_mom[inucl]) { switch (nucleus_max_ang_mom[inucl]) {
case 0: case 0:
ao_vgl_1[0] = s1 * f[0];
ao_vgl_2[0] = s2 * f[0];
ao_vgl_3[0] = s3 * f[0];
ao_vgl_4[0] = s4 * f[0];
ao_vgl_5[0] = s5;
break; break;
case 1: case 1:
poly_vgl_1 = &(poly_vgl_l1[0][idx]); poly_vgl_1 = &(poly_vgl_l1[0][idx]);
poly_vgl_2 = &(poly_vgl_l1[1][idx]); poly_vgl_2 = &(poly_vgl_l1[1][idx]);
poly_vgl_3 = &(poly_vgl_l1[2][idx]); poly_vgl_3 = &(poly_vgl_l1[2][idx]);
poly_vgl_4 = &(poly_vgl_l1[3][idx]); poly_vgl_4 = &(poly_vgl_l1[3][idx]);
for (int64_t il=0 ; il<n ; ++il) { break;
case 2:
poly_vgl_1 = &(poly_vgl_l2[0][idx]);
poly_vgl_2 = &(poly_vgl_l2[1][idx]);
poly_vgl_3 = &(poly_vgl_l2[2][idx]);
poly_vgl_4 = &(poly_vgl_l2[3][idx]);
break;
default:
poly_vgl_1 = &(poly_vgl[0][idx]);
poly_vgl_2 = &(poly_vgl[1][idx]);
poly_vgl_3 = &(poly_vgl[2][idx]);
poly_vgl_4 = &(poly_vgl[3][idx]);
poly_vgl_5 = &(poly_vgl[4][idx]);
}
switch (n) {
case(1):
ao_vgl_1[0] = s1 * f[0];
ao_vgl_2[0] = s2 * f[0];
ao_vgl_3[0] = s3 * f[0];
ao_vgl_4[0] = s4 * f[0];
ao_vgl_5[0] = s5;
break;
case (3):
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int il=0 ; il<3 ; ++il) {
ao_vgl_1[il] = poly_vgl_1[il] * s1 * f[il]; ao_vgl_1[il] = poly_vgl_1[il] * s1 * f[il];
ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il]; ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il]; ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
@ -5467,30 +5609,26 @@ qmckl_compute_ao_vgl_hpc_gaussian (
poly_vgl_4[il] * s4 )) * f[il]; poly_vgl_4[il] * s4 )) * f[il];
} }
break; break;
case 2: case(5):
poly_vgl_1 = &(poly_vgl_l2[0][idx]); #ifdef HAVE_OPENMP
poly_vgl_2 = &(poly_vgl_l2[1][idx]); #pragma omp simd
poly_vgl_3 = &(poly_vgl_l2[2][idx]); #endif
poly_vgl_4 = &(poly_vgl_l2[3][idx]); for (int il=0 ; il<5 ; ++il) {
poly_vgl_5 = &(poly_vgl_l2[4][idx]);
for (int64_t il=0 ; il<n ; ++il) {
ao_vgl_1[il] = poly_vgl_1[il] * s1 * f[il]; ao_vgl_1[il] = poly_vgl_1[il] * s1 * f[il];
ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il]; ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il]; ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il]; ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
ao_vgl_5[il] = (poly_vgl_5[il] * s1 + poly_vgl_1[il] * s5 + ao_vgl_5[il] = (poly_vgl_1[il] * s5 +
2.0*(poly_vgl_2[il] * s2 + 2.0*(poly_vgl_2[il] * s2 +
poly_vgl_3[il] * s3 + poly_vgl_3[il] * s3 +
poly_vgl_4[il] * s4 )) * f[il]; poly_vgl_4[il] * s4 )) * f[il];
} }
break; break;
default: default:
poly_vgl_1 = &(poly_vgl[0][idx]); #ifdef HAVE_OPENMP
poly_vgl_2 = &(poly_vgl[1][idx]); #pragma omp simd simdlen(8)
poly_vgl_3 = &(poly_vgl[2][idx]); #endif
poly_vgl_4 = &(poly_vgl[3][idx]); for (int il=0 ; il<n ; ++il) {
poly_vgl_5 = &(poly_vgl[4][idx]);
for (int64_t il=0 ; il<n ; ++il) {
ao_vgl_1[il] = poly_vgl_1[il] * s1 * f[il]; ao_vgl_1[il] = poly_vgl_1[il] * s1 * f[il];
ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il]; ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il]; ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
@ -5502,7 +5640,6 @@ qmckl_compute_ao_vgl_hpc_gaussian (
} }
break; break;
} }
} else { } else {
for (int64_t il=0 ; il<n ; ++il) { for (int64_t il=0 ; il<n ; ++il) {
ao_vgl_1[il] = 0.0; ao_vgl_1[il] = 0.0;
@ -5519,6 +5656,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
} }
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
} }
#endif
#+end_src #+end_src
** Interfaces ** Interfaces
@ -5545,6 +5683,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
double* const ao_vgl ); double* const ao_vgl );
#+end_src #+end_src
#+begin_src c :tangle (eval h_private_func) :comments org #+begin_src c :tangle (eval h_private_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian ( qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian (
const qmckl_context context, const qmckl_context context,
const int64_t ao_num, const int64_t ao_num,
@ -5565,6 +5704,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
const double* expo, const double* expo,
const double* coef_normalized, const double* coef_normalized,
double* const ao_vgl ); double* const ao_vgl );
#endif
#+end_src #+end_src
#+CALL: generate_c_interface(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_doc")) #+CALL: generate_c_interface(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_doc"))

253
org/qmckl_blas.org
View File

@ -46,7 +46,8 @@
#+begin_src c :comments link :tangle (eval c_test) :noweb yes #+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h" #include "qmckl.h"
#include "assert.h" #include <stdio.h>
#include <assert.h>
#ifdef HAVE_CONFIG_H #ifdef HAVE_CONFIG_H
#include "config.h" #include "config.h"
#endif #endif
@ -60,6 +61,18 @@ int main() {
#+end_src #+end_src
* -
:PROPERTIES:
:UNNUMBERED: t
:END:
Basic linear algebra data types and operations are described in this file.
The data types are private, so that HPC implementations can use
whatever data structures they prefer.
These data types are expected to be used internally in QMCkl. They
are not intended to be passed to external codes.
* Data types * Data types
** Vector ** Vector
@ -69,7 +82,7 @@ int main() {
| ~size~ | ~int64_t~ | Dimension of the vector | | ~size~ | ~int64_t~ | Dimension of the vector |
| ~data~ | ~double*~ | Elements | | ~data~ | ~double*~ | Elements |
#+begin_src c :comments org :tangle (eval h_private_type) #+begin_src c :comments org :tangle (eval h_private_type) :exports none
typedef struct qmckl_vector { typedef struct qmckl_vector {
int64_t size; int64_t size;
double* data; double* data;
@ -85,7 +98,7 @@ qmckl_vector_alloc( qmckl_context context,
Allocates a new vector. If the allocation failed the size is zero. Allocates a new vector. If the allocation failed the size is zero.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_vector qmckl_vector
qmckl_vector_alloc( qmckl_context context, qmckl_vector_alloc( qmckl_context context,
const int64_t size) const int64_t size)
@ -114,7 +127,7 @@ qmckl_vector_free( qmckl_context context,
qmckl_vector* vector); qmckl_vector* vector);
#+end_src #+end_src
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_vector_free( qmckl_context context, qmckl_vector_free( qmckl_context context,
qmckl_vector* vector) qmckl_vector* vector)
@ -145,7 +158,7 @@ qmckl_vector_free( qmckl_context context,
The dimensions use Fortran ordering: two elements differing by one The dimensions use Fortran ordering: two elements differing by one
in the first dimension are consecutive in memory. in the first dimension are consecutive in memory.
#+begin_src c :comments org :tangle (eval h_private_type) #+begin_src c :comments org :tangle (eval h_private_type) :exports none
typedef struct qmckl_matrix { typedef struct qmckl_matrix {
int64_t size[2]; int64_t size[2];
double* data; double* data;
@ -162,7 +175,7 @@ qmckl_matrix_alloc( qmckl_context context,
Allocates a new matrix. If the allocation failed the sizes are zero. Allocates a new matrix. If the allocation failed the sizes are zero.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_matrix qmckl_matrix
qmckl_matrix_alloc( qmckl_context context, qmckl_matrix_alloc( qmckl_context context,
const int64_t size1, const int64_t size1,
@ -195,7 +208,7 @@ qmckl_matrix_free( qmckl_context context,
qmckl_matrix* matrix); qmckl_matrix* matrix);
#+end_src #+end_src
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_matrix_free( qmckl_context context, qmckl_matrix_free( qmckl_context context,
qmckl_matrix* matrix) qmckl_matrix* matrix)
@ -228,7 +241,7 @@ qmckl_matrix_free( qmckl_context context,
The dimensions use Fortran ordering: two elements differing by one The dimensions use Fortran ordering: two elements differing by one
in the first dimension are consecutive in memory. in the first dimension are consecutive in memory.
#+begin_src c :comments org :tangle (eval h_private_type) #+begin_src c :comments org :tangle (eval h_private_type) :exports none
#define QMCKL_TENSOR_ORDER_MAX 16 #define QMCKL_TENSOR_ORDER_MAX 16
typedef struct qmckl_tensor { typedef struct qmckl_tensor {
@ -249,7 +262,7 @@ qmckl_tensor_alloc( qmckl_context context,
Allocates memory for a tensor. If the allocation failed, the size Allocates memory for a tensor. If the allocation failed, the size
is zero. is zero.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_tensor qmckl_tensor
qmckl_tensor_alloc( qmckl_context context, qmckl_tensor_alloc( qmckl_context context,
const int64_t order, const int64_t order,
@ -289,7 +302,7 @@ qmckl_tensor_free( qmckl_context context,
qmckl_tensor* tensor); qmckl_tensor* tensor);
#+end_src #+end_src
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_tensor_free( qmckl_context context, qmckl_tensor_free( qmckl_context context,
qmckl_tensor* tensor) qmckl_tensor* tensor)
@ -325,7 +338,7 @@ qmckl_matrix_of_vector(const qmckl_vector vector,
Reshapes a vector into a matrix. Reshapes a vector into a matrix.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_matrix qmckl_matrix
qmckl_matrix_of_vector(const qmckl_vector vector, qmckl_matrix_of_vector(const qmckl_vector vector,
const int64_t size1, const int64_t size1,
@ -355,7 +368,7 @@ qmckl_tensor_of_vector(const qmckl_vector vector,
Reshapes a vector into a tensor. Reshapes a vector into a tensor.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_tensor qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector, qmckl_tensor_of_vector(const qmckl_vector vector,
const int64_t order, const int64_t order,
@ -385,7 +398,7 @@ qmckl_vector_of_matrix(const qmckl_matrix matrix);
Reshapes a matrix into a vector. Reshapes a matrix into a vector.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_vector qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix) qmckl_vector_of_matrix(const qmckl_matrix matrix)
{ {
@ -409,7 +422,7 @@ qmckl_tensor_of_matrix(const qmckl_matrix matrix,
Reshapes a matrix into a tensor. Reshapes a matrix into a tensor.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_tensor qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix, qmckl_tensor_of_matrix(const qmckl_matrix matrix,
const int64_t order, const int64_t order,
@ -439,7 +452,7 @@ qmckl_vector_of_tensor(const qmckl_tensor tensor);
Reshapes a tensor into a vector. Reshapes a tensor into a vector.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_vector qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor) qmckl_vector_of_tensor(const qmckl_tensor tensor)
{ {
@ -468,7 +481,7 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
Reshapes a tensor into a vector. Reshapes a tensor into a vector.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_matrix qmckl_matrix
qmckl_matrix_of_tensor(const qmckl_tensor tensor, qmckl_matrix_of_tensor(const qmckl_tensor tensor,
const int64_t size1, const int64_t size1,
@ -493,13 +506,86 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
** Access macros ** Access macros
#+begin_src c :comments org :tangle (eval h_private_func) Macros are provided to ease the access to vectors, matrices and
tensors. Matrices use column-major ordering, as in Fortran.
#+begin_src c :tangle no
double qmckl_vec (qmckl_vector v, int64_t i); // v[i]
double qmckl_mat (qmckl_matrix m, int64_t i, int64_t j) // m[j][i]
double qmckl_ten3(qmckl_tensor t, int64_t i, int64_t j, int64_t k); // t[k][j][i]
double qmckl_ten4(qmckl_tensor t, int64_t i, int64_t j, int64_t k, int64_t l);
double qmckl_ten5(qmckl_tensor t, int64_t i, int64_t j, int64_t k, int64_t l, int64_t m);
...
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
#define qmckl_vec(v, i) v.data[i] #define qmckl_vec(v, i) v.data[i]
#define qmckl_mat(m, i, j) m.data[(i) + (j)*m.size[0]] #define qmckl_mat(m, i, j) m.data[(i) + (j)*m.size[0]]
#define qmckl_ten3(t, i, j, k) t.data[(i) + m.size[0]*((j) + size[1]*(k))] #define qmckl_ten3(t, i, j, k) t.data[(i) + t.size[0]*((j) + t.size[1]*(k))]
#define qmckl_ten4(t, i, j, k, l) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*(l)))] #define qmckl_ten4(t, i, j, k, l) t.data[(i) + t.size[0]*((j) + t.size[1]*((k) + t.size[2]*(l)))]
#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))))] #define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + t.size[0]*((j) + t.size[1]*((k) + t.size[2]*((l) + t.size[3]*(m))))]
#+end_src
For example:
** Set all elements
*** Vector
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_set(qmckl_vector vector, double value);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_vector
qmckl_vector_set(qmckl_vector vector, double value)
{
for (int64_t i=0 ; i<vector.size ; ++i) {
qmckl_vec(vector, i) = value;
}
return vector;
}
#+end_src
*** Matrix
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_set(qmckl_matrix matrix, double value);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_matrix
qmckl_matrix_set(qmckl_matrix matrix, double value)
{
qmckl_vector vector = qmckl_vector_of_matrix(matrix);
for (int64_t i=0 ; i<vector.size ; ++i) {
qmckl_vec(vector, i) = value;
}
return qmckl_matrix_of_vector(vector, matrix.size[0], matrix.size[1]);
}
#+end_src
*** Tensor
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_set(qmckl_tensor tensor, double value);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_tensor
qmckl_tensor_set(qmckl_tensor tensor, double value)
{
qmckl_vector vector = qmckl_vector_of_tensor(tensor);
for (int64_t i=0 ; i<vector.size ; ++i) {
qmckl_vec(vector, i) = value;
}
return qmckl_tensor_of_vector(vector, tensor.order, tensor.size);
}
#+end_src #+end_src
** Copy to/from to ~double*~ ** Copy to/from to ~double*~
@ -514,7 +600,7 @@ qmckl_double_of_vector(const qmckl_context context,
Converts a vector to a ~double*~. Converts a vector to a ~double*~.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_double_of_vector(const qmckl_context context, qmckl_double_of_vector(const qmckl_context context,
const qmckl_vector vector, const qmckl_vector vector,
@ -546,7 +632,7 @@ qmckl_double_of_matrix(const qmckl_context context,
Converts a matrix to a ~double*~. Converts a matrix to a ~double*~.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_double_of_matrix(const qmckl_context context, qmckl_double_of_matrix(const qmckl_context context,
const qmckl_matrix matrix, const qmckl_matrix matrix,
@ -569,7 +655,7 @@ qmckl_double_of_tensor(const qmckl_context context,
Converts a tensor to a ~double*~. Converts a tensor to a ~double*~.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_double_of_tensor(const qmckl_context context, qmckl_double_of_tensor(const qmckl_context context,
const qmckl_tensor tensor, const qmckl_tensor tensor,
@ -592,7 +678,7 @@ qmckl_vector_of_double(const qmckl_context context,
Converts a ~double*~ to a vector. Converts a ~double*~ to a vector.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_vector_of_double(const qmckl_context context, qmckl_vector_of_double(const qmckl_context context,
const double* target, const double* target,
@ -636,9 +722,9 @@ qmckl_matrix_of_double(const qmckl_context context,
qmckl_matrix* matrix); qmckl_matrix* matrix);
#+end_src #+end_src
Converts a matrix to a ~double*~. Converts a ~double*~ to a matrix.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_matrix_of_double(const qmckl_context context, qmckl_matrix_of_double(const qmckl_context context,
const double* target, const double* target,
@ -662,9 +748,9 @@ qmckl_tensor_of_double(const qmckl_context context,
qmckl_tensor* tensor); qmckl_tensor* tensor);
#+end_src #+end_src
Converts a matrix to a ~double*~. Converts a ~double*~ to a tensor.
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_tensor_of_double(const qmckl_context context, qmckl_tensor_of_double(const qmckl_context context,
const double* target, const double* target,
@ -679,12 +765,9 @@ qmckl_tensor_of_double(const qmckl_context context,
} }
#+end_src #+end_src
** Tests :noexport:
#+begin_src c :comments link :tangle (eval c_test) :exports none
** Tests
#+begin_src c :comments link :tangle (eval c_test)
{ {
int64_t m = 3; int64_t m = 3;
int64_t n = 4; int64_t n = 4;
@ -720,7 +803,7 @@ qmckl_tensor_of_double(const qmckl_context context,
** ~qmckl_dgemm~ ** ~qmckl_dgemm~
Matrix multiplication: Matrix multiplication with a BLAS interface:
\[ \[
C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj} C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
@ -780,7 +863,7 @@ qmckl_tensor_of_double(const qmckl_context context,
const int64_t ldc ); const int64_t ldc );
#+end_src #+end_src
#+begin_src f90 :tangle (eval f) #+begin_src f90 :tangle (eval f) :exports none
integer function qmckl_dgemm_f(context, TransA, TransB, & integer function qmckl_dgemm_f(context, TransA, TransB, &
m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) & m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
result(info) result(info)
@ -975,14 +1058,14 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
end function test_qmckl_dgemm end function test_qmckl_dgemm
#+end_src #+end_src
#+begin_src c :comments link :tangle (eval c_test) #+begin_src c :comments link :tangle (eval c_test) :exports none
qmckl_exit_code test_qmckl_dgemm(qmckl_context context); qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
#+end_src #+end_src
** ~qmckl_matmul~ ** ~qmckl_matmul~
Matrix multiplication: Matrix multiplication using the =qmckl_matrix= data type:
\[ \[
C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj} C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
@ -1017,7 +1100,7 @@ qmckl_matmul (const qmckl_context context,
qmckl_matrix* const C ); qmckl_matrix* const C );
#+end_src #+end_src
#+begin_src c :tangle (eval c) :comments org #+begin_src c :tangle (eval c) :comments org :exports none
qmckl_exit_code qmckl_exit_code
qmckl_matmul (const qmckl_context context, qmckl_matmul (const qmckl_context context,
const char TransA, const char TransA,
@ -1155,9 +1238,87 @@ qmckl_matmul (const qmckl_context context,
} }
#+end_src #+end_src
*** Test :noexport:
*** TODO Test :noexport: #+begin_src python :exports none :results output
import numpy as np
A = np.array([[ 1., 2., 3., 4. ],
[ 5., 6., 7., 8. ],
[ 9., 10., 11., 12. ]])
B = np.array([[ 1., -2., 3., 4., 5. ],
[ 5., -6., 7., 8., 9. ],
[ 9., 10., 11., 12., 13. ],
[10., 11., 12., 15., 14. ]])
C = 0.5 * A @ B
print(A.T)
print(B.T)
print(C.T)
#+end_src
#+RESULTS:
#+begin_example
[[ 1. 5. 9.]
[ 2. 6. 10.]
[ 3. 7. 11.]
[ 4. 8. 12.]]
[[ 1. 5. 9. 10.]
[-2. -6. 10. 11.]
[ 3. 7. 11. 12.]
[ 4. 8. 12. 15.]
[ 5. 9. 13. 14.]]
[[ 39. 89. 139.]
[ 30. 56. 82.]
[ 49. 115. 181.]
[ 58. 136. 214.]
[ 59. 141. 223.]]
#+end_example
#+begin_src c :comments link :tangle (eval c_test) :exports none
{
double a[12] = { 1., 5., 9.,
2., 6., 10.,
3., 7., 11.,
4., 8., 12. };
double b[20] = { 1., 5., 9., 10.,
-2., -6., 10., 11.,
3., 7., 11., 12.,
4., 8., 12., 15.,
5., 9., 13., 14. };
double c[15] = { 39., 89., 139.,
30., 56., 82.,
49., 115., 181.,
58., 136., 214.,
59., 141., 223. };
double cnew[15];
qmckl_exit_code rc;
qmckl_matrix A = qmckl_matrix_alloc(context, 3, 4);
rc = qmckl_matrix_of_double(context, a, 12, &A);
assert(rc == QMCKL_SUCCESS);
qmckl_matrix B = qmckl_matrix_alloc(context, 4, 5);
rc = qmckl_matrix_of_double(context, b, 20, &B);
assert(rc == QMCKL_SUCCESS);
qmckl_matrix C = qmckl_matrix_alloc(context, 3, 5);
rc = qmckl_matmul(context, 'N', 'N', 0.5, A, B, 0., &C);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_double_of_matrix(context, C, cnew, 15);
assert(rc == QMCKL_SUCCESS);
for (int i=0 ; i<15 ; ++i) {
printf("%f %f\n", cnew[i], c[i]);
assert (c[i] == cnew[i]);
}
}
#+end_src
** ~qmckl_adjugate~ ** ~qmckl_adjugate~
Given a matrix $\mathbf{A}$, the adjugate matrix Given a matrix $\mathbf{A}$, the adjugate matrix
@ -1208,7 +1369,7 @@ qmckl_matmul (const qmckl_context context,
for performance motivations. For larger sizes, we rely on the for performance motivations. For larger sizes, we rely on the
LAPACK library. LAPACK library.
#+begin_src f90 :tangle (eval f) #+begin_src f90 :tangle (eval f) :exports none
integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) & integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) &
result(info) result(info)
use qmckl use qmckl
@ -2048,10 +2209,10 @@ assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
Transposes a matrix: $A^\dagger_{ji} = A_{ij}$. Transposes a matrix: $A^\dagger_{ji} = A_{ij}$.
| Variable | Type | In/Out | Description | | Variable | Type | In/Out | Description |
|----------+---------------+--------+-------------------| |-----------+-----------------+--------+-------------------|
| context | qmckl_context | in | Global state | | ~context~ | ~qmckl_context~ | in | Global state |
| A | qmckl_matrix | in | Input matrix | | ~A~ | ~qmckl_matrix~ | in | Input matrix |
| At | qmckl_matrix | out | Transposed matrix | | ~At~ | ~qmckl_matrix~ | out | Transposed matrix |
#+begin_src c :tangle (eval h_private_func) :comments org #+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_exit_code
@ -2061,7 +2222,7 @@ qmckl_transpose (qmckl_context context,
#+end_src #+end_src
#+begin_src c :tangle (eval c) :comments org #+begin_src c :tangle (eval c) :comments org :exports none
qmckl_exit_code qmckl_exit_code
qmckl_transpose (qmckl_context context, qmckl_transpose (qmckl_context context,
const qmckl_matrix A, const qmckl_matrix A,
@ -2100,7 +2261,7 @@ qmckl_transpose (qmckl_context context,
} }
#+end_src #+end_src
*** Test *** Test :noexport:
#+begin_src c :comments link :tangle (eval c_test) #+begin_src c :comments link :tangle (eval c_test)
{ {

75
org/qmckl_context.org
View File

@ -62,6 +62,7 @@ int main() {
#include <stdio.h> #include <stdio.h>
#include <string.h> #include <string.h>
#include <errno.h> #include <errno.h>
#include <pthread.h>
#include "qmckl.h" #include "qmckl.h"
#include "qmckl_context_private_type.h" #include "qmckl_context_private_type.h"
@ -69,6 +70,8 @@ int main() {
#+end_src #+end_src
* Context handling * Context handling
The context variable is a handle for the state of the library, The context variable is a handle for the state of the library,
@ -146,7 +149,7 @@ typedef struct qmckl_context_struct {
the pointer associated with a context is non-null, we can still the pointer associated with a context is non-null, we can still
verify that it points to the expected data structure. verify that it points to the expected data structure.
#+begin_src c :comments org :tangle (eval h_private_type) :noweb yes #+begin_src c :comments org :tangle (eval h_private_type) :noweb yes :exports none
#define VALID_TAG 0xBEEFFACE #define VALID_TAG 0xBEEFFACE
#define INVALID_TAG 0xDEADBEEF #define INVALID_TAG 0xDEADBEEF
#+end_src #+end_src
@ -156,10 +159,11 @@ typedef struct qmckl_context_struct {
if the context is valid, ~QMCKL_NULL_CONTEXT~ otherwise. if the context is valid, ~QMCKL_NULL_CONTEXT~ otherwise.
#+begin_src c :comments org :tangle (eval h_func) :noexport #+begin_src c :comments org :tangle (eval h_func) :noexport
qmckl_context qmckl_context_check(const qmckl_context context) ; qmckl_context
qmckl_context_check (const qmckl_context context) ;
#+end_src #+end_src
#+begin_src c :tangle (eval c) #+begin_src c :tangle (eval c) :exports none
qmckl_context qmckl_context_check(const qmckl_context context) { qmckl_context qmckl_context_check(const qmckl_context context) {
if (context == QMCKL_NULL_CONTEXT) if (context == QMCKL_NULL_CONTEXT)
@ -176,19 +180,27 @@ qmckl_context qmckl_context_check(const qmckl_context context) {
} }
#+end_src #+end_src
The context keeps a ``date'' that allows to check which data needs The context keeps a /date/ that allows to check which data needs
to be recomputed. The date is incremented when the context is touched. to be recomputed. The date is incremented when the context is touched.
When a new element is added to the context, the functions When a new element is added to the context, the functions
[[Creation][qmckl_context_create]], [[Destroy][qmckl_context_destroy]] and [[Copy][qmckl_context_copy]] [[Creation][=qmckl_context_create=]] [[Destroy][=qmckl_context_destroy=]] and [[Copy][=qmckl_context_copy=]]
should be updated in order to make deep copies. should be updated in order to make deep copies.
#+begin_src c :comments org :tangle (eval h_func) :noexport When the electron coordinates have changed, the context is touched
qmckl_exit_code qmckl_context_touch(const qmckl_context context) ; using the following function.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_context_touch (const qmckl_context context);
#+end_src #+end_src
#+begin_src c :tangle (eval c) This has the effect to increment the date of the context.
qmckl_exit_code qmckl_context_touch(const qmckl_context context) {
#+begin_src c :tangle (eval c) :exports none
qmckl_exit_code
qmckl_context_touch(const qmckl_context context)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context, return qmckl_failwith( context,
@ -213,12 +225,12 @@ qmckl_exit_code qmckl_context_touch(const qmckl_context context) {
- A new context always has all its members initialized with a NULL value - A new context always has all its members initialized with a NULL value
# Header # Header
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func)
qmckl_context qmckl_context_create(); qmckl_context qmckl_context_create();
#+end_src #+end_src
# Source # Source
#+begin_src c :tangle (eval c) #+begin_src c :tangle (eval c) :exports none
qmckl_context qmckl_context_create() { qmckl_context qmckl_context_create() {
qmckl_context_struct* const ctx = qmckl_context_struct* const ctx =
@ -241,7 +253,9 @@ qmckl_context qmckl_context_create() {
rc = pthread_mutexattr_init(&attr); rc = pthread_mutexattr_init(&attr);
assert (rc == 0); assert (rc == 0);
#ifdef PTHREAD_MUTEX_RECURSIVE
(void) pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE); (void) pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE);
#endif
rc = pthread_mutex_init ( &(ctx->mutex), &attr); rc = pthread_mutex_init ( &(ctx->mutex), &attr);
assert (rc == 0); assert (rc == 0);
@ -326,13 +340,13 @@ assert( qmckl_context_check(context) == context );
~lock_count~ attribute. ~lock_count~ attribute.
# Header # Header
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func)
void qmckl_lock (qmckl_context context); void qmckl_lock (qmckl_context context);
void qmckl_unlock(qmckl_context context); void qmckl_unlock(qmckl_context context);
#+end_src #+end_src
# Source # Source
#+begin_src c :tangle (eval c) #+begin_src c :tangle (eval c) :exports none
void qmckl_lock(qmckl_context context) { void qmckl_lock(qmckl_context context) {
if (context == QMCKL_NULL_CONTEXT) if (context == QMCKL_NULL_CONTEXT)
return ; return ;
@ -345,9 +359,6 @@ void qmckl_lock(qmckl_context context) {
} }
assert (rc == 0); assert (rc == 0);
ctx->lock_count += 1; ctx->lock_count += 1;
/*
printf(" lock : %d\n", ctx->lock_count);
,*/
} }
void qmckl_unlock(const qmckl_context context) { void qmckl_unlock(const qmckl_context context) {
@ -359,9 +370,6 @@ void qmckl_unlock(const qmckl_context context) {
} }
assert (rc == 0); assert (rc == 0);
ctx->lock_count -= 1; ctx->lock_count -= 1;
/*
printf("unlock : %d\n", ctx->lock_count);
,*/
} }
#+end_src #+end_src
@ -372,12 +380,14 @@ void qmckl_unlock(const qmckl_context context) {
# Header # Header
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func) :exports none
/*
qmckl_context qmckl_context_copy(const qmckl_context context); qmckl_context qmckl_context_copy(const qmckl_context context);
*/
#+end_src #+end_src
# Source # Source
#+begin_src c :tangle (eval c) #+begin_src c :tangle (eval c) :exports none
qmckl_context qmckl_context_copy(const qmckl_context context) { qmckl_context qmckl_context_copy(const qmckl_context context) {
const qmckl_context checked_context = qmckl_context_check(context); const qmckl_context checked_context = qmckl_context_check(context);
@ -418,13 +428,13 @@ qmckl_context qmckl_context_copy(const qmckl_context context) {
# Fortran interface # Fortran interface
#+begin_src f90 :tangle (eval fh_func) :exports none #+begin_src f90 :tangle (eval fh_func) :exports none
interface ! interface
integer (qmckl_context) function qmckl_context_copy(context) bind(C) ! integer (qmckl_context) function qmckl_context_copy(context) bind(C)
use, intrinsic :: iso_c_binding ! use, intrinsic :: iso_c_binding
import ! import
integer (qmckl_context), intent(in), value :: context ! integer (qmckl_context), intent(in), value :: context
end function qmckl_context_copy ! end function qmckl_context_copy
end interface ! end interface
#+end_src #+end_src
# Test # Test
@ -443,13 +453,16 @@ munit_assert_int64(qmckl_context_check(new_context), ==, new_context);
It frees the context, and returns the previous context. It frees the context, and returns the previous context.
# Header # Header
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_context_destroy(const qmckl_context context); qmckl_exit_code
qmckl_context_destroy (const qmckl_context context);
#+end_src #+end_src
# Source # Source
#+begin_src c :tangle (eval c) #+begin_src c :tangle (eval c) :exports none
qmckl_exit_code qmckl_context_destroy(const qmckl_context context) { qmckl_exit_code
qmckl_context_destroy (const qmckl_context context)
{
const qmckl_context checked_context = qmckl_context_check(context); const qmckl_context checked_context = qmckl_context_check(context);
if (checked_context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT; if (checked_context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;

49
org/qmckl_error.org
View File

@ -46,7 +46,7 @@ int main() {
#+end_src #+end_src
* * -
:PROPERTIES: :PROPERTIES:
:UNNUMBERED: t :UNNUMBERED: t
:END: :END:
@ -204,17 +204,24 @@ return '\n'.join(result)
string is assumed to be large enough to contain the error message string is assumed to be large enough to contain the error message
(typically 128 characters). (typically 128 characters).
* hidden :noexport:
#+NAME: MAX_STRING_LENGTH
: 128
* Decoding errors * Decoding errors
To decode the error messages, ~qmckl_string_of_error~ converts an To decode the error messages, ~qmckl_string_of_error~ converts an
error code into a string. error code into a string.
#+NAME: MAX_STRING_LENGTH #+begin_src c :comments org :tangle (eval h_func) :noweb yes
: 128 const char*
qmckl_string_of_error (const qmckl_exit_code error);
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none :noweb yes #+begin_src c :comments org :tangle (eval h_private_func) :exports none :noweb yes
const char* qmckl_string_of_error(const qmckl_exit_code error); void
void qmckl_string_of_error_f(const qmckl_exit_code error, qmckl_string_of_error_f(const qmckl_exit_code error,
char result[<<MAX_STRING_LENGTH()>>]); char result[<<MAX_STRING_LENGTH()>>]);
#+end_src #+end_src
@ -326,7 +333,7 @@ return '\n'.join(result)
#+end_example #+end_example
# Source # Source
#+begin_src c :comments org :tangle (eval c) :noweb yes #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
const char* qmckl_string_of_error(const qmckl_exit_code error) { const char* qmckl_string_of_error(const qmckl_exit_code error) {
switch (error) { switch (error) {
<<cases()>> <<cases()>>
@ -340,7 +347,7 @@ void qmckl_string_of_error_f(const qmckl_exit_code error, char result[<<MAX_STRI
#+end_src #+end_src
# Fortran interface # Fortran interface
#+begin_src f90 :tangle (eval fh_func) :noexport :noweb yes #+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes
interface interface
subroutine qmckl_string_of_error (error, string) bind(C, name='qmckl_string_of_error_f') subroutine qmckl_string_of_error (error, string) bind(C, name='qmckl_string_of_error_f')
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
@ -353,7 +360,7 @@ void qmckl_string_of_error_f(const qmckl_exit_code error, char result[<<MAX_STRI
* Data structure in context * Data structure in context
The strings are declared with a maximum fixed size to avoid The strings are declared internally with a maximum fixed size to avoid
dynamic memory allocation. dynamic memory allocation.
#+begin_src c :comments org :tangle (eval h_private_type) #+begin_src c :comments org :tangle (eval h_private_type)
@ -377,7 +384,7 @@ typedef struct qmckl_error_struct {
explaining the error. The exit code can't be ~QMCKL_SUCCESS~. explaining the error. The exit code can't be ~QMCKL_SUCCESS~.
# Header # Header
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_exit_code
qmckl_set_error(qmckl_context context, qmckl_set_error(qmckl_context context,
const qmckl_exit_code exit_code, const qmckl_exit_code exit_code,
@ -386,7 +393,7 @@ qmckl_set_error(qmckl_context context,
#+end_src #+end_src
# Source # Source
#+begin_src c :tangle (eval c) #+begin_src c :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_set_error(qmckl_context context, qmckl_set_error(qmckl_context context,
const qmckl_exit_code exit_code, const qmckl_exit_code exit_code,
@ -428,7 +435,7 @@ qmckl_set_error(qmckl_context context,
function name and message is mandatory. function name and message is mandatory.
# Header # Header
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_exit_code
qmckl_get_error(qmckl_context context, qmckl_get_error(qmckl_context context,
qmckl_exit_code *exit_code, qmckl_exit_code *exit_code,
@ -437,7 +444,7 @@ qmckl_get_error(qmckl_context context,
#+end_src #+end_src
# Source # Source
#+begin_src c :tangle (eval c) #+begin_src c :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_get_error(qmckl_context context, qmckl_get_error(qmckl_context context,
qmckl_exit_code *exit_code, qmckl_exit_code *exit_code,
@ -488,19 +495,21 @@ qmckl_get_error(qmckl_context context,
the desired return code. the desired return code.
Upon failure, a ~QMCKL_NULL_CONTEXT~ is returned. Upon failure, a ~QMCKL_NULL_CONTEXT~ is returned.
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_failwith(qmckl_context context, qmckl_exit_code
qmckl_failwith(qmckl_context context,
const qmckl_exit_code exit_code, const qmckl_exit_code exit_code,
const char* function, const char* function,
const char* message) ; const char* message) ;
#+end_src #+end_src
#+begin_src c :comments org :tangle (eval c) #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_failwith(qmckl_context context, qmckl_exit_code
qmckl_failwith(qmckl_context context,
const qmckl_exit_code exit_code, const qmckl_exit_code exit_code,
const char* function, const char* function,
const char* message) { const char* message)
{
assert (exit_code > 0); assert (exit_code > 0);
assert (exit_code < QMCKL_INVALID_EXIT_CODE); assert (exit_code < QMCKL_INVALID_EXIT_CODE);
assert (function != NULL); assert (function != NULL);
@ -544,8 +553,6 @@ if (x < 0) {
#endif #endif
#+end_src #+end_src
** Test
#+begin_src c :comments link :tangle (eval c_test) #+begin_src c :comments link :tangle (eval c_test)
/* Initialize the variables */ /* Initialize the variables */
char function_name[QMCKL_MAX_FUN_LEN]=""; char function_name[QMCKL_MAX_FUN_LEN]="";

62
org/qmckl_nucleus.org
View File

@ -187,6 +187,19 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
} }
#+end_src #+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_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(out) :: num
end function qmckl_get_nucleus_num
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_exit_code
@ -242,6 +255,20 @@ qmckl_get_nucleus_charge (const qmckl_context context,
#+end_src #+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_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(out) :: charge(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_nucleus_charge
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_exit_code
@ -277,6 +304,18 @@ qmckl_get_nucleus_rescale_factor (const qmckl_context context,
} }
#+end_src #+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_nucleus_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(out) :: kappa
end function qmckl_get_nucleus_rescale_factor
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_exit_code
@ -342,6 +381,21 @@ qmckl_get_nucleus_coord (const qmckl_context context,
} }
#+end_src #+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_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(out) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_nucleus_coord
end interface
#+end_src
When all the data relative to nuclei have been set, the following When all the data relative to nuclei have been set, the following
function returns ~true~. function returns ~true~.
@ -424,7 +478,7 @@ interface
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: num integer (c_int64_t) , intent(in) , value :: num
end function end function qmckl_set_nucleus_num
end interface end interface
#+end_src #+end_src
@ -570,7 +624,7 @@ interface
character(c_char) , intent(in) , value :: transp character(c_char) , intent(in) , value :: transp
real (c_double) , intent(in) :: coord(*) real (c_double) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function end function qmckl_set_nucleus_coord
end interface end interface
#+end_src #+end_src
@ -604,14 +658,14 @@ qmckl_set_nucleus_rescale_factor(qmckl_context context,
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface interface
integer(c_int32_t) function qmckl_set_rescale_factor(context, kappa) & integer(c_int32_t) function qmckl_set_nucleus_rescale_factor(context, kappa) &
bind(C) bind(C)
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
import import
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) , value :: kappa real (c_double) , intent(in) , value :: kappa
end function end function qmckl_set_nucleus_rescale_factor
end interface end interface
#+end_src #+end_src