TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-07-18 08:53:47 +02:00
Code Issues Projects Releases Wiki Activity

Merge branch 'master' into wf_det_grad_cof

Browse Source
This commit is contained in:
vijay 2021-10-13 16:18:21 +02:00 committed by GitHub
commit 19b4f93a0b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
14 changed files with 767 additions and 143 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
.github/workflows/vfc_test_workflow.yml vendored Normal file
View File

@ -0,0 +1,51 @@
# This workflow will be executed when master is updated:
# it will run the configured tests and upload the results on vfc_ci_master.
name: "Verificarlo CI (master)"
on:
# Triggers the workflow when master is updated
push:
branches: [ master ]
workflow_dispatch:
jobs:
run_verificarlo_tests:
runs-on: ubuntu-latest
container: verificarlo/verificarlo
steps:
- uses: actions/checkout@v2
with:
fetch-depth: 0
- name: Install dependencies
run: |
ln -s /usr/bin/python3 /usr/bin/python
apt update
apt -y install emacs pkg-config
- name: Run tests
run: vfc_ci test -g -r
- name: Commit test results
run: |
git_hash=$(git rev-parse --short "$GITHUB_SHA")
git config --local user.email "action@github.com"
git config --local user.name "GitHub Action"
git checkout vfc_ci_master
mkdir -p vfcruns
mv *.vfcrun.h5 vfcruns
git add vfcruns/*
git commit -m "[auto] New test results for commit ${git_hash}"
git push
- name: Upload raw results as artifacts
uses: actions/upload-artifact@v2
with:
name: ${{github.sha}}.vfcraw
path: ./*.vfcraw.h5

3
.gitignore vendored
View File

@ -13,6 +13,8 @@ config.status
config.sub
configure
generated.mk
.vfcwrapper.o
libtool
m4/libtool.m4
m4/ltoptions.m4
m4/ltsugar.m4
@ -20,6 +22,7 @@ m4/ltversion.m4
m4/lt~obsolete.m4
qmckl-*.tar.gz
qmckl.mod
qmckl_probes_f.mod
qmckl.pc
stamp-h1
tools/compile

26
Makefile.am
View File

@ -30,6 +30,9 @@
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
if VFC_CI
AM_LDFLAGS=-lvfc_probes -lvfc_probes_f
endif
ACLOCAL_AMFLAGS = -I m4
@ -38,7 +41,7 @@ VERSION_MINOR = @VERSION_MINOR@
VERSION_PATCH = @VERSION_PATCH@
SUBDIRS =
CLEANFILES = qmckl.mod
CLEANFILES = qmckl.mod qmckl_probes_f.mod
EXTRA_DIST =
pkgconfigdir = $(libdir)/pkgconfig
@ -53,6 +56,8 @@ test_qmckl_fo = tests/qmckl_f.o
src_qmckl_f = src/qmckl_f.f90
src_qmckl_fo = src/qmckl_f.o
header_tests = tests/chbrclf.h tests/n2.h
qmckl_probes_src = src/qmckl_probes.h src/qmckl_probes.c src/qmckl_probes_f.f90
fortrandir = $(datadir)/$(PACKAGE_NAME)/fortran/
dist_fortran_DATA = $(qmckl_f)
@ -60,11 +65,11 @@ dist_fortran_DATA = $(qmckl_f)
AM_CPPFLAGS = -I$(srcdir)/src -I$(srcdir)/include
lib_LTLIBRARIES = src/libqmckl.la
src_libqmckl_la_SOURCES = $(qmckl_h) $(src_qmckl_f) $(C_FILES) $(F_FILES) $(H_PRIVATE_FUNC_FILES) $(H_PRIVATE_TYPE_FILES) $(header_tests)
src_libqmckl_la_SOURCES = $(qmckl_h) $(src_qmckl_f) $(C_FILES) $(F_FILES) $(H_PRIVATE_FUNC_FILES) $(H_PRIVATE_TYPE_FILES) $(header_tests) $(qmckl_probes_src)
export qmckl_f qmckl_h srcdir
CLEANFILES+=$(test_qmckl_f) $(src_qmckl_f) $(test_qmckl_o) $(src_qmckl_o)
CLEANFILES+=$(test_qmckl_f) $(src_qmckl_f) $(test_qmckl_o) $(src_qmckl_o)
htmlize_el=share/doc/qmckl/html/htmlize.el
@ -91,6 +96,15 @@ $(src_qmckl_fo): $(src_qmckl_f)
$(src_qmckl_f): $(srcdir)/$(qmckl_f)
cp $(srcdir)/$(qmckl_f) $(src_qmckl_f)
src/qmckl_probes.c:
cp $(srcdir)/tools/qmckl_probes.c $(srcdir)/src/qmckl_probes.c
src/qmckl_probes.h:
cp $(srcdir)/tools/qmckl_probes.h $(srcdir)/src/qmckl_probes.h
src/qmckl_probes_f.f90:
cp $(srcdir)/tools/qmckl_probes_f.f90 $(srcdir)/src/qmckl_probes_f.f90
share/doc/qmckl/html/index.html: share/doc/qmckl/html/README.html
$(ln_s_verbose)cd share/doc/qmckl/html/ && \
rm -rf index.html && \
@ -109,7 +123,7 @@ dist_src_DATA = $(ORG_FILES) $(TANGLED_FILES) $(EXPORTED_FILES)
BUILT_SOURCES = $(C_FILES) $(F_FILES) $(FH_FUNC_FILES) $(FH_TYPE_FILES) $(H_FUNC_FILES) $(H_TYPE_FILES) $(H_PRIVATE_FUNC_FILES) $(H_PRIVATE_TYPE_FILES) $(qmckl_f) $(qmckl_h) $(header_tests)
CLEANFILES += $(BUILT_SOURCES) $(C_TEST_FILES) $(F_TEST_FILES) $(TANGLED_FILES) $(C_TEST_FILES) $(F_TEST_FILES) $(qmckl_f) $(qmckl_h) $(HTML_FILES) $(TEXT_FILES) share/doc/qmckl/html/index.html $(EXPORTED_FILES) $(header_tests)
CLEANFILES += $(BUILT_SOURCES) $(C_TEST_FILES) $(F_TEST_FILES) $(TANGLED_FILES) $(C_TEST_FILES) $(F_TEST_FILES) $(qmckl_f) $(qmckl_h) $(HTML_FILES) $(TEXT_FILES) share/doc/qmckl/html/index.html $(EXPORTED_FILES) $(header_tests)
EXTRA_DIST += \
tools/build_doc.sh \
@ -156,14 +170,14 @@ $(htmlize_el):
$(srcdir)/tools/install_htmlize.sh $(htmlize_el)
tests/chbrclf.h: $(qmckl_h)
tests/n2.h: $(qmckl_h)
generated.mk: $(ORG_FILES)
$(PYTHON) $(srcdir)/tools/build_makefile.py
cppcheck: cppcheck.out
cppcheck.out: $(qmckl_h)
cd src/ && \
cppcheck --addon=cert -q --error-exitcode=0 \

18
ci_install.sh Executable file
View File

@ -0,0 +1,18 @@
#!/bin/bash
# This scripts is meant to be used by Verificarlo CI to automatically install
# the library dependencies and build it with Verificarlo with vfc_probes support
# enabled.
./autogen.sh
QMCKL_DEVEL=1 ./configure --prefix=$PWD/_install \
--enable-silent-rules --enable-maintainer-mode --enable-vfc_ci --host=x86_64 \
CC="verificarlo-f" FC="verificarlo-f"
make all
# Here we build the test suite, but expect it to fail because it is run without
# specifying VFC_BACKENDS. However, the generated executables will be reused
# individually by the CI.
make check
exit 0

38
configure.ac
View File

@ -4,24 +4,24 @@
# QMCkl - Quantum Monte Carlo kernel library
#
# BSD 3-Clause License
#
#
# Copyright (c) 2020, TREX Center of Excellence
# All rights reserved.
#
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#
# 1. Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimer.
#
#
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.
#
#
# 3. Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from
# this software without specific prior written permission.
#
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@ -167,7 +167,7 @@ AC_TYPE_UINT64_T
# Checks for library functions.
## qmckl
AC_FUNC_MALLOC
# AC_FUNC_MALLOC
AC_CHECK_FUNCS([memset strerror])
# Development mode
@ -196,6 +196,29 @@ if test "x${QMCKL_DEVEL}" != "x"; then
fi
# Enable Verificarlo tests
AC_ARG_ENABLE([vfc_ci],
[ --enable-vfc_ci Build the library with vfc_ci support],
[case "${enableval}" in
yes) vfc_ci=true && FCFLAGS="-D VFC_CI $FCFLAGS" && CFLAGS="-D VFC_CI $CFLAGS";;
no) vfc_ci=false ;;
*) AC_MSG_ERROR([bad value ${enableval} for --enable_vfc_ci]) ;;
esac],[vfc_ci=false])
AM_CONDITIONAL([VFC_CI], [test x$vfc_ci = xtrue])
# Enable Fortran preprocessor
if test "$FC" = "gfortran"; then
AC_MSG_NOTICE(gfortran detected)
# Arguments order is important here
FCFLAGS="-cpp $FCFLAGS"
fi
if test "$FC" = "verificarlo-f"; then
AC_MSG_NOTICE(verificarlo-f detected)
# Arguments order is important here
FCFLAGS="-Mpreprocess $FCFLAGS"
fi
#PKG-CONFIG
#mkl-dynamic-lp64-seq
@ -237,4 +260,3 @@ where the optional <target> is:
check - run tests
install - install ${PACKAGE_NAME}
--------------------------------------------------"

7
include/config.h.in
View File

@ -30,10 +30,6 @@
/* Define to 1 if you have the `pthread' library (-lpthread). */
#undef HAVE_LIBPTHREAD
/* Define to 1 if your system has a GNU libc compatible `malloc' function, and
to 0 otherwise. */
#undef HAVE_MALLOC
/* Define to 1 if you have the <malloc.h> header file. */
#undef HAVE_MALLOC_H
@ -145,9 +141,6 @@
such a type exists and the standard includes do not define it. */
#undef int64_t
/* Define to rpl_malloc if the replacement function should be used. */
#undef malloc
/* Define to `unsigned int' if <sys/types.h> does not define. */
#undef size_t

48
org/qmckl.org
View File

@ -256,6 +256,53 @@ cppcheck --addon=cert --enable=all *.c &> cppcheck.out
functions will use the precision specified in the =context=
variable.
In order to automatize numerical accuracy tests, QMCkl uses
[[https://github.com/verificarlo/verificarlo][Verificarlo]] and
its CI functionality. You can read Verificarlo CI's documentation
at the [[https://github.com/verificarlo/verificarlo/blob/master/doc/06-Postprocessing.md#verificarlo-ci][following link]].
Reading it is advised to understand the remainder of this section.
To enable support for Verificarlo CI tests when building the
library, use the following configure command :
#+begin_src bash
QMCKL_DEVEL=1 ./configure --prefix=$PWD/_install --enable-silent-rules --enable-maintainer-mode CC=verificarlo-f FC=verificarlo-f --host=x86_64 --enable-vfc_ci
#+end_src
Note that this does require an install of Verificarlo *with
Fortran support*. Enabling the support for CI will define the
~VFC_CI~ preprocessor variable which use will be explained now.
As explained in the documentation, Verificarlo CI uses a probes
system to export variables from test programs to the tools itself.
To make tests easier to use, QMCkl has its own interface to the
probes system. To make the system very easy to use, these functions
are always defined, but will behave differently depending on the
~VFC_CI~ preprocessor variable. There are 3 functions at your
disposal. When the CI is enabled, they will place a ~vfc_ci~ probe
as if calling ~vfc_probes~ directly. Otherwise, they will either do
nothing or perform a check on the tested value and return its result
as a boolean that you are free to use or ignore.
Here are these 3 functions :
- ~qmckl_probe~ : place a normal probe witout any check. Won't do anything when ~vfc_ci~ is disabled (false is returned).
- ~qmckl_probe_check~ : place a probe with an absolute check. If ~vfc_ci~ is disabled, this will return the result of an absolute check (|val - ref| < accuracy target ?). If the check fails, true is returned (false otherwise).
- ~qmckl_probe_check_relative~ : place a probe with a relative check. If ~vfc_ci~ is disabled, this will return the result of a relative check (|val - ref| / ref < accuracy target?). If the check fails, true is returned (false otherwise).
If you need more details on these functions or their Fortran
interfaces, have a look at the ~tools/qmckl_probes~ files.
Finally, if you need to add a QMCkl kernel to the CI tests
or modify an existing one, you should pay attention to the
following points :
- you should add the new kernel to the ~vfc_tests_config.json~ file, which controls the backends and repetitions for each executable. More details can be found in the ~vfc_ci~ documentation.
- in order to call the ~qmckl_probes~ functions from Fortran, import the ~qmckl_probes_f~ module.
- if your tests include some asserts that rely on accurate FP values, you should probably wrap them inside a ~#ifndef VFC_CI~ statement, as the asserts would otherwise risk to fail when executed with the different Verificarlo backends.
** Algorithms
Reducing the scaling of an algorithm usually implies also reducing
@ -264,4 +311,3 @@ cppcheck --addon=cert --enable=all *.c &> cppcheck.out
algorithms are better adapted than linear scaling algorithms. As
QMCkl is a general purpose library, multiple algorithms should be
implemented adapted to different problem sizes.

197
org/qmckl_ao.org
View File

@ -117,7 +117,7 @@ int main() {
| ~ao_shell~ | ~[ao_num]~ | For each AO, specify to which shell it belongs |
Computed data:
|--------------------------+----------------------------+-----------------------------------------------------------------------------------------------|
| ~coefficient_normalized~ | ~[prim_num]~ | Normalized primitive coefficients |
| ~nucleus_prim_index~ | ~[nucl_num]~ | Index of the first primitive for each nucleus |
@ -193,7 +193,7 @@ prim_factor = [ 1.0006253235944540e+01, 2.4169531573445120e+00, 7.96109248497664
typedef struct qmckl_ao_basis_struct {
int64_t shell_num;
int64_t prim_num;
int64_t ao_num;
int64_t ao_num;
int64_t * nucleus_index;
int64_t * nucleus_shell_num;
int32_t * shell_ang_mom;
@ -230,10 +230,10 @@ typedef struct qmckl_ao_basis_struct {
Some values are initialized by default, and are not concerned by
this mechanism.
#+begin_src c :comments org :tangle (eval h_private_func)
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code qmckl_init_ao_basis(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_ao_basis(qmckl_context context) {
@ -252,7 +252,7 @@ qmckl_exit_code qmckl_init_ao_basis(qmckl_context context) {
return QMCKL_SUCCESS;
}
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
@ -1180,7 +1180,7 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
}
ctx->ao_basis.nucleus_prim_index[nucl_num] = ctx->ao_basis.prim_num;
}
/* Normalize coefficients */
{
@ -1200,34 +1200,34 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
for (int64_t iprim=ctx->ao_basis.shell_prim_index[ishell] ;
iprim < ctx->ao_basis.shell_prim_index[ishell]+ctx->ao_basis.shell_prim_num[ishell] ;
++iprim) {
ctx->ao_basis.coefficient_normalized[iprim] =
ctx->ao_basis.coefficient_normalized[iprim] =
ctx->ao_basis.coefficient[iprim] * ctx->ao_basis.prim_factor[iprim] *
ctx->ao_basis.shell_factor[ishell];
}
}
}
/* Find max angular momentum on each nucleus */
{
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->nucleus.num * sizeof(int32_t);
ctx->ao_basis.nucleus_max_ang_mom = (int32_t *) qmckl_malloc(context, mem_info);
if (ctx->ao_basis.nucleus_max_ang_mom == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"ao_basis.nucleus_max_ang_mom",
NULL);
}
for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) {
ctx->ao_basis.nucleus_max_ang_mom[inucl] = 0;
for (int64_t ishell=ctx->ao_basis.nucleus_index[inucl] ;
ishell < ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl] ;
++ishell) {
ctx->ao_basis.nucleus_max_ang_mom[inucl] =
ctx->ao_basis.nucleus_max_ang_mom[inucl] =
ctx->ao_basis.nucleus_max_ang_mom[inucl] > ctx->ao_basis.shell_ang_mom[ishell] ?
ctx->ao_basis.nucleus_max_ang_mom[inucl] : ctx->ao_basis.shell_ang_mom[ishell] ;
}
@ -1259,7 +1259,7 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
iprim < ctx->ao_basis.shell_prim_index[ishell] + ctx->ao_basis.shell_prim_num[ishell] ;
++iprim) {
double range = 1./ctx->ao_basis.exponent[iprim];
ctx->ao_basis.nucleus_range[inucl] =
ctx->ao_basis.nucleus_range[inucl] =
ctx->ao_basis.nucleus_range[inucl] > range ?
ctx->ao_basis.nucleus_range[inucl] : range;
}
@ -1271,10 +1271,10 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
return QMCKL_SUCCESS;
}
#+end_src
** Fortran interfaces
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_ao_basis_type (context, t) &
bind(C)
@ -1593,15 +1593,15 @@ for (int64_t i=0 ; i < prim_num ; ++i) {
ao_num_test = qmckl_get_ao_basis_ao_num(context);
assert(ao_num == ao_num_test);
ao_factor_test = qmckl_get_ao_basis_ao_factor (context);
ao_factor_test = qmckl_get_ao_basis_ao_factor (context);
for (int64_t i=0 ; i < ao_num ; ++i) {
assert(ao_factor_test[i] == ao_factor[i]);
}
#+end_src
* Radial part
** TODO Helper functions to accelerate calculations
** General functions for Gaussian basis functions
@ -1745,9 +1745,12 @@ end function qmckl_ao_gaussian_vgl
#+begin_src f90 :tangle (eval f_test)
integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
use qmckl
use qmckl_probes_f
implicit none
integer(c_int64_t), intent(in), value :: context
logical(C_BOOL) :: vfc_err
integer*8 :: n, ldv, j, i
double precision :: X(3), R(3), Y(3), r2
@ -1756,6 +1759,13 @@ integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
epsilon = qmckl_get_numprec_epsilon(context)
#ifdef VFC_CI
! Multplying epsilon by 16 = 2^4 is equivalent to asking 4 significant digits
! less. This makes sense because we are adding noise with MCA so we can't be
! as strict on the accuracy target.
epsilon = epsilon * 16
#endif
X = (/ 1.1 , 2.2 , 3.3 /)
R = (/ 0.1 , 1.2 , -2.3 /)
Y(:) = X(:) - R(:)
@ -1772,10 +1782,29 @@ integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
test_qmckl_ao_gaussian_vgl = &
qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv)
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "gaussian_vgl_2_1"//C_NULL_CHAR, &
DBLE(VGL(2,1)), DBLE(0), DBLE(epsilon))
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "gaussian_vgl_2_2"//C_NULL_CHAR, &
DBLE(VGL(2,2)), DBLE(0), DBLE(epsilon))
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "gaussian_vgl_2_3"//C_NULL_CHAR, &
DBLE(VGL(2,3)), DBLE(0), DBLE(epsilon))
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "gaussian_vgl_2_4"//C_NULL_CHAR, &
DBLE(VGL(2,4)), DBLE(0), DBLE(epsilon))
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "gaussian_vgl_2_5"//C_NULL_CHAR, &
DBLE(VGL(2,5)), DBLE(0), DBLE(epsilon))
#ifndef VFC_CI
if (test_qmckl_ao_gaussian_vgl /= 0) return
#endif
test_qmckl_ao_gaussian_vgl = -1
#ifndef VFC_CI
do i=1,n
test_qmckl_ao_gaussian_vgl = -11
if (dabs(1.d0 - VGL(i,1) / (&
@ -1802,6 +1831,7 @@ integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) &
)) > epsilon ) return
end do
#endif
test_qmckl_ao_gaussian_vgl = 0
@ -1826,7 +1856,7 @@ qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double*
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double* const primitive_vgl) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -1869,7 +1899,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
"qmckl_ao_basis_primitive_vgl",
NULL);
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->ao_basis.primitive_vgl_date) {
@ -1890,7 +1920,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
ctx->ao_basis.primitive_vgl = primitive_vgl;
}
qmckl_exit_code rc;
qmckl_exit_code rc;
if (ctx->ao_basis.type == 'G') {
rc = qmckl_compute_ao_basis_primitive_gaussian_vgl(context,
ctx->ao_basis.prim_num,
@ -1906,7 +1936,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
QMCKL_FAILURE,
"compute_ao_basis_primitive_vgl",
"Not yet implemented");
}
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -1935,7 +1965,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
| double | nucl_coord[3][elec_num] | in | Nuclear coordinates |
| double | expo[prim_num] | in | Exponents of the primitives |
| double | primitive_vgl[5][elec_num][prim_num] | out | Value, gradients and Laplacian of the primitives |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
prim_num, elec_num, nucl_num, &
@ -1965,9 +1995,9 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
! C is zero-based, so shift bounds by one
do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1)
do ielec = 1, elec_num
x = elec_coord(ielec,1) - nucl_coord(inucl,1)
y = elec_coord(ielec,2) - nucl_coord(inucl,2)
z = elec_coord(ielec,3) - nucl_coord(inucl,3)
x = elec_coord(ielec,1) - nucl_coord(inucl,1)
y = elec_coord(ielec,2) - nucl_coord(inucl,2)
z = elec_coord(ielec,3) - nucl_coord(inucl,3)
r2 = x*x + y*y + z*z
ar2 = expo(iprim)*r2
@ -1977,7 +2007,7 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
two_a = -2.d0 * expo(iprim) * v
primitive_vgl(iprim, ielec, 1) = v
primitive_vgl(iprim, ielec, 2) = two_a * x
primitive_vgl(iprim, ielec, 2) = two_a * x
primitive_vgl(iprim, ielec, 3) = two_a * y
primitive_vgl(iprim, ielec, 4) = two_a * z
primitive_vgl(iprim, ielec, 5) = two_a * (3.d0 - 2.d0*ar2)
@ -2050,7 +2080,7 @@ qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl(
import numpy as np
def f(a,x,y):
return np.exp( -a*(np.linalg.norm(x-y))**2 )
return np.exp( -a*(np.linalg.norm(x-y))**2 )
def df(a,x,y,n):
h0 = 1.e-6
@ -2102,7 +2132,7 @@ int64_t elec_up_num = chbrclf_elec_up_num;
int64_t elec_dn_num = chbrclf_elec_dn_num;
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
@ -2110,7 +2140,7 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
assert(rc == QMCKL_SUCCESS);
@ -2125,7 +2155,7 @@ assert( fabs(prim_vgl[1][26][7] - (-7.5014974095310560E-004)) < 1.e-14 );
assert( fabs(prim_vgl[2][26][7] - (-3.8250692897610380E-003)) < 1.e-14 );
assert( fabs(prim_vgl[3][26][7] - ( 3.4950559194080275E-003)) < 1.e-14 );
assert( fabs(prim_vgl[4][26][7] - ( 2.0392163767356572E-002)) < 1.e-14 );
}
@ -2140,11 +2170,11 @@ assert( fabs(prim_vgl[4][26][7] - ( 2.0392163767356572E-002)) < 1.e-14 );
k=0;
for (j=0 ; j<elec_num ; ++j) {
for (i=0 ; i<nucl_num ; ++i) {
r2 = nucl_elec_dist[i][j];
if (r2 < nucl_radius2[i]) {
for (l=0 ; l<prim_num ; ++l) {
tmp[k].i = i;
tmp[k].j = j;
@ -2171,7 +2201,7 @@ qmckl_exit_code qmckl_get_ao_basis_shell_vgl(qmckl_context context, double* cons
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_ao_basis_shell_vgl(qmckl_context context, double* const shell_vgl) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -2235,7 +2265,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
"qmckl_electron",
NULL);
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->ao_basis.shell_vgl_date) {
@ -2255,7 +2285,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
ctx->ao_basis.shell_vgl = shell_vgl;
}
qmckl_exit_code rc;
qmckl_exit_code rc;
if (ctx->ao_basis.type == 'G') {
rc = qmckl_compute_ao_basis_shell_gaussian_vgl(context,
ctx->ao_basis.prim_num,
@ -2276,7 +2306,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
QMCKL_FAILURE,
"compute_ao_basis_shell_vgl",
"Not yet implemented");
}
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -2310,7 +2340,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
| ~double~ | ~expo[prim_num]~ | in | Exponents of the primitives |
| ~double~ | ~coef_normalized[prim_num]~ | in | Coefficients of the primitives |
| ~double~ | ~shell_vgl[5][elec_num][shell_num]~ | out | Value, gradients and Laplacian of the shells |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
prim_num, shell_num, elec_num, nucl_num, &
@ -2347,9 +2377,9 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
do ielec = 1, elec_num
x = elec_coord(ielec,1) - nucl_coord(inucl,1)
y = elec_coord(ielec,2) - nucl_coord(inucl,2)
z = elec_coord(ielec,3) - nucl_coord(inucl,3)
x = elec_coord(ielec,1) - nucl_coord(inucl,1)
y = elec_coord(ielec,2) - nucl_coord(inucl,2)
z = elec_coord(ielec,3) - nucl_coord(inucl,3)
r2 = x*x + y*y + z*z
@ -2376,13 +2406,13 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
shell_vgl(ishell, ielec, 1) + v
shell_vgl(ishell, ielec, 2) = &
shell_vgl(ishell, ielec, 2) + two_a * x
shell_vgl(ishell, ielec, 2) + two_a * x
shell_vgl(ishell, ielec, 3) = &
shell_vgl(ishell, ielec, 3) + two_a * y
shell_vgl(ishell, ielec, 3) + two_a * y
shell_vgl(ishell, ielec, 4) = &
shell_vgl(ishell, ielec, 4) + two_a * z
shell_vgl(ishell, ielec, 4) + two_a * z
shell_vgl(ishell, ielec, 5) = &
shell_vgl(ishell, ielec, 5) + two_a * (3.d0 - 2.d0*ar2)
@ -2415,7 +2445,7 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
const double* nucl_coord,
const double* expo,
const double* coef_normalized,
double* const shell_vgl );
double* const shell_vgl );
#+end_src
#+CALL: generate_c_interface(table=qmckl_ao_basis_shell_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_shell_gaussian_vgl"))
@ -2481,7 +2511,7 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
import numpy as np
def f(a,x,y):
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
def df(a,x,y,n):
h0 = 1.e-6
@ -2545,7 +2575,7 @@ int64_t elec_up_num = chbrclf_elec_up_num;
int64_t elec_dn_num = chbrclf_elec_dn_num;
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
@ -2553,7 +2583,7 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
assert(rc == QMCKL_SUCCESS);
@ -2567,15 +2597,15 @@ printf(" shell_vgl[1][1][26] %25.15e\n", shell_vgl[1][26][1]);
printf(" shell_vgl[1][2][26] %25.15e\n", shell_vgl[2][26][1]);
printf(" shell_vgl[1][3][26] %25.15e\n", shell_vgl[3][26][1]);
printf(" shell_vgl[1][4][26] %25.15e\n", shell_vgl[4][26][1]);
assert( fabs(shell_vgl[0][26][1] - ( 3.564393437193868e-02)) < 1.e-14 );
assert( fabs(shell_vgl[1][26][1] - (-6.030177987072189e-03)) < 1.e-14 );
assert( fabs(shell_vgl[2][26][1] - (-3.074832579537582e-02)) < 1.e-14 );
assert( fabs(shell_vgl[3][26][1] - ( 2.809546963519935e-02)) < 1.e-14 );
assert( fabs(shell_vgl[4][26][1] - ( 1.896046117183968e-02)) < 1.e-14 );
}
}
#+end_src
* Polynomial part
@ -2637,7 +2667,7 @@ assert( fabs(shell_vgl[4][26][1] - ( 1.896046117183968e-02)) < 1.e-14 );
const double* X,
const int32_t* LMAX,
double* const P,
const int64_t ldp );
const int64_t ldp );
#+end_src
*** Source
@ -2744,8 +2774,12 @@ end function qmckl_ao_power_f
#+begin_src f90 :tangle (eval f_test)
integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
use qmckl
use qmckl_probes_f
implicit none
logical(C_BOOL) :: vfc_err
integer(qmckl_context), intent(in), value :: context
integer*8 :: n, LDP
@ -2756,6 +2790,13 @@ integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
epsilon = qmckl_get_numprec_epsilon(context)
#ifdef VFC_CI
! Multplying epsilon by 16 = 2^4 is equivalent to asking 4 significant digits
! less. This makes sense because we are adding noise with MCA so we can't be
! as strict on the accuracy target.
epsilon = epsilon * 16
#endif
n = 100;
LDP = 10;
@ -2767,10 +2808,15 @@ integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
end do
test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP)
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "power_2_2"//C_NULL_CHAR, &
DBLE(P(2,2)), DBLE(0), DBLE(epsilon))
if (test_qmckl_ao_power /= QMCKL_SUCCESS) return
test_qmckl_ao_power = QMCKL_FAILURE
#ifndef VFC_CI
do j=1,n
do i=1,LMAX(j)
if ( X(j)**i == 0.d0 ) then
@ -2780,6 +2826,7 @@ integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
end if
end do
end do
#endif
test_qmckl_ao_power = QMCKL_SUCCESS
deallocate(X,P,LMAX)
@ -3072,9 +3119,12 @@ end function qmckl_ao_polynomial_vgl_f
#+begin_src f90 :tangle (eval f_test)
integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
use qmckl
use qmckl_probes_f
implicit none
integer(c_int64_t), intent(in), value :: context
logical(C_BOOL) :: vfc_err
integer :: lmax, d, i
integer, allocatable :: L(:,:)
@ -3101,9 +3151,25 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
test_qmckl_ao_polynomial_vgl = &
qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv)
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "polynomial_vgl_1_2"//C_NULL_CHAR, &
DBLE(VGL(1,2)), DBLE(0), DBLE(epsilon))
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "polynomial_vgl_2_2"//C_NULL_CHAR, &
DBLE(VGL(2,2)), DBLE(0), DBLE(epsilon))
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "polynomial_vgl_3_2"//C_NULL_CHAR, &
DBLE(VGL(3,2)), DBLE(0), DBLE(epsilon))
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "polynomial_vgl_4_2"//C_NULL_CHAR, &
DBLE(VGL(4,2)), DBLE(0), DBLE(epsilon))
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "polynomial_vgl_5_2"//C_NULL_CHAR, &
DBLE(VGL(5,2)), DBLE(0), DBLE(epsilon))
if (test_qmckl_ao_polynomial_vgl /= QMCKL_SUCCESS) return
if (n /= d) return
#ifndef VFC_CI
do j=1,n
test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
do i=1,3
@ -3154,6 +3220,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
end if
if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return
end do
#endif
test_qmckl_ao_polynomial_vgl = QMCKL_SUCCESS
@ -3176,7 +3243,7 @@ qmckl_exit_code qmckl_get_ao_vgl(qmckl_context context, double* const ao_vgl);
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_ao_vgl(qmckl_context context, double* const ao_vgl) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
@ -3240,7 +3307,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
"qmckl_electron",
NULL);
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->ao_basis.ao_vgl_date) {
@ -3317,7 +3384,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
| ~double~ | ~ao_factor[ao_num]~ | in | Normalization factor of the AOs |
| ~double~ | ~shell_vgl[5][elec_num][shell_num]~ | in | Value, gradients and Laplacian of the shells |
| ~double~ | ~ao_vgl[5][elec_num][ao_num]~ | out | Value, gradients and Laplacian of the AOs |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ao_vgl_f(context, &
ao_num, shell_num, elec_num, nucl_num, &
@ -3459,7 +3526,7 @@ end function qmckl_compute_ao_vgl_f
const int32_t* shell_ang_mom,
const double* ao_factor,
const double* shell_vgl,
double* const ao_vgl );
double* const ao_vgl );
#+end_src
#+CALL: generate_c_interface(table=qmckl_ao_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl"))
@ -3529,7 +3596,7 @@ import numpy as np
from math import sqrt
def f(a,x,y):
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
def df(a,x,y,n):
h0 = 1.e-6
@ -3567,7 +3634,7 @@ norm = sqrt(3.)
print ( "[0][26][219] : %25.15e"%(f(a,x,y) * (x[0] - y[0])**2) )
print ( "[1][26][219] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + 2.*f(a,x,y) * (x[0] - y[0])) )
print ( "[0][26][220] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1]) ))
print ( "[0][26][220] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1]) ))
print ( "[1][26][220] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + norm*f(a,x,y) * (x[1] - y[1])) )
print ( "[0][26][221] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) )
@ -3613,7 +3680,7 @@ int64_t elec_up_num = chbrclf_elec_up_num;
int64_t elec_dn_num = chbrclf_elec_dn_num;
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
@ -3621,7 +3688,7 @@ assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
assert(rc == QMCKL_SUCCESS);
@ -3657,8 +3724,8 @@ assert( fabs(ao_vgl[0][26][223] - (-4.021908374204471e-09)) < 1.e-14 );
assert( fabs(ao_vgl[1][26][223] - ( 2.154644255710413e-08)) < 1.e-14 );
assert( fabs(ao_vgl[0][26][224] - ( 7.175045873560788e-10)) < 1.e-14 );
assert( fabs(ao_vgl[1][26][224] - (-3.843864637762753e-09)) < 1.e-14 );
}
}
#+end_src
* End of files :noexport:
@ -3705,5 +3772,3 @@ assert( fabs(ao_vgl[1][26][224] - (-3.843864637762753e-09)) < 1.e-14 );
# -*- mode: org -*-
# vim: syntax=c

109
org/qmckl_distance.org
View File

@ -12,6 +12,7 @@ Functions for the computation of distances between particles.
#+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#include <stdio.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
@ -19,6 +20,7 @@ int main() {
qmckl_context context;
context = qmckl_context_create();
#+end_src
* Squared distance
@ -282,10 +284,17 @@ end function qmckl_distance_sq_f
*** Test :noexport:
#+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_distance_sq(context) bind(C)
use qmckl
use qmckl_probes_f
use iso_c_binding
implicit none
integer(qmckl_context), intent(in), value :: context
logical(C_BOOL) :: vfc_err
double precision, allocatable :: A(:,:), B(:,:), C(:,:)
integer*8 :: m, n, LDA, LDB, LDC
@ -299,7 +308,6 @@ integer(qmckl_exit_code) function test_qmckl_distance_sq(context) bind(C)
LDC = 5
allocate( A(LDA,m), B(LDB,n), C(LDC,n) )
do j=1,m
do i=1,m
A(i,j) = -10.d0 + dble(i+j)
@ -314,17 +322,26 @@ integer(qmckl_exit_code) function test_qmckl_distance_sq(context) bind(C)
test_qmckl_distance_sq = &
qmckl_distance_sq(context, 'X', 't', m, n, A, LDA, B, LDB, C, LDC)
vfc_err = qmckl_probe("distance"//C_NULL_CHAR, "distance_sq_Xt_2_2"//C_NULL_CHAR, DBLE(C(2,2)))
if (test_qmckl_distance_sq == 0) return
test_qmckl_distance_sq = &
qmckl_distance_sq(context, 't', 'X', m, n, A, LDA, B, LDB, C, LDC)
vfc_err = qmckl_probe("distance"//C_NULL_CHAR, "distance_sq_tX_2_2"//C_NULL_CHAR, DBLE(C(2,2)))
if (test_qmckl_distance_sq == 0) return
test_qmckl_distance_sq = &
qmckl_distance_sq(context, 'T', 't', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_distance_sq /= 0) return
vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_sq_Tt_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))
if (test_qmckl_distance_sq == 0) return
test_qmckl_distance_sq = QMCKL_FAILURE
@ -333,14 +350,17 @@ integer(qmckl_exit_code) function test_qmckl_distance_sq(context) bind(C)
x = (A(i,1)-B(j,1))**2 + &
(A(i,2)-B(j,2))**2 + &
(A(i,3)-B(j,3))**2
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
#ifndef VFC_CI
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14) return
#endif
end do
end do
test_qmckl_distance_sq = &
qmckl_distance_sq(context, 'n', 'T', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_distance_sq /= 0) return
vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_sq_nT_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))
test_qmckl_distance_sq = QMCKL_FAILURE
@ -349,14 +369,18 @@ integer(qmckl_exit_code) function test_qmckl_distance_sq(context) bind(C)
x = (A(1,i)-B(j,1))**2 + &
(A(2,i)-B(j,2))**2 + &
(A(3,i)-B(j,3))**2
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
#ifndef VFC_CI
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14) return
#endif
end do
end do
test_qmckl_distance_sq = &
qmckl_distance_sq(context, 'T', 'n', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_distance_sq /= 0) return
vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_sq_Tn_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))
if (test_qmckl_distance_sq == 0) return
test_qmckl_distance_sq = QMCKL_FAILURE
@ -365,14 +389,16 @@ integer(qmckl_exit_code) function test_qmckl_distance_sq(context) bind(C)
x = (A(i,1)-B(1,j))**2 + &
(A(i,2)-B(2,j))**2 + &
(A(i,3)-B(3,j))**2
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
#ifndef VFC_CI
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14) return
#endif
end do
end do
test_qmckl_distance_sq = &
qmckl_distance_sq(context, 'n', 'N', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_distance_sq /= 0) return
vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_sq_nN_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))
test_qmckl_distance_sq = QMCKL_FAILURE
@ -392,8 +418,8 @@ end function test_qmckl_distance_sq
#+end_src
#+begin_src c :comments link :tangle (eval c_test)
qmckl_exit_code test_qmckl_distance_sq(qmckl_context context);
assert(test_qmckl_distance_sq(context) == QMCKL_SUCCESS);
qmckl_exit_code test_qmckl_distance_sq(qmckl_context context);
assert(test_qmckl_distance_sq(context) == QMCKL_SUCCESS);
#+end_src
* Distance
@ -690,10 +716,17 @@ end function qmckl_distance_f
*** Test :noexport:
#+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_dist(context) bind(C)
use qmckl
use qmckl_probes_f
use iso_c_binding
implicit none
integer(qmckl_context), intent(in), value :: context
logical(C_BOOL) :: vfc_err
double precision, allocatable :: A(:,:), B(:,:), C(:,:)
integer*8 :: m, n, LDA, LDB, LDC
@ -722,17 +755,24 @@ integer(qmckl_exit_code) function test_qmckl_dist(context) bind(C)
test_qmckl_dist = &
qmckl_distance(context, 'X', 't', m, n, A, LDA, B, LDB, C, LDC)
vfc_err = qmckl_probe("distance"//C_NULL_CHAR, "distance_Xt_2_2"//C_NULL_CHAR, DBLE(C(2,2)))
if (test_qmckl_dist == 0) return
test_qmckl_dist = &
qmckl_distance(context, 't', 'X', m, n, A, LDA, B, LDB, C, LDC)
vfc_err = qmckl_probe("distance"//C_NULL_CHAR, "distance_tX_2_2"//C_NULL_CHAR, DBLE(C(2,2)))
if (test_qmckl_dist == 0) return
test_qmckl_dist = &
qmckl_distance(context, 'T', 't', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_dist /= 0) return
vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_Tt_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))
if (test_qmckl_dist == 0) return
test_qmckl_dist = QMCKL_FAILURE
@ -741,14 +781,19 @@ integer(qmckl_exit_code) function test_qmckl_dist(context) bind(C)
x = dsqrt((A(i,1)-B(j,1))**2 + &
(A(i,2)-B(j,2))**2 + &
(A(i,3)-B(j,3))**2)
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
#ifndef VFC_CI
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14) return
#endif
end do
end do
test_qmckl_dist = &
qmckl_distance(context, 'n', 'T', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_dist /= 0) return
vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_nT_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))
if (test_qmckl_dist == 0) return
test_qmckl_dist = QMCKL_FAILURE
@ -757,14 +802,19 @@ integer(qmckl_exit_code) function test_qmckl_dist(context) bind(C)
x = dsqrt((A(1,i)-B(j,1))**2 + &
(A(2,i)-B(j,2))**2 + &
(A(3,i)-B(j,3))**2)
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
#ifndef VFC_CI
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14) return
#endif
end do
end do
test_qmckl_dist = &
qmckl_distance(context, 'T', 'n', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_dist /= 0) return
vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_Tn_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))
if (test_qmckl_dist == 0) return
test_qmckl_dist = QMCKL_FAILURE
@ -773,14 +823,19 @@ integer(qmckl_exit_code) function test_qmckl_dist(context) bind(C)
x = dsqrt((A(i,1)-B(1,j))**2 + &
(A(i,2)-B(2,j))**2 + &
(A(i,3)-B(3,j))**2)
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
#ifndef VFC_CI
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14) return
#endif
end do
end do
test_qmckl_dist = &
qmckl_distance(context, 'n', 'N', m, n, A, LDA, B, LDB, C, LDC)
if (test_qmckl_dist /= 0) return
vfc_err = qmckl_probe_check("distance"//C_NULL_CHAR, "distance_nN_2_2"//C_NULL_CHAR, DBLE(C(2,2)), DBLE(0), DBLE(1.d-14))
if (test_qmckl_dist == 0) return
test_qmckl_dist = QMCKL_FAILURE
@ -789,7 +844,9 @@ integer(qmckl_exit_code) function test_qmckl_dist(context) bind(C)
x = dsqrt((A(1,i)-B(1,j))**2 + &
(A(2,i)-B(2,j))**2 + &
(A(3,i)-B(3,j))**2)
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
#ifndef VFC_CI
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14) return
#endif
end do
end do
@ -800,8 +857,8 @@ end function test_qmckl_dist
#+end_src
#+begin_src c :comments link :tangle (eval c_test)
qmckl_exit_code test_qmckl_dist(qmckl_context context);
assert(test_qmckl_dist(context) == QMCKL_SUCCESS);
qmckl_exit_code test_qmckl_dist(qmckl_context context);
assert(test_qmckl_dist(context) == QMCKL_SUCCESS);
#+end_src
* Rescaled Distance
@ -1114,12 +1171,12 @@ end function qmckl_distance_rescaled_f
:FRetType: qmckl_exit_code
:END:
~qmckl_distance_rescaled_deriv_e~ computes the matrix of the gradient and laplacian of the
~qmckl_distance_rescaled_deriv_e~ computes the matrix of the gradient and laplacian of the
rescaled distance with respect to the electron coordinates. The derivative is a rank 3 tensor.
The first dimension has a dimension of 4 of which the first three coordinates
contains the gradient vector and the last index is the laplacian.
\[
C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa
\]
@ -1130,12 +1187,12 @@ end function qmckl_distance_rescaled_f
\nabla (C_{ij}(\mathbf{r}_{ee})) = \left(\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta x},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta y},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta z} \right)
\]
and the laplacian is defined as follows:
\[
\triangle (C_{ij}(r_{ee})) = \frac{\delta^2}{\delta x^2} + \frac{\delta^2}{\delta y^2} + \frac{\delta^2}{\delta z^2}
\]
Using the above three formulae, the expression for the gradient and laplacian is
Using the above three formulae, the expression for the gradient and laplacian is
as follows:
\[
@ -1462,7 +1519,9 @@ end function qmckl_distance_rescaled_deriv_e_f
* End of files :noexport:
#+begin_src c :comments link :tangle (eval c_test)
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0;
}

139
org/qmckl_sherman_morrison_woodbury.org
View File

@ -27,7 +27,7 @@ int main() {
#+end_src
* Naïve Sherman-Morrison
** ~qmckl_sherman_morrison~
:PROPERTIES:
:Name: qmckl_sherman_morrison
@ -60,6 +60,9 @@ int main() {
the update is applied as usual and the kernel exits with return code \texttt{QMCKL_SUCCESS}.
If $1+v_j^TS^{-1}u_j \leq \epsilon$ the update is rejected and the kernel exits with return code \texttt{QMCKL_FAILURE}.
If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting
from applying the updates to the original matrix.
#+NAME: qmckl_sherman_morrison_args
| qmckl_context | context | in | Global state |
| uint64_t | Dim | in | Leading dimension of Slater_inv |
@ -68,6 +71,7 @@ int main() {
| uint64_t | Updates_index[N_updates] | in | Array containing the rank-1 updates |
| double | breakdown | in | Break-down parameter on which to fail or not |
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
| double* | determinant | inout | Determinant of the Slater-matrix |
*** Requirements
@ -92,7 +96,8 @@ int main() {
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv);
double* Slater_inv,
double* determinant);
#+end_src
*** C source
@ -108,7 +113,8 @@ qmckl_exit_code qmckl_sherman_morrison(const qmckl_context context,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv) {
double* Slater_inv,
double* determinant) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -130,11 +136,16 @@ qmckl_exit_code qmckl_sherman_morrison(const qmckl_context context,
// Denominator
double den = 1 + C[Updates_index[l] - 1];
if (fabs(den) < breakdown) {
return QMCKL_FAILURE;
}
double iden = 1 / den;
// Update det(A)
if (determinant != NULL)
*determinant *= den;
// D = v^T x A^{-1}
for (uint64_t j = 0; j < Dim; j++) {
D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
@ -174,7 +185,7 @@ qmckl_exit_code qmckl_sherman_morrison(const qmckl_context context,
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_sherman_morrison &
(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv) &
(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
@ -188,6 +199,7 @@ qmckl_exit_code qmckl_sherman_morrison(const qmckl_context context,
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
real (c_double ) , intent(inout) :: determinant
end function qmckl_sherman_morrison
end interface
@ -234,7 +246,14 @@ double Slater_inv5[441] = {-0.054189244668834902, -105.426713929607, -88.4584964
assert(Updates1 != NULL);
assert(Updates_index1 != NULL);
assert(Slater_inv1 != NULL);
rc = qmckl_sherman_morrison(context, Dim, N_updates1, Updates1, Updates_index1, breakdown, Slater_inv1);
// original determinant of Slater1 (before applying updates)
double det = 3.407025646103221e-10;
rc = qmckl_sherman_morrison(context, Dim, N_updates1, Updates1, Updates_index1, breakdown, Slater_inv1, &det);
// Check that the determinant is updated properly
assert(fabs(det + 4.120398385068217e-10) < 1e-15);
for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) {
res[i * Dim + j] = 0;
@ -279,6 +298,10 @@ assert(rc == QMCKL_SUCCESS);
$B := 1 + VC$, the $2 \times 2$ matrix that is going to be inverted
$D := VS^{-1}$, a $2 \times Dim$ matrix
If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting
from applying the updates to the original matrix.
#+NAME: qmckl_woodbury_2_args
| qmckl_context | context | in | Global state |
@ -287,6 +310,7 @@ assert(rc == QMCKL_SUCCESS);
| uint64_t | Updates_index[2] | in | Array containing the rank-1 updates |
| double | breakdown | in | Break-down parameter on which to fail or not |
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
| double* | determinant | inout | Determinant of Slater-matrix |
*** Requirements
@ -309,7 +333,8 @@ assert(rc == QMCKL_SUCCESS);
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv);
double* Slater_inv,
double* determinant);
#+end_src
*** C source
@ -324,7 +349,8 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv) {
double* Slater_inv,
double* determinant) {
/*
C := S^{-1} * U, dim x 2
B := 1 + V * C, 2 x 2
@ -362,6 +388,10 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context,
return QMCKL_FAILURE;
}
// Update det(S) when passed
if (determinant != NULL)
*determinant *= det;
// Compute B^{-1} with explicit formula for 2x2 inversion
double Binv[4], idet = 1.0 / det;
Binv[0] = idet * B3;
@ -408,7 +438,7 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context,
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_woodbury_2 &
(context, Dim, Updates, Updates_index, breakdown, Slater_inv) &
(context, Dim, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
@ -420,6 +450,7 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context,
integer (c_int64_t) , intent(in) :: Updates_index(2)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
real (c_double ) , intent(inout) :: determinant
end function qmckl_woodbury_2
end interface
@ -431,7 +462,9 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context,
assert(Updates2 != NULL);
assert(Updates_index2 != NULL);
assert(Slater_inv2 != NULL);
rc = qmckl_woodbury_2(context, Dim, Updates2, Updates_index2, breakdown, Slater_inv2);
det = -1.4432116661319376e-11;
rc = qmckl_woodbury_2(context, Dim, Updates2, Updates_index2, breakdown, Slater_inv2, &det);
assert(fabs(det-2.367058141251457e-10) < 1e-15);
for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) {
res[i * Dim + j] = 0;
@ -471,6 +504,11 @@ assert(rc == QMCKL_SUCCESS);
$B := 1 + VC$, the $3 \times 3$ matrix that is going to be inverted
$D := VS^{-1}$, a $3 \times Dim$ matrix
If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting
from applying the updates to the original matrix.
#+NAME: qmckl_woodbury_3_args
| qmckl_context | context | in | Global state |
| uint64_t | Dim | in | Leading dimension of Slater_inv |
@ -478,6 +516,7 @@ assert(rc == QMCKL_SUCCESS);
| uint64_t | Updates_index[3] | in | Array containing the rank-1 updates |
| double | breakdown | in | Break-down parameter on which to fail or not |
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
| double* | determinant | inout | Determinant of Slater-matrix |
*** Requirements
@ -500,7 +539,8 @@ assert(rc == QMCKL_SUCCESS);
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv);
double* Slater_inv,
double* determinant);
#+end_src
*** C source
@ -515,7 +555,8 @@ qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv) {
double* Slater_inv,
double* determinant) {
/*
C := S^{-1} * U, dim x 3
B := 1 + V * C, 3 x 3
@ -561,6 +602,10 @@ qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
return QMCKL_FAILURE;
}
// Update det(Slater) if passed
if (determinant != NULL)
*determinant *= det;
// Compute B^{-1} with explicit formula for 3x3 inversion
double Binv[9], idet = 1.0 / det;
Binv[0] = (B4 * B8 - B7 * B5) * idet;
@ -614,7 +659,7 @@ qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_woodbury_3 &
(context, Dim, Updates, Updates_index, breakdown, Slater_inv) &
(context, Dim, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
@ -626,6 +671,7 @@ qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
integer (c_int64_t) , intent(in) :: Updates_index(3)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
real (c_double ) , intent(inout) :: determinant
end function qmckl_woodbury_3
end interface
@ -637,7 +683,9 @@ qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
assert(Updates3 != NULL);
assert(Updates_index3 != NULL);
assert(Slater_inv3_1 != NULL);
rc = qmckl_woodbury_3(context, Dim, Updates3, Updates_index3, breakdown, Slater_inv3_1);
det = -1.23743195512859e-09;
rc = qmckl_woodbury_3(context, Dim, Updates3, Updates_index3, breakdown, Slater_inv3_1, &det);
assert(fabs(det - 1.602708950725074e-10) < 1e-15);
for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) {
res[i * Dim + j] = 0;
@ -680,6 +728,9 @@ assert(rc == QMCKL_SUCCESS);
case the Slater-matrix that would have resulted from applying the updates is singular and therefore the
kernel exits with an exit code.
If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting
from applying the updates to the original matrix.
#+NAME: qmckl_sherman_morrison_splitting_args
| qmckl_context | context | in | Global state |
| uint64_t | Dim | in | Leading dimension of Slater_inv |
@ -688,6 +739,10 @@ assert(rc == QMCKL_SUCCESS);
| uint64_t | Updates_index[N_updates] | in | Array containing the rank-1 updates |
| double | breakdown | in | Break-down parameter on which to fail or not |
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
| double* | determinant | inout | Determinant of the Slater-matrix |
If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting
from applying the updates to the original matrix.
*** Requirements
@ -712,7 +767,8 @@ assert(rc == QMCKL_SUCCESS);
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv);
double* Slater_inv,
double* determinant);
#+end_src
*** C source
@ -727,7 +783,8 @@ qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv) {
double* Slater_inv,
double* determinant) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -738,11 +795,11 @@ qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context,
uint64_t later = 0;
(void) qmckl_slagel_splitting(Dim, N_updates, Updates, Updates_index,
breakdown, Slater_inv, later_updates, later_index, &later);
breakdown, Slater_inv, later_updates, later_index, &later, determinant);
if (later > 0) {
(void) qmckl_sherman_morrison_splitting(context, Dim, later,
later_updates, later_index, breakdown, Slater_inv);
later_updates, later_index, breakdown, Slater_inv, determinant);
}
return QMCKL_SUCCESS;
@ -767,7 +824,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context,
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_sherman_morrison_splitting &
(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv) &
(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
@ -780,6 +837,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context,
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
real (c_double ) , intent(inout) :: determinant
end function qmckl_sherman_morrison_splitting
end interface
@ -791,7 +849,9 @@ qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context,
assert(Updates3 != NULL);
assert(Updates_index3 != NULL);
assert(Slater_inv3_2 != NULL);
rc = qmckl_sherman_morrison_splitting(context, Dim, N_updates3, Updates3, Updates_index3, breakdown, Slater_inv3_2);
det = -1.23743195512859e-09;
rc = qmckl_sherman_morrison_splitting(context, Dim, N_updates3, Updates3, Updates_index3, breakdown, Slater_inv3_2, &det);
assert(fabs(det - 1.602708950725074e-10) < 1e-15);
for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) {
res[i * Dim + j] = 0;
@ -829,6 +889,9 @@ assert(rc == QMCKL_SUCCESS);
with Woodbury 2x2 instead of sending them all to Sherman-Morrison with update splitting. For example, in the case of
5 updates the updates are applied in 1 block of 3 updates end 1 block of 2 updates.
If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting
from applying the updates to the original matrix.
#+NAME: qmckl_sherman_morrison_smw32s_args
| qmckl_context | context | in | Global state |
| uint64_t | Dim | in | Leading dimension of Slater_inv |
@ -837,6 +900,8 @@ assert(rc == QMCKL_SUCCESS);
| uint64_t | Updates_index[N_updates] | in | Array containing the rank-1 updates |
| double | breakdown | in | Break-down parameter on which to fail or not |
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
| double* | determinant | inout | Determinant of the Slater-matrix |
*** Requirements
@ -861,7 +926,8 @@ assert(rc == QMCKL_SUCCESS);
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv);
double* Slater_inv,
double* determinant);
#+end_src
*** C source
@ -876,7 +942,8 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv) {
double* Slater_inv,
double* determinant) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -897,11 +964,11 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
for (uint64_t i = 0; i < n_of_3blocks; i++) {
const double *Updates_3block = &Updates[i * length_3block];
const uint64_t *Updates_index_3block = &Updates_index[i * 3];
rc = qmckl_woodbury_3(context, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv);
rc = qmckl_woodbury_3(context, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv, determinant);
if (rc != 0) { // Send the entire block to slagel_splitting
uint64_t l = 0;
rc = qmckl_slagel_splitting(Dim, 3, Updates_3block, Updates_index_3block,
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l);
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant);
later = later + l;
}
}
@ -911,11 +978,11 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
if (remainder == 2) {
const double *Updates_2block = &Updates[n_of_3blocks * length_3block];
const uint64_t *Updates_index_2block = &Updates_index[3 * n_of_3blocks];
rc = qmckl_woodbury_2(context, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv);
rc = qmckl_woodbury_2(context, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv, determinant);
if (rc != 0) { // Send the entire block to slagel_splitting
uint64_t l = 0;
(void) qmckl_slagel_splitting(Dim, 2, Updates_2block, Updates_index_2block,
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l);
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant);
later = later + l;
}
}
@ -925,12 +992,12 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
const uint64_t *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
uint64_t l = 0;
(void) qmckl_slagel_splitting(Dim, 1, Updates_1block, Updates_index_1block,
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l);
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant);
later = later + l;
}
if (later > 0) {
(void) qmckl_sherman_morrison_splitting(context, Dim, later, later_updates, later_index, breakdown, Slater_inv);
(void) qmckl_sherman_morrison_splitting(context, Dim, later, later_updates, later_index, breakdown, Slater_inv, determinant);
}
return QMCKL_SUCCESS;
}
@ -978,7 +1045,10 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
assert(Updates5 != NULL);
assert(Updates_index5 != NULL);
assert(Slater_inv5 != NULL);
rc = qmckl_sherman_morrison_smw32s(context, Dim, N_updates5, Updates5, Updates_index5, breakdown, Slater_inv5);
det = -3.186005284713128e-10;
rc = qmckl_sherman_morrison_smw32s(context, Dim, N_updates5, Updates5, Updates_index5, breakdown, Slater_inv5, &det);
assert(fabs(det + 5.260200118412903e-10) < 1e-15);
for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) {
res[i * Dim + j] = 0;
@ -1022,6 +1092,9 @@ These functions can only be used internally by the kernels in this module.
as $\frac{1}{2}u_j$. One half is applied immediately, the other half will be applied at the end of the algorithm, using vectors
$u_{j'}=\frac{1}{2}u_j$ and $v_{j'}^T=v_{j}^T$, which are stored in the array \texttt{later_updates}.
If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting
from applying the updates to the original matrix.
#+NAME: qmckl_slagel_splitting_args
| uint64_t | Dim | in | Leading dimension of Slater_inv |
| uint64_t | N_updates | in | Number of rank-1 updates to be applied to Slater_inv |
@ -1032,6 +1105,7 @@ These functions can only be used internally by the kernels in this module.
| double | later_updates[Dim * N_updates] | inout | Array containing the split updates for later |
| uint64_t | later_index[N_updates] | inout | Array containing the positions of the split updates for later |
| uint64_t | later | inout | Number of split updates for later |
| double* | determinant | inout | Determinant of the Slater-matrix |
*** Requirements
@ -1061,7 +1135,8 @@ These functions can only be used internally by the kernels in this module.
double* Slater_inv,
double* later_updates,
uint64_t* later_index,
uint64_t* later );
uint64_t* later,
double* determinant);
#+end_src
*** C source
@ -1079,7 +1154,8 @@ qmckl_exit_code qmckl_slagel_splitting(uint64_t Dim,
double *Slater_inv,
double *later_updates,
uint64_t *later_index,
uint64_t *later) {
uint64_t *later,
double *determinant) {
// #ifdef DEBUG // Leave commented out since debugging information is not yet implemented in QMCkl.
// std::cerr << "Called slagel_splitting with " << N_updates << " updates" << std::endl;
// #endif
@ -1114,6 +1190,9 @@ qmckl_exit_code qmckl_slagel_splitting(uint64_t Dim,
} // From here onwards we continue with applying the first havel of the update to Slater_inv
double iden = 1 / den;
if (determinant != NULL)
*determinant *= den;
// D = v^T x S^{-1}
for (uint64_t j = 0; j < Dim; j++) {
D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];

141
tools/qmckl_probes.c Normal file
View File

@ -0,0 +1,141 @@
#include <stdbool.h>
#include <stdlib.h>
#include <stdio.h>
#ifdef VFC_CI
#include <vfc_probes.h>
vfc_probes probes;
#endif
// QMCkl is a wrapper to Verificarlo's vfc_probes system. The goal of QMCkl
// probes isto simplify the use of vfc_probes, and to provide functions that
// can be called either wit or without vfc_ci support by using #ifndef
// statements :
//
// - when vfc_ci is disabled, qmckl_probes functions will either return false
// (no error) or perform a check based on a reference value
// - when vfc_ci is enabled, qmckl_probe functions will simply encapsulate
// calls to vfc_probe
//
// Moreover, one does not have to worry about the life cycle of the probes
// structure, as it is automatically created, dumped and freed by this wrapper.
//
// vfc_ci support can be enabled by using the following configure command :
// QMCKL_DEVEL=1 ./configure --prefix=$PWD/_install --enable-silent-rules
// --enable-maintainer-mode CC=verificarlo-f FC=verificarlo-f --host=x86_64
//
// Finally, this wrapper also comes with a Fortran interface (in its dedicated
// file).
//
// To learn more about Verificarlo CI :
// https://github.com/verificarlo/verificarlo/blob/master/doc/06-Postprocessing.md#verificarlo-ci
// Automatically initialize the vfc_probe object if VFC_CI is defined
#ifdef VFC_CI
void __attribute__((constructor)) qmckl_init_probes(){
probes = vfc_init_probes();
}
#endif
// Standard probe, without check
// - if VFC_CI is defined, place a standard probe
// - if VFC_CI is undefined, return false (no error)
bool qmckl_probe(
char * testName,
char * varName,
double value
) {
#ifdef VFC_CI
return vfc_probe(&probes, testName, varName, value);
#else
return false;
#endif
}
// Probe with absolute check
// - if VFC_CI is defined, place a probe with an absolute check
// - if VFC_CI is undefined, perform an absolute check based on target value
// and accuracy
bool qmckl_probe_check(
char * testName,
char * varName,
double value,
double expectedValue,
double accuracyTarget
) {
#ifdef VFC_CI
return vfc_probe_check(&probes, testName, varName, value, accuracyTarget);
#else
return !(abs(value - expectedValue) < accuracyTarget);
#endif
}
// Probe with relative check
// - if VFC_CI is defined, place a probe with a relative check
// - if VFC_CI is undefined, perform a relative check based on target value
// and accuracy
bool qmckl_probe_check_relative (
char * testName,
char * varName,
double value,
double expectedValue,
double accuracyTarget
) {
#ifdef VFC_CI
return vfc_probe_check_relative(&probes, testName, varName, value, accuracyTarget);
#else
return !(abs(value - expectedValue) / abs(expectedValue) < accuracyTarget);
#endif
}
// Automatically delete and dump the vfc_probe object if VFC_CI is defined
#ifdef VFC_CI
void __attribute__((destructor)) qmckl_dump_probes(){
vfc_dump_probes(&probes);
}
#endif
// Fortran wrappers
bool qmckl_probe_f(
char * testName,
char * varName,
double * value
) {
return qmckl_probe(testName, varName, *value);
}
bool qmckl_probe_check_f(
char * testName,
char * varName,
double * value,
double * expectedValue,
double * accuracyTarget
) {
return qmckl_probe_check(
testName, varName,
*value, *expectedValue, *accuracyTarget
);
}
bool qmckl_probe_check_relative_f(
char * testName,
char * varName,
double * value,
double * expectedValue,
double * accuracyTarget
) {
return qmckl_probe_check_relative(
testName, varName,
*value, *expectedValue, *accuracyTarget
);
}

65
tools/qmckl_probes.h Normal file
View File

@ -0,0 +1,65 @@
#include <stdbool.h>
#ifdef VFC_CI
#include <vfc_probes.h>
extern vfc_probes * probes;
#endif
// Wrappers to Verificarlo functions
#ifdef VFC_CI
void qmckl_init_probes() __attribute__((constructor));
#endif
bool qmckl_probe(
char * testName,
char * varName,
double value
);
bool qmckl_probe_check(
char * testName,
char * varName,
double value,
double expectedValue,
double accuracyTarget
);
bool qmckl_probe_check_relative(
char * testName,
char * varName,
double value,
double expectedValue,
double accuracyTarget
);
#ifdef VFC_CI
void qmckl_dump_probes() __attribute__((destructor));
#endif
// Fortran wrappers
bool qmckl_probe_f(
char * testName,
char * varName,
double * value
);
bool qmckl_probe_check_f(
char * testName,
char * varName,
double * value,
double * expectedValue,
double * accuracyTarget
);
bool qmckl_probe_check_relative_f(
char * testName,
char * varName,
double * value,
double * expectedValue,
double * accuracyTarget
);

49
tools/qmckl_probes_f.f90 Normal file
View File

@ -0,0 +1,49 @@
module qmckl_probes_f
interface
logical(c_bool) function qmckl_probe &
(testName, varName, val) &
bind(C, name="qmckl_probe_f")
use, intrinsic :: iso_c_binding
import
implicit none
character(C_CHAR), dimension(*) :: testName
character(C_CHAR), dimension(*) :: varName
real(C_DOUBLE) :: val
end function qmckl_probe
logical(c_bool) function qmckl_probe_check &
(testName, varName, val, expectedValue, accuracyTarget) &
bind(C, name="qmckl_probe_check_f")
use, intrinsic :: iso_c_binding
import
implicit none
character(C_CHAR), dimension(*) :: testName
character(C_CHAR), dimension(*) :: varName
real(C_DOUBLE) :: val
real(C_DOUBLE) :: expectedValue
real(C_DOUBLE) :: accuracyTarget
end function qmckl_probe_check
logical(c_bool) function qmckl_probe_check_relative &
(testName, varName, val, expectedValue, accuracyTarget) &
bind(C, name="qmckl_probe_check_relative_f")
use, intrinsic :: iso_c_binding
import
implicit none
character(C_CHAR), dimension(*) :: testName
character(C_CHAR), dimension(*) :: varName
real(C_DOUBLE) :: val
real(C_DOUBLE) :: expectedValue
real(C_DOUBLE) :: accuracyTarget
end function qmckl_probe_check_relative
end interface
end module qmckl_probes_f

19
vfc_tests_config.json Normal file
View File

@ -0,0 +1,19 @@
{
"make_command": "./ci_install.sh",
"executables": [
{
"executable": "tests/test_qmckl_distance",
"vfc_backends": [{
"name": "libinterflop_mca.so",
"repetitions": 25
}]
},
{
"executable": "tests/test_qmckl_ao",
"vfc_backends": [{
"name": "libinterflop_mca.so",
"repetitions": 25
}]
}
]
}