diff --git a/.gitignore b/.gitignore index 44474a6..0a4ba27 100644 --- a/.gitignore +++ b/.gitignore @@ -8,4 +8,4 @@ make.config src/ZMQ/f77_zmq.h ezfio_config/properties.config bin/ - +*.mod diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..ec36f9f --- /dev/null +++ b/.gitmodules @@ -0,0 +1,6 @@ +[submodule "EZFIO"] + path = EZFIO + url = https://gitlab.com/scemama/EZFIO.git +[submodule "f77_zmq"] + path = f77_zmq + url = https://github.com/zeromq/f77_zmq.git diff --git a/EZFIO b/EZFIO new file mode 160000 index 0000000..ed1df9f --- /dev/null +++ b/EZFIO @@ -0,0 +1 @@ +Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c diff --git a/Makefile b/Makefile index f988cee..aad6084 100644 --- a/Makefile +++ b/Makefile @@ -1,3 +1,14 @@ +# QMC=Chem + +FC = gfortran +LIBS = +CPPFLAGS = +LDFLAGS= + +ZMQ_LIBS = -lzmq +ZMQ_LDFLAGS = +ZMQ_CPPFLAGS = + include make.config diff --git a/Makefile.am b/Makefile.am new file mode 100644 index 0000000..6a452da --- /dev/null +++ b/Makefile.am @@ -0,0 +1,76 @@ +# QMC=Chem + +FC = @FC@ +LIBS = @LIBS@ +CPPFLAGS = @CPPFLAGS@ +LDFLAGS= @LDFLAGS@ + +ZMQ_LIBS = @ZMQ_LIBS@ +ZMQ_LDFLAGS = @ZMQ_LDFLAGS@ +ZMQ_CPPFLAGS = @ZMQ_CPPFLAGS@ + + + +default_target: all + +.phony: clean always + +clean: + ./scripts/clean.sh + +# put all files of PROPERTIES here +# -------------------------------- + +ezfio_config/properties.config ocaml/Property.ml: scripts/create_properties_python.py src/PROPERTIES/properties.irp.f src/PROPERTIES/properties_energy.irp.f src/PROPERTIES/properties_general.irp.f src/irpf90.make + bash -c "source qmcchemrc ; exec python2 ./scripts/create_properties_ezfio.py" + +# EZFIO +# ----- + +build lib/libezfio.a lib/libezfio_irp.a EZFIO/lib/libezfio.a EZFIO/lib/libezfio_irp.a EZFIO/Ocaml/ezfio.ml EZFIO/Python/ezfio.py: ezfio_config/qmc.config ezfio_config/properties.config make.config scripts/create_properties_ezfio.py src/tags src/irpf90_entities src/irpf90.make src/irpf90.make + ./scripts/compile_ezfio.sh + + +# Fortran executables +# ------------------- + + +always: /dev/null + +src/tags src/irpf90_entities src/irpf90.make: make.config always + ./scripts/compile_irpf90.sh + +src/MAIN/qmc src/MAIN/vmc_test src/MAIN/qmc_create_walkers src/MAIN/qmcchem_info: lib/libezfio.a lib/libezfio_irp.a src/tags src/irpf90_entities src/irpf90.make src/irpf90.make + ./scripts/compile_src.sh + +# OCaml +# ----- + +ocaml/qmcchem: EZFIO/Ocaml/ezfio.ml $(wildcard ocaml/*.ml) + ./scripts/compile_ocaml.sh + +# Archive +# ------- + +qmcchem.tar.gz: all + git archive --format tar.gz HEAD --prefix "QmcChem/" -7 -o qmcchem.tar.gz + + +# Binaries +# -------- + +bin/qmc: src/MAIN/qmc + cp src/MAIN/qmc bin/ + +bin/qmcchem_info: src/MAIN/qmcchem_info + cp src/MAIN/qmcchem_info bin/ + +bin/qmc_create_walkers: src/MAIN/qmc_create_walkers + cp src/MAIN/qmc_create_walkers bin/ + +bin/qmcchem: ocaml/qmcchem + cp ocaml/qmcchem bin/ + +all: bin/qmc bin/qmcchem_info bin/qmc_create_walkers bin/qmcchem + + diff --git a/configure.ac b/configure.ac new file mode 100644 index 0000000..7f42e88 --- /dev/null +++ b/configure.ac @@ -0,0 +1,57 @@ +# -*- Autoconf -*- +# Process this file with autoconf to produce a configure script. + +AC_PREREQ([2.69]) + +AC_INIT([QMC=Chem], [2.0.0], [https://gitlab.com/scemama/qmcchem/-/issues/new]) +AM_INIT_AUTOMAKE([foreign subdir-objects silent-rules]) + +AC_CONFIG_SRCDIR([README.md]) +AC_CONFIG_MACRO_DIR([m4]) +AC_CONFIG_HEADERS([include/config.h]) + +# Checks for programs. +AC_PROG_AWK +AC_PROG_CC +AC_PROG_INSTALL +AC_PROG_LN_S +AC_PROG_MAKE_SET +AX_PROG_IRPF90 + +# Fortran compiler checks +AC_PROG_FC([ifort gfortran]) +AC_FC_LINE_LENGTH([512]) + + +# Checks for libraries. +AX_ZMQ([4.0], [], [ AC_MSG_ERROR([Please install ZeroMQ >= 4.0]) ]) +AX_BLAS() +AX_LAPACK() + +# Required by EZFIO +AC_CHECK_FUNCS([mkdir strerror]) +AC_TYPE_SIZE_T + +# Required by OCaml C bindings +AC_CHECK_FUNCS([inet_ntoa]) +AC_CHECK_FUNCS([memset]) +AC_CHECK_FUNCS([socket]) +AC_CHECK_HEADERS([arpa/inet.h netinet/in.h sys/ioctl.h sys/socket.h]) + +# Required by ZeroMQ +AC_CHECK_HEADERS([stddef.h]) +AC_TYPE_INT32_T +AC_TYPE_UINT16_T +AC_TYPE_UINT32_T +AC_TYPE_UINT8_T + +# Required by QMCkl +#AC_CHECK_HEADER_STDBOOL +#AC_TYPE_INT64_T +#AC_TYPE_UINT64_T + + +AC_LANG([Fortran]) + +AC_CONFIG_FILES([Makefile]) +AC_OUTPUT diff --git a/configure.sh b/configure.sh index 1da73d7..9c36e9a 100755 --- a/configure.sh +++ b/configure.sh @@ -41,6 +41,7 @@ export LD_LIBRARY_PATH="\${QMCCHEM_PATH}/lib:\${LD_LIBRARY_PATH}" export LIBRARY_PATH="\${QMCCHEM_PATH}/lib:\${LIBRARY_PATH}" export QMCCHEM_MPIRUN="mpirun" export QMCCHEM_MPIRUN_FLAGS="" +#export QMCCHEM_IO="B" export C_INCLUDE_PATH="\${QMCCHEM_PATH}/include:\${C_INCLUDE_PATH}" #export QMCCHEM_NIC=ib0 source \${QMCCHEM_PATH}/irpf90/bin/irpman diff --git a/f77_zmq b/f77_zmq new file mode 160000 index 0000000..934e063 --- /dev/null +++ b/f77_zmq @@ -0,0 +1 @@ +Subproject commit 934e063881553e42d81c5d1aaed35822988529f0 diff --git a/m4/ax_blas.m4 b/m4/ax_blas.m4 new file mode 100644 index 0000000..9571905 --- /dev/null +++ b/m4/ax_blas.m4 @@ -0,0 +1,241 @@ +# =========================================================================== +# https://www.gnu.org/software/autoconf-archive/ax_blas.html +# =========================================================================== +# +# SYNOPSIS +# +# AX_BLAS([ACTION-IF-FOUND[, ACTION-IF-NOT-FOUND]]) +# +# DESCRIPTION +# +# This macro looks for a library that implements the BLAS linear-algebra +# interface (see http://www.netlib.org/blas/). On success, it sets the +# BLAS_LIBS output variable to hold the requisite library linkages. +# +# To link with BLAS, you should link with: +# +# $BLAS_LIBS $LIBS $FLIBS +# +# in that order. FLIBS is the output variable of the +# AC_F77_LIBRARY_LDFLAGS macro (called if necessary by AX_BLAS), and is +# sometimes necessary in order to link with F77 libraries. Users will also +# need to use AC_F77_DUMMY_MAIN (see the autoconf manual), for the same +# reason. +# +# Many libraries are searched for, from ATLAS to CXML to ESSL. The user +# may also use --with-blas= in order to use some specific BLAS +# library . In order to link successfully, however, be aware that you +# will probably need to use the same Fortran compiler (which can be set +# via the F77 env. var.) as was used to compile the BLAS library. +# +# ACTION-IF-FOUND is a list of shell commands to run if a BLAS library is +# found, and ACTION-IF-NOT-FOUND is a list of commands to run it if it is +# not found. If ACTION-IF-FOUND is not specified, the default action will +# define HAVE_BLAS. +# +# LICENSE +# +# Copyright (c) 2008 Steven G. Johnson +# Copyright (c) 2019 Geoffrey M. Oxberry +# +# This program is free software: you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program. If not, see . +# +# As a special exception, the respective Autoconf Macro's copyright owner +# gives unlimited permission to copy, distribute and modify the configure +# scripts that are the output of Autoconf when processing the Macro. You +# need not follow the terms of the GNU General Public License when using +# or distributing such scripts, even though portions of the text of the +# Macro appear in them. The GNU General Public License (GPL) does govern +# all other use of the material that constitutes the Autoconf Macro. +# +# This special exception to the GPL applies to versions of the Autoconf +# Macro released by the Autoconf Archive. When you make and distribute a +# modified version of the Autoconf Macro, you may extend this special +# exception to the GPL to apply to your modified version as well. + +#serial 17 + +AU_ALIAS([ACX_BLAS], [AX_BLAS]) +AC_DEFUN([AX_BLAS], [ +AC_PREREQ([2.55]) +AC_REQUIRE([AC_F77_LIBRARY_LDFLAGS]) +AC_REQUIRE([AC_CANONICAL_HOST]) +ax_blas_ok=no + +AC_ARG_WITH(blas, + [AS_HELP_STRING([--with-blas=], [use BLAS library ])]) +case $with_blas in + yes | "") ;; + no) ax_blas_ok=disable ;; + -* | */* | *.a | *.so | *.so.* | *.dylib | *.dylib.* | *.o) + BLAS_LIBS="$with_blas" + ;; + *) BLAS_LIBS="-l$with_blas" ;; +esac + +# Get fortran linker names of BLAS functions to check for. +AC_F77_FUNC(sgemm) +AC_F77_FUNC(dgemm) + +ax_blas_save_LIBS="$LIBS" +LIBS="$LIBS $FLIBS" + +# First, check BLAS_LIBS environment variable +if test $ax_blas_ok = no; then +if test "x$BLAS_LIBS" != x; then + save_LIBS="$LIBS"; LIBS="$BLAS_LIBS $LIBS" + AC_MSG_CHECKING([for $sgemm in $BLAS_LIBS]) + AC_LINK_IFELSE([AC_LANG_CALL([], [$sgemm])], [ax_blas_ok=yes], [BLAS_LIBS=""]) + AC_MSG_RESULT($ax_blas_ok) + LIBS="$save_LIBS" +fi +fi + +# BLAS linked to by default? (happens on some supercomputers) +if test $ax_blas_ok = no; then + save_LIBS="$LIBS"; LIBS="$LIBS" + AC_MSG_CHECKING([if $sgemm is being linked in already]) + AC_LINK_IFELSE([AC_LANG_CALL([], [$sgemm])], [ax_blas_ok=yes]) + AC_MSG_RESULT($ax_blas_ok) + LIBS="$save_LIBS" +fi + +# BLAS in OpenBLAS library? (http://xianyi.github.com/OpenBLAS/) +if test $ax_blas_ok = no; then + AC_CHECK_LIB(openblas, $sgemm, [ax_blas_ok=yes + BLAS_LIBS="-lopenblas"]) +fi + +# BLAS in ATLAS library? (http://math-atlas.sourceforge.net/) +if test $ax_blas_ok = no; then + AC_CHECK_LIB(atlas, ATL_xerbla, + [AC_CHECK_LIB(f77blas, $sgemm, + [AC_CHECK_LIB(cblas, cblas_dgemm, + [ax_blas_ok=yes + BLAS_LIBS="-lcblas -lf77blas -latlas"], + [], [-lf77blas -latlas])], + [], [-latlas])]) +fi + +# BLAS in PhiPACK libraries? (requires generic BLAS lib, too) +if test $ax_blas_ok = no; then + AC_CHECK_LIB(blas, $sgemm, + [AC_CHECK_LIB(dgemm, $dgemm, + [AC_CHECK_LIB(sgemm, $sgemm, + [ax_blas_ok=yes; BLAS_LIBS="-lsgemm -ldgemm -lblas"], + [], [-lblas])], + [], [-lblas])]) +fi + +# BLAS in Intel MKL library? +if test $ax_blas_ok = no; then + # MKL for gfortran + if test x"$ac_cv_fc_compiler_gnu" = xyes; then + # 64 bit + if test $host_cpu = x86_64; then + AC_CHECK_LIB(mkl_gf_lp64, $sgemm, + [ax_blas_ok=yes;BLAS_LIBS="-lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread"],, + [-lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread]) + # 32 bit + elif test $host_cpu = i686; then + AC_CHECK_LIB(mkl_gf, $sgemm, + [ax_blas_ok=yes;BLAS_LIBS="-lmkl_gf -lmkl_sequential -lmkl_core -lpthread"],, + [-lmkl_gf -lmkl_sequential -lmkl_core -lpthread]) + fi + # MKL for other compilers (Intel, PGI, ...?) + else + # 64-bit + if test $host_cpu = x86_64; then + AC_CHECK_LIB(mkl_intel_lp64, $sgemm, + [ax_blas_ok=yes;BLAS_LIBS="-lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread"],, + [-lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread]) + # 32-bit + elif test $host_cpu = i686; then + AC_CHECK_LIB(mkl_intel, $sgemm, + [ax_blas_ok=yes;BLAS_LIBS="-lmkl_intel -lmkl_sequential -lmkl_core -lpthread"],, + [-lmkl_intel -lmkl_sequential -lmkl_core -lpthread]) + fi + fi +fi +# Old versions of MKL +if test $ax_blas_ok = no; then + AC_CHECK_LIB(mkl, $sgemm, [ax_blas_ok=yes;BLAS_LIBS="-lmkl -lguide -lpthread"],,[-lguide -lpthread]) +fi + +# BLAS in Apple vecLib library? +if test $ax_blas_ok = no; then + save_LIBS="$LIBS"; LIBS="-framework vecLib $LIBS" + AC_MSG_CHECKING([for $sgemm in -framework vecLib]) + AC_LINK_IFELSE([AC_LANG_CALL([], [$sgemm])], [ax_blas_ok=yes;BLAS_LIBS="-framework vecLib"]) + AC_MSG_RESULT($ax_blas_ok) + LIBS="$save_LIBS" +fi + +# BLAS in Alpha CXML library? +if test $ax_blas_ok = no; then + AC_CHECK_LIB(cxml, $sgemm, [ax_blas_ok=yes;BLAS_LIBS="-lcxml"]) +fi + +# BLAS in Alpha DXML library? (now called CXML, see above) +if test $ax_blas_ok = no; then + AC_CHECK_LIB(dxml, $sgemm, [ax_blas_ok=yes;BLAS_LIBS="-ldxml"]) +fi + +# BLAS in Sun Performance library? +if test $ax_blas_ok = no; then + if test "x$GCC" != xyes; then # only works with Sun CC + AC_CHECK_LIB(sunmath, acosp, + [AC_CHECK_LIB(sunperf, $sgemm, + [BLAS_LIBS="-xlic_lib=sunperf -lsunmath" + ax_blas_ok=yes],[],[-lsunmath])]) + fi +fi + +# BLAS in SCSL library? (SGI/Cray Scientific Library) +if test $ax_blas_ok = no; then + AC_CHECK_LIB(scs, $sgemm, [ax_blas_ok=yes; BLAS_LIBS="-lscs"]) +fi + +# BLAS in SGIMATH library? +if test $ax_blas_ok = no; then + AC_CHECK_LIB(complib.sgimath, $sgemm, + [ax_blas_ok=yes; BLAS_LIBS="-lcomplib.sgimath"]) +fi + +# BLAS in IBM ESSL library? (requires generic BLAS lib, too) +if test $ax_blas_ok = no; then + AC_CHECK_LIB(blas, $sgemm, + [AC_CHECK_LIB(essl, $sgemm, + [ax_blas_ok=yes; BLAS_LIBS="-lessl -lblas"], + [], [-lblas $FLIBS])]) +fi + +# Generic BLAS library? +if test $ax_blas_ok = no; then + AC_CHECK_LIB(blas, $sgemm, [ax_blas_ok=yes; BLAS_LIBS="-lblas"]) +fi + +AC_SUBST(BLAS_LIBS) + +LIBS="$ax_blas_save_LIBS" + +# Finally, execute ACTION-IF-FOUND/ACTION-IF-NOT-FOUND: +if test x"$ax_blas_ok" = xyes; then + ifelse([$1],,AC_DEFINE(HAVE_BLAS,1,[Define if you have a BLAS library.]),[$1]) + : +else + ax_blas_ok=no + $2 +fi +])dnl AX_BLAS diff --git a/m4/ax_check_gnu_make.m4 b/m4/ax_check_gnu_make.m4 new file mode 100644 index 0000000..71ffe34 --- /dev/null +++ b/m4/ax_check_gnu_make.m4 @@ -0,0 +1,96 @@ +# =========================================================================== +# https://www.gnu.org/software/autoconf-archive/ax_check_gnu_make.html +# =========================================================================== +# +# SYNOPSIS +# +# AX_CHECK_GNU_MAKE([run-if-true],[run-if-false]) +# +# DESCRIPTION +# +# This macro searches for a GNU version of make. If a match is found: +# +# * The makefile variable `ifGNUmake' is set to the empty string, otherwise +# it is set to "#". This is useful for including a special features in a +# Makefile, which cannot be handled by other versions of make. +# * The makefile variable `ifnGNUmake' is set to #, otherwise +# it is set to the empty string. This is useful for including a special +# features in a Makefile, which can be handled +# by other versions of make or to specify else like clause. +# * The variable `_cv_gnu_make_command` is set to the command to invoke +# GNU make if it exists, the empty string otherwise. +# * The variable `ax_cv_gnu_make_command` is set to the command to invoke +# GNU make by copying `_cv_gnu_make_command`, otherwise it is unset. +# * If GNU Make is found, its version is extracted from the output of +# `make --version` as the last field of a record of space-separated +# columns and saved into the variable `ax_check_gnu_make_version`. +# * Additionally if GNU Make is found, run shell code run-if-true +# else run shell code run-if-false. +# +# Here is an example of its use: +# +# Makefile.in might contain: +# +# # A failsafe way of putting a dependency rule into a makefile +# $(DEPEND): +# $(CC) -MM $(srcdir)/*.c > $(DEPEND) +# +# @ifGNUmake@ ifeq ($(DEPEND),$(wildcard $(DEPEND))) +# @ifGNUmake@ include $(DEPEND) +# @ifGNUmake@ else +# fallback code +# @ifGNUmake@ endif +# +# Then configure.in would normally contain: +# +# AX_CHECK_GNU_MAKE() +# AC_OUTPUT(Makefile) +# +# Then perhaps to cause gnu make to override any other make, we could do +# something like this (note that GNU make always looks for GNUmakefile +# first): +# +# if ! test x$_cv_gnu_make_command = x ; then +# mv Makefile GNUmakefile +# echo .DEFAULT: > Makefile ; +# echo \ $_cv_gnu_make_command \$@ >> Makefile; +# fi +# +# Then, if any (well almost any) other make is called, and GNU make also +# exists, then the other make wraps the GNU make. +# +# LICENSE +# +# Copyright (c) 2008 John Darrington +# Copyright (c) 2015 Enrico M. Crisostomo +# +# Copying and distribution of this file, with or without modification, are +# permitted in any medium without royalty provided the copyright notice +# and this notice are preserved. This file is offered as-is, without any +# warranty. + +#serial 12 + +AC_DEFUN([AX_CHECK_GNU_MAKE],dnl + [AC_PROG_AWK + AC_CACHE_CHECK([for GNU make],[_cv_gnu_make_command],[dnl + _cv_gnu_make_command="" ; +dnl Search all the common names for GNU make + for a in "$MAKE" make gmake gnumake ; do + if test -z "$a" ; then continue ; fi ; + if "$a" --version 2> /dev/null | grep GNU 2>&1 > /dev/null ; then + _cv_gnu_make_command=$a ; + AX_CHECK_GNU_MAKE_HEADLINE=$("$a" --version 2> /dev/null | grep "GNU Make") + ax_check_gnu_make_version=$(echo ${AX_CHECK_GNU_MAKE_HEADLINE} | ${AWK} -F " " '{ print $(NF); }') + break ; + fi + done ;]) +dnl If there was a GNU version, then set @ifGNUmake@ to the empty string, '#' otherwise + AS_VAR_IF([_cv_gnu_make_command], [""], [AS_VAR_SET([ifGNUmake], ["#"])], [AS_VAR_SET([ifGNUmake], [""])]) + AS_VAR_IF([_cv_gnu_make_command], [""], [AS_VAR_SET([ifnGNUmake], [""])], [AS_VAR_SET([ifnGNUmake], ["#"])]) + AS_VAR_IF([_cv_gnu_make_command], [""], [AS_UNSET(ax_cv_gnu_make_command)], [AS_VAR_SET([ax_cv_gnu_make_command], [${_cv_gnu_make_command}])]) + AS_VAR_IF([_cv_gnu_make_command], [""],[$2],[$1]) + AC_SUBST([ifGNUmake]) + AC_SUBST([ifnGNUmake]) +]) + diff --git a/m4/ax_irpf90.m4 b/m4/ax_irpf90.m4 new file mode 100644 index 0000000..50942a2 --- /dev/null +++ b/m4/ax_irpf90.m4 @@ -0,0 +1,73 @@ +# =========================================================================== +# https://www.gnu.org/software/autoconf-archive/ax_irpf90.html +# =========================================================================== +# +# SYNOPSIS +# +# AX_PROG_IRPF90 +# +# DESCRIPTION +# +# Check for the program 'irpf90', let script continue if exists, pops up +# error message if not. +# +# Besides checking existence, this macro also set these environment +# variables upon completion: +# +# IRPF90 = which irpf90 +# +# DEPENDENCIES +# +# AX_CHECK_GNU_MAKE +# +# LICENSE +# +# Copyright (c) 2021 Anthony Scemama +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 2 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program. If not, see . +# +# As a special exception, the respective Autoconf Macro's copyright owner +# gives unlimited permission to copy, distribute and modify the configure +# scripts that are the output of Autoconf when processing the Macro. You +# need not follow the terms of the GNU General Public License when using +# or distributing such scripts, even though portions of the text of the +# Macro appear in them. The GNU General Public License (GPL) does govern +# all other use of the material that constitutes the Autoconf Macro. +# +# This special exception to the GPL applies to versions of the Autoconf +# Macro released by the Autoconf Archive. When you make and distribute a +# modified version of the Autoconf Macro, you may extend this special +# exception to the GPL to apply to your modified version as well. + + +AU_ALIAS([AC_PROG_IRPF90], [AX_PROG_IRPF90]) + +AC_DEFUN([AX_PROG_IRPF90], [ + +# Requirements +AC_REQUIRE([AX_CHECK_GNU_MAKE]) + +AS_IF([test "x$ifGNUmake" = "x#"], [ AC_MSG_ERROR([GNU Make (gmake) is required with IRPF90]) ]) + +# IRPF90 +AC_PATH_PROG([IRPF90], [irpf90], [nocommand]) +AS_IF([test "x$IRPF90" = xnocommand], [ + AC_MSG_ERROR([irpf90 not found in $PATH]) ]) + +AC_FC_FREEFORM +AC_FC_LINE_LENGTH +AC_FC_MODULE_EXTENSION + +]) + diff --git a/m4/ax_lapack.m4 b/m4/ax_lapack.m4 new file mode 100644 index 0000000..86265ff --- /dev/null +++ b/m4/ax_lapack.m4 @@ -0,0 +1,135 @@ +# =========================================================================== +# https://www.gnu.org/software/autoconf-archive/ax_lapack.html +# =========================================================================== +# +# SYNOPSIS +# +# AX_LAPACK([ACTION-IF-FOUND[, ACTION-IF-NOT-FOUND]]) +# +# DESCRIPTION +# +# This macro looks for a library that implements the LAPACK linear-algebra +# interface (see http://www.netlib.org/lapack/). On success, it sets the +# LAPACK_LIBS output variable to hold the requisite library linkages. +# +# To link with LAPACK, you should link with: +# +# $LAPACK_LIBS $BLAS_LIBS $LIBS $FLIBS +# +# in that order. BLAS_LIBS is the output variable of the AX_BLAS macro, +# called automatically. FLIBS is the output variable of the +# AC_F77_LIBRARY_LDFLAGS macro (called if necessary by AX_BLAS), and is +# sometimes necessary in order to link with F77 libraries. Users will also +# need to use AC_F77_DUMMY_MAIN (see the autoconf manual), for the same +# reason. +# +# The user may also use --with-lapack= in order to use some specific +# LAPACK library . In order to link successfully, however, be aware +# that you will probably need to use the same Fortran compiler (which can +# be set via the F77 env. var.) as was used to compile the LAPACK and BLAS +# libraries. +# +# ACTION-IF-FOUND is a list of shell commands to run if a LAPACK library +# is found, and ACTION-IF-NOT-FOUND is a list of commands to run it if it +# is not found. If ACTION-IF-FOUND is not specified, the default action +# will define HAVE_LAPACK. +# +# LICENSE +# +# Copyright (c) 2009 Steven G. Johnson +# Copyright (c) 2019 Geoffrey M. Oxberry +# +# This program is free software: you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program. If not, see . +# +# As a special exception, the respective Autoconf Macro's copyright owner +# gives unlimited permission to copy, distribute and modify the configure +# scripts that are the output of Autoconf when processing the Macro. You +# need not follow the terms of the GNU General Public License when using +# or distributing such scripts, even though portions of the text of the +# Macro appear in them. The GNU General Public License (GPL) does govern +# all other use of the material that constitutes the Autoconf Macro. +# +# This special exception to the GPL applies to versions of the Autoconf +# Macro released by the Autoconf Archive. When you make and distribute a +# modified version of the Autoconf Macro, you may extend this special +# exception to the GPL to apply to your modified version as well. + +#serial 10 + +AU_ALIAS([ACX_LAPACK], [AX_LAPACK]) +AC_DEFUN([AX_LAPACK], [ +AC_REQUIRE([AX_BLAS]) +ax_lapack_ok=no + +AC_ARG_WITH(lapack, + [AS_HELP_STRING([--with-lapack=], [use LAPACK library ])]) +case $with_lapack in + yes | "") ;; + no) ax_lapack_ok=disable ;; + -* | */* | *.a | *.so | *.so.* | *.dylib | *.dylib.* | *.o) + LAPACK_LIBS="$with_lapack" + ;; + *) LAPACK_LIBS="-l$with_lapack" ;; +esac + +# Get fortran linker name of LAPACK function to check for. +AC_F77_FUNC(cheev) + +# We cannot use LAPACK if BLAS is not found +if test "x$ax_blas_ok" != xyes; then + ax_lapack_ok=noblas + LAPACK_LIBS="" +fi + +# First, check LAPACK_LIBS environment variable +if test "x$LAPACK_LIBS" != x; then + save_LIBS="$LIBS"; LIBS="$LAPACK_LIBS $BLAS_LIBS $LIBS $FLIBS" + AC_MSG_CHECKING([for $cheev in $LAPACK_LIBS]) + AC_LINK_IFELSE([AC_LANG_CALL([], [$cheev])], [ax_lapack_ok=yes], [LAPACK_LIBS=""]) + AC_MSG_RESULT($ax_lapack_ok) + LIBS="$save_LIBS" + if test $ax_lapack_ok = no; then + LAPACK_LIBS="" + fi +fi + +# LAPACK linked to by default? (is sometimes included in BLAS lib) +if test $ax_lapack_ok = no; then + save_LIBS="$LIBS"; LIBS="$LIBS $BLAS_LIBS $FLIBS" + AC_CHECK_FUNC($cheev, [ax_lapack_ok=yes]) + LIBS="$save_LIBS" +fi + +# Generic LAPACK library? +for lapack in lapack lapack_rs6k; do + if test $ax_lapack_ok = no; then + save_LIBS="$LIBS"; LIBS="$BLAS_LIBS $LIBS" + AC_CHECK_LIB($lapack, $cheev, + [ax_lapack_ok=yes; LAPACK_LIBS="-l$lapack"], [], [$FLIBS]) + LIBS="$save_LIBS" + fi +done + +AC_SUBST(LAPACK_LIBS) + +# Finally, execute ACTION-IF-FOUND/ACTION-IF-NOT-FOUND: +if test x"$ax_lapack_ok" = xyes; then + ifelse([$1],,AC_DEFINE(HAVE_LAPACK,1,[Define if you have LAPACK library.]),[$1]) + : +else + ax_lapack_ok=no + $2 +fi +])dnl AX_LAPACK + diff --git a/m4/ax_zmq.m4 b/m4/ax_zmq.m4 new file mode 100644 index 0000000..83e7b78 --- /dev/null +++ b/m4/ax_zmq.m4 @@ -0,0 +1,91 @@ +# =========================================================================== +# https://www.gnu.org/software/autoconf-archive/ax_zmq.html +# =========================================================================== +# +# SYNOPSIS +# +# AX_ZMQ([MINIMUM-VERSION], [ACTION-IF-FOUND], [ACTION-IF-NOT-FOUND]) +# +# DESCRIPTION +# +# Test for the ZMQ libraries of a particular version (or newer). The +# default version tested for is 4.0.0. +# +# The macro tests for ZMQ libraries in the library/include path, and, when +# provided, also in the path given by --with-zmq. +# +# This macro calls: +# +# AC_SUBST(ZMQ_CPPFLAGS) / AC_SUBST(ZMQ_LDFLAGS) / AC_SUBST(ZMQ_LIBS) +# +# And sets: +# +# HAVE_ZMQ +# +# LICENSE +# +# Copyright (c) 2016 Jeroen Meijer +# +# Copying and distribution of this file, with or without modification, are +# permitted in any medium without royalty provided the copyright notice +# and this notice are preserved. This file is offered as-is, without any +# warranty. + +#serial 3 + +AC_DEFUN([AX_ZMQ], [ + AC_ARG_WITH([zmq], [AS_HELP_STRING([--with-zmq=],[ZMQ prefix directory])], [ + ZMQ_LDFLAGS="-L${with_zmq}/lib" + ZMQ_CPPFLAGS="-I${with_zmq}/include" + ]) + + HAVE_ZMQ=0 + if test "$with_zmq" != "no"; then + + LD_FLAGS="$LDFLAGS $ZMQ_LDFLAGS" + CPPFLAGS="$CPPFLAGS $ZMQ_CPPFLAGS" + + AC_LANG_PUSH([C]) + AC_CHECK_HEADER(zmq.h, [zmq_h=yes], [zmq_h=no]) + AC_LANG_POP([C]) + + if test "$zmq_h" = "yes"; then + version=ifelse([$1], ,4.0.0,$1) + AC_MSG_CHECKING([for ZMQ version >= $version]) + version=$(echo $version | tr '.' ',') + AC_EGREP_CPP([version_ok], [ +#include +#if defined(ZMQ_VERSION) && ZMQ_VERSION >= ZMQ_MAKE_VERSION($version) + version_ok +#endif + ],[ + AC_MSG_RESULT(yes) + AC_CHECK_LIB([zmq],[zmq_send],[ + HAVE_ZMQ=1 + ZMQ_LIBS="-lzmq" + ],[ + AC_MSG_ERROR([libzmq.so not in LIBRARY_PATH]) + ]) + + AC_SUBST(ZMQ_LDFLAGS) + AC_SUBST(ZMQ_CPPFLAGS) + AC_SUBST(ZMQ_LIBS) + ], AC_MSG_RESULT([no valid ZMQ version was found])) + else + AC_MSG_WARN([no valid ZMQ installation was found]) + fi + + if test $HAVE_ZMQ = 1; then + # execute ACTION-IF-FOUND (if present): + ifelse([$2], , :, [$2]) + else + # execute ACTION-IF-NOT-FOUND (if present): + ifelse([$3], , :, [$3]) + fi + else + AC_MSG_NOTICE([not checking for ZMQ]) + fi + + AC_DEFINE(HAVE_ZMQ,,[define if the ZMQ library is available]) +]) + diff --git a/ocaml/Block.ml b/ocaml/Block.ml index 21cbf62..6790df7 100644 --- a/ocaml/Block.ml +++ b/ocaml/Block.ml @@ -1,6 +1,6 @@ open Qptypes -type t = +type t = { property : Property.t ; value : Sample.t ; weight : Weight.t ; @@ -12,7 +12,7 @@ type t = let re = Str.regexp "[ |#|\n]+" -let of_string s = +let of_string s = try let lst = @@ -23,24 +23,24 @@ let of_string s = | b :: pid :: c:: p :: w :: v :: [] -> Some { property = Property.of_string p ; value = Sample.of_float (float_of_string v) ; - weight = Weight.of_float (float_of_string w) ; + weight = Weight.of_float (float_of_string w) ; compute_node = Compute_node.of_string c; pid = int_of_string pid; block_id = Block_id.of_int (int_of_string b) ; } - | b :: pid :: c:: p :: w :: v -> - let v = + | b :: pid :: c:: p :: w :: v -> + let v = List.rev v - |> Array.of_list + |> Array.of_list |> Array.map float_of_string in - let dim = + let dim = Array.length v in Some { property = Property.of_string p ; value = Sample.of_float_array ~dim v ; - weight = Weight.of_float (float_of_string w) ; + weight = Weight.of_float (float_of_string w) ; compute_node = Compute_node.of_string c; pid = int_of_string pid; block_id = Block_id.of_int (int_of_string b) ; @@ -50,20 +50,105 @@ let of_string s = | _ -> None - + +let to_short_string b = + Printf.sprintf "%s # %s %d %d" + (Property.to_string b.property) + (Compute_node.to_string b.compute_node) + b.pid + (Block_id.to_int b.block_id) + + let to_string b = - Printf.sprintf "%s %s # %s %s %s %d" - (Sample.to_string b.value ) - (Weight.to_float b.weight |> string_of_float) - (Property.to_string b.property) - (Compute_node.to_string b.compute_node) - (string_of_int b.pid) - (Block_id.to_int b.block_id) + Printf.sprintf "%s %s # %s %s %s %d" + (Sample.to_string b.value ) + (Weight.to_float b.weight |> string_of_float) + (Property.to_string b.property) + (Compute_node.to_string b.compute_node) + (string_of_int b.pid) + (Block_id.to_int b.block_id) + +let zero = + bytes_of_int 0 + +let to_bytes b = + (* [ Length of b + [ Length of value ; + Value ; + Length of weight ; + Weight ; + ... ] ] *) + let l = + [ Property.to_bytes b.property ; + Sample.to_bytes b.value ; + Weight.to_bytes b.weight ; + bytes_of_int b.pid ; + Block_id.to_bytes b.block_id ; + Compute_node.to_bytes b.compute_node ] + |> List.map (fun x -> [ bytes_of_int (Bytes.length x) ; x ] ) + |> List.concat + in + let result = + Bytes.concat Bytes.empty (zero :: l) + in + Bytes.set_int64_le result 0 (Int64.of_int ((Bytes.length result) - 8)); + result + + +let read_bytes b = + (* Reads m, the first 8 bytes as an int64 containing the number of bytes to read. + Then, read the next m bytes and return a tuple containing the decoded data and the rest. + *) + let l = Bytes.length b in + if l < 8 then + failwith "Zero-sized bytes" + else + let m = + Bytes.get_int64_le b 0 + |> Int64.to_int + in + let nl = l-m-8 in + if nl > 0 then + (Bytes.sub b 8 m, Some (Bytes.sub b (8+m) nl)) + else + (Bytes.sub b 8 m, None) + + + +let of_bytes b = + let rec loop accu s = + match read_bytes s with + | data, None -> List.rev (data :: accu) + | data, (Some rest) -> loop (data :: accu) rest + in + let result = + match loop [] b with + | property :: value :: weight :: pid :: block_id :: compute_node :: [] -> + Some + { property = Property.of_bytes property; + value = Sample.of_bytes value; + weight = Weight.of_bytes weight; + pid = int_of_bytes pid; + block_id = Block_id.of_bytes block_id; + compute_node = Compute_node.of_bytes compute_node; + } + | _ -> None + in + result + + +let of_string_or_bytes s = + if Qmcchem_config.binary_io then + Bytes.of_string s + |> of_bytes + else + of_string s + let dir_name = lazy( - let ezfio_filename = + let ezfio_filename = Lazy.force Qputils.ezfio_filename in let md5 = @@ -72,8 +157,8 @@ let dir_name = lazy( let d = Filename.concat ezfio_filename "blocks" in if not ( Sys.file_exists d ) then Unix.mkdir d 0o755; - List.fold_right Filename.concat - [ ezfio_filename ; "blocks" ; md5 ; Filename.dir_sep ] "" + List.fold_right Filename.concat + [ ezfio_filename ; "blocks" ; md5 ; Filename.dir_sep ] "" ) @@ -84,11 +169,11 @@ let _raw_data = let update_raw_data ?(locked=true) () = (* Create array of files to read *) - let dir_name = + let dir_name = Lazy.force dir_name in - let files = - let result = + let files = + let result = if Sys.file_exists dir_name && Sys.is_directory dir_name then begin Sys.readdir dir_name @@ -102,7 +187,7 @@ let update_raw_data ?(locked=true) () = else List.filter (fun x -> try - let _ = + let _ = Str.search_backward (Str.regexp "locked") x ((String.length x) - 1) in false with @@ -110,7 +195,7 @@ let update_raw_data ?(locked=true) () = ) result in - let rec transform new_list = function + let rec transform new_list = function | [] -> new_list | head :: tail -> let head = String.trim head in @@ -120,33 +205,74 @@ let update_raw_data ?(locked=true) () = | Some x -> transform (x::new_list) tail in - let result = - let rec aux ic accu = - let l = - try - Some (input_line ic) - with - | End_of_file -> None + if Qmcchem_config.binary_io then + begin + let result = + let rec aux buf accu = + (* Read one block *) + let item, rest = + read_bytes buf + in + match of_bytes item with + | None -> [] + | Some item -> + match rest with + | None -> List.rev (item::accu) + | Some rest -> (aux [@tailcall]) rest (item::accu) + in + List.concat_map (fun filename -> + let ic = open_in filename in + let length = in_channel_length ic in + let result = + if length > 0 then + let buf = Bytes.create length in + really_input ic buf 0 length; + aux buf [] + else [] + in + close_in ic; + result ) files in - match l with - | None -> List.rev accu - | Some l -> (aux [@tailcall]) ic (l::accu) - in - List.concat_map (fun filename -> - let ic = open_in filename in - let result = aux ic [] in - close_in ic; - result ) files - |> transform [] - in - result + result + end + else + begin + let result = + let rec aux ic accu = + let l = + try + Some (input_line ic) + with + | End_of_file -> None + in + match l with + | None -> List.rev accu + | Some l -> (aux [@tailcall]) ic (l::accu) + in + List.concat_map (fun filename -> + let ic = open_in filename in + let result = aux ic [] in + close_in ic; + result ) files + |> transform [] + in + result + end -let raw_data ?(locked=true) () = +let to_string_or_bytes b = + if Qmcchem_config.binary_io then + to_bytes b + |> Bytes.to_string + else + to_string b + + +let raw_data ?(locked=true) () = match !_raw_data with | Some x -> x | None -> - let result = + let result = update_raw_data ~locked () in _raw_data := Some result; @@ -156,7 +282,7 @@ let raw_data ?(locked=true) () = let properties = lazy ( let h = Hashtbl.create 63 in - List.iter (fun x -> + List.iter (fun x -> Hashtbl.replace h (Property.to_string x.property) x.property) (raw_data ()); Hashtbl.fold (fun k v a -> v :: a) h [] diff --git a/ocaml/Message.ml b/ocaml/Message.ml index 1c34cea..8f6bb11 100644 --- a/ocaml/Message.ml +++ b/ocaml/Message.ml @@ -11,7 +11,7 @@ type t = | Error of string -let create m = +let create m = try match m with | [ "cpu" ; c ; pid ; b ; "1" ; v ] -> @@ -23,7 +23,7 @@ let create m = compute_node = Compute_node.of_string c; pid = int_of_string pid; block_id = Block_id.of_int (int_of_string b); - } + } | [ "accep" ; c ; pid ; b ; "1" ; v ] -> let open Block in Property @@ -33,8 +33,8 @@ let create m = compute_node = Compute_node.of_string c; pid = int_of_string pid; block_id = Block_id.of_int (int_of_string b); - } - | [ prop ; c ; pid ; b ; w ; v ] -> + } + | [ prop ; c ; pid ; b ; w ; v ] -> let open Block in Property { property = Property.of_string prop; @@ -43,19 +43,19 @@ let create m = compute_node = Compute_node.of_string c; pid = int_of_string pid; block_id = Block_id.of_int (int_of_string b); - } + } | "elec_coord" :: c :: pid :: _ :: n ::walkers -> begin let elec_num = Lazy.force Qputils.elec_num - and n = + and n = int_of_string n in assert (n = List.length walkers); let rec build_walker accu = function - | (0,tail) -> - let result = - List.rev accu + | (0,tail) -> + let result = + List.rev accu |> List.rev_map float_of_string |> List.rev |> Array.of_list @@ -65,10 +65,10 @@ let create m = build_walker (head::accu) (n-1, tail) | _ -> failwith "Bad walkers" in - let rec build accu = function + let rec build accu = function | [] -> Array.of_list accu - | w -> - let (result, tail) = + | w -> + let (result, tail) = build_walker [] (3*elec_num+3, w) in build (result::accu) tail @@ -80,13 +80,13 @@ let create m = | [ "unregister" ; c ; pid ] -> Unregister (Compute_node.of_string c, int_of_string pid) | [ "Test" ] -> Test | [ "Ezfio" ; ezfio_msg ] -> Ezfio ezfio_msg - | prop :: c :: pid :: b :: d :: w :: l -> + | prop :: c :: pid :: b :: d :: w :: l -> let property = Property.of_string prop in begin assert (not (Property.is_scalar property)); - let a = + let a = Array.of_list l |> Array.map float_of_string and dim = @@ -101,7 +101,7 @@ let create m = compute_node = Compute_node.of_string c; pid = int_of_string pid; block_id = Block_id.of_int (int_of_string b); - } + } end | l -> Error (String.concat ":" l) with @@ -114,9 +114,9 @@ let to_string = function | Walkers (h,p,w) -> Printf.sprintf "Walkers : %s %d : %d walkers" (Compute_node.to_string h) p (Array.length w) | GetWalkers n -> Printf.sprintf "GetWalkers %d" (Strictly_positive_int.to_int n) - | Register (h,p) -> Printf.sprintf "Register : %s %d" + | Register (h,p) -> Printf.sprintf "Register : %s %d" (Compute_node.to_string h) p - | Unregister (h,p) -> Printf.sprintf "Unregister : %s %d" + | Unregister (h,p) -> Printf.sprintf "Unregister : %s %d" (Compute_node.to_string h) p | Test -> "Test" | Ezfio msg -> "Ezfio "^msg diff --git a/ocaml/Qmcchem_config.ml b/ocaml/Qmcchem_config.ml index ac88031..dd3679e 100644 --- a/ocaml/Qmcchem_config.ml +++ b/ocaml/Qmcchem_config.ml @@ -8,10 +8,10 @@ let root = lazy ( (* PATH environment variable as a list of strings *) let path = lazy ( - let p = + let p = try Sys.getenv "PATH" with | Not_found -> failwith "PATH environment variable is not set" - in + in String.split_on_char ':' p ) @@ -23,8 +23,8 @@ let full_path exe = | [] -> None | head :: tail -> begin - let fp = - Filename.concat head exe + let fp = + Filename.concat head exe in if Sys.file_exists fp then Some fp @@ -82,8 +82,8 @@ let hostname = lazy ( ) -external get_ipv4_address_for_interface : string -> string = - "get_ipv4_address_for_interface" +external get_ipv4_address_for_interface : string -> string = + "get_ipv4_address_for_interface" let ip_address = lazy ( @@ -106,7 +106,7 @@ let ip_address = lazy ( end | Some interface -> let result = get_ipv4_address_for_interface interface in - if String.sub result 0 5 = "error" then + if String.sub result 0 5 = "error" then Printf.sprintf "Unable to use network interface %s" interface |> failwith else @@ -115,3 +115,8 @@ let ip_address = lazy ( +let binary_io = + try + Sys.getenv "QMCCHEM_IO" = "B" + with Not_found -> false + diff --git a/ocaml/Qmcchem_dataserver.ml b/ocaml/Qmcchem_dataserver.ml index c009bd2..c46c3e5 100644 --- a/ocaml/Qmcchem_dataserver.ml +++ b/ocaml/Qmcchem_dataserver.ml @@ -17,9 +17,9 @@ let initialization_timeout = 600. let bind_socket ~socket_type ~socket ~address = try Zmq.Socket.bind socket address - with - | Unix.Unix_error (_, message, f) -> - failwith @@ Printf.sprintf + with + | Unix.Unix_error (_, message, f) -> + failwith @@ Printf.sprintf "\n%s\nUnable to bind the dataserver's %s socket :\n %s\n%s" f socket_type address message | other_exception -> raise other_exception @@ -30,7 +30,7 @@ let run ?(daemon=true) ezfio_filename = Qputils.set_ezfio_filename ezfio_filename; (** Check if walkers need to be created. *) - let () = + let () = if ( not(Ezfio.has_electrons_elec_coord_pool ()) ) then begin Printf.printf "Generating initial walkers...\n%!"; @@ -39,12 +39,12 @@ let run ?(daemon=true) ezfio_filename = Unix.execvp (Lazy.force Qmcchem_config.qmc_create_walkers) [|"qmc_create_walkers" ; ezfio_filename|] - | pid -> + | pid -> begin ignore @@ Unix.waitpid [] pid; Printf.printf "Initial walkers ready\n%!" end - end + end in @@ -74,7 +74,7 @@ let run ?(daemon=true) ezfio_filename = in (** Status variable (mutable) *) - let status = + let status = ref Status.Queued in @@ -91,7 +91,7 @@ let run ?(daemon=true) ezfio_filename = let adress_prefix = "tcp://*:" in - let result = + let result = List.fold_left (fun accu i -> let address = adress_prefix ^ (string_of_int (n+i)) @@ -99,14 +99,14 @@ let run ?(daemon=true) ezfio_filename = let socket = Zmq.Socket.create zmq_context Zmq.Socket.rep in - let result = + let result = try Zmq.Socket.bind socket address; accu with | _ -> false in - Zmq.Socket.close socket; + Zmq.Socket.close socket; result ) true [0;1;2;3] in @@ -120,7 +120,7 @@ let run ?(daemon=true) ezfio_filename = (** Random port number between 49152 and 65535 *) let port = - let newport = + let newport = ref ( 1024 + (Random.int (49151-1024))) in while ((check_port !newport) = `Unavailable) @@ -131,7 +131,7 @@ let run ?(daemon=true) ezfio_filename = in - let debug_socket = + let debug_socket = Zmq.Socket.create zmq_context Zmq.Socket.xpub and address = Printf.sprintf "tcp://*:%d" (port+4) @@ -140,17 +140,17 @@ let run ?(daemon=true) ezfio_filename = Zmq.Socket.set_linger_period debug_socket 100 ; let close_debug_socket () = - Zmq.Socket.close debug_socket + Zmq.Socket.close debug_socket in (** Sends a log text to the debug socket *) let send_log socket size t0 text = - let dt = + let dt = delta_t t0 in - let message = - Printf.sprintf "%20s : %8d : %10s : %e" - socket size text dt + let message = + Printf.sprintf "%20s : %8d : %10s : %e" + socket size text dt in Zmq.Socket.send debug_socket message in @@ -159,7 +159,7 @@ let run ?(daemon=true) ezfio_filename = (** {2 Walkers} *) (** Number of electrons *) - let elec_num = + let elec_num = Lazy.force Qputils.elec_num in @@ -181,14 +181,14 @@ let run ?(daemon=true) ezfio_filename = in (** Array of walkers. The size is [walk_num_tot]. *) - let walkers_array = - let t0 = + let walkers_array = + let t0 = Unix.gettimeofday () in let j = 3*elec_num + 3 in - let result = + let result = let size = Ezfio.get_electrons_elec_coord_pool_size () and ez = @@ -204,7 +204,7 @@ let run ?(daemon=true) ezfio_filename = failwith "Walkers file is broken." in String.concat " " [ "Read" ; string_of_int (Array.length result) ; - "walkers"] + "walkers"] |> send_log "status" 0 t0 ; result in @@ -217,7 +217,7 @@ let run ?(daemon=true) ezfio_filename = (** Last time when the walkers were saved to disk. *) - let last_save_walkers = + let last_save_walkers = ref 0. in @@ -227,14 +227,12 @@ let run ?(daemon=true) ezfio_filename = if (delta_t !last_save_walkers > 10. ) then begin Ezfio.set_electrons_elec_coord_pool_size walk_num_tot ; - let walkers_list = - Array.map Array.to_list walkers_array - |> Array.to_list - |> List.concat - |> List.rev_map float_of_string - |> List.rev + let walkers_list = + Array.to_list walkers_array + |> Array.concat + |> Array.map float_of_string in - Ezfio.set_electrons_elec_coord_pool (Ezfio.ezfio_array_of_list + Ezfio.set_electrons_elec_coord_pool (Ezfio.ezfio_array_of_array ~rank:3 ~dim:[| elec_num+1 ; 3 ; walk_num_tot |] ~data:walkers_list); let t0 = Unix.gettimeofday () @@ -263,7 +261,7 @@ let run ?(daemon=true) ezfio_filename = (** {2 Set of workers} *) - (** A hash table is kept to track the running workers. The keys are the + (** A hash table is kept to track the running workers. The keys are the built as string containing the couple ([Compute_node], [PID]), and the values are the last communication time. *) @@ -271,16 +269,16 @@ let run ?(daemon=true) ezfio_filename = (** The hash table for workers *) - let workers_hash = + let workers_hash = Hashtbl.create 63 in (** Creates a key using the couple ([Compute_node], [PID]) *) - let key compute_node pid = - String.concat " " [ - (Compute_node.to_string compute_node); + let key compute_node pid = + String.concat " " [ + (Compute_node.to_string compute_node); (string_of_int pid) ] in @@ -294,7 +292,7 @@ let run ?(daemon=true) ezfio_filename = in match Hashtbl.find_opt workers_hash s with | Some _ -> failwith (s^" already registered") - | None -> Hashtbl.add workers_hash s (Unix.gettimeofday ()) + | None -> Hashtbl.add workers_hash s (Unix.gettimeofday ()) in @@ -307,14 +305,14 @@ let run ?(daemon=true) ezfio_filename = in match Hashtbl.find_opt workers_hash s with | Some x -> Hashtbl.remove workers_hash s - | None -> failwith (s^" not registered") + | None -> failwith (s^" not registered") in (** Sets the last access of the worker to [Unix.gettimeofday ()] *) - let touch_worker w pid = - let s = + let touch_worker w pid = + let s = key w pid in Hashtbl.replace workers_hash s (Unix.gettimeofday ()) @@ -326,11 +324,11 @@ let run ?(daemon=true) ezfio_filename = let delta = initialization_timeout +. block_time *. 2. in - Hashtbl.fold (fun k v accu -> + Hashtbl.fold (fun k v accu -> if (now -. v) <= delta then v :: accu else - accu ) hash [] + accu ) hash [] |> List.length in @@ -344,14 +342,14 @@ let run ?(daemon=true) ezfio_filename = (** Name of the blocks file written by the current process. *) - let block_channel_filename = - let dirname = + let block_channel_filename = + let dirname = Lazy.force Block.dir_name in - let () = + let () = if not ( Sys.file_exists dirname ) then Unix.mkdir dirname 0o755 - in + in Filename.concat dirname ( hostname ^ "." ^ (string_of_int dataserver_pid) ) @@ -387,7 +385,7 @@ let run ?(daemon=true) ezfio_filename = the compute nodes. Happens when [max_file_size] is reached. *) let compress_block_file filename = - let t0 = + let t0 = Unix.gettimeofday () in close_out !block_channel; @@ -406,16 +404,16 @@ let run ?(daemon=true) ezfio_filename = (** {3 Status thread} *) - let start_status_thread = + let start_status_thread = let t0 = Unix.gettimeofday () in - Thread.create (fun () -> + Thread.create (fun () -> send_log "status" 0 t0 "Starting status thread"; let socket = Zmq.Socket.create zmq_context Zmq.Socket.pub - and address = + and address = Printf.sprintf "tcp://*:%d" (port+1) in bind_socket "PUB" socket address; @@ -423,15 +421,15 @@ let run ?(daemon=true) ezfio_filename = and delay_read = 2. in - let start_time = + let start_time = Unix.gettimeofday () - and stop_time = + and stop_time = ref (Input.Stop_time.(read () |> to_float) ) in - + let last_update = ref start_time - in + in while (!status <> Status.Stopped) do @@ -439,46 +437,46 @@ let run ?(daemon=true) ezfio_filename = let now = Unix.gettimeofday () in - let status_string = + let status_string = Status.to_string !status in Zmq.Socket.send socket status_string; send_log "status" (String.length status_string) now status_string; - let test = + let test = if (now -. !last_update > delay_read) then - let n_connect = - n_connected workers_hash now + let n_connect = + n_connected workers_hash now in `Update n_connect else if (now -. start_time > !stop_time) then - `Terminate + `Terminate else if (now -. start_time > initialization_timeout) then - `Timeout + `Timeout else - `None + `None in match (daemon, !status, test) with | (_ , _ , `None ) -> () | (_ , Status.Running , `Terminate ) -> change_status Status.Stopping | (false, Status.Running , `Update 0 ) -> change_status Status.Stopped - | (true , Status.Running , `Update 0 ) -> change_status Status.Queued - | (_ , _ , `Update i ) -> + | (true , Status.Running , `Update 0 ) -> change_status Status.Queued + | (_ , _ , `Update i ) -> begin status := Status.read (); last_update := now; stop_time := Input.Stop_time.(read () |> to_float) ; let n_tot = - Hashtbl.length workers_hash + Hashtbl.length workers_hash in - if (i <> n_tot) then + if (i <> n_tot) then begin Printf.sprintf "Connected workers : %d / %d" i n_tot - |> send_log "status" 0 now + |> send_log "status" 0 now end end - | (false, Status.Queued , `Timeout ) -> change_status Status.Stopped + | (false, Status.Queued , `Timeout ) -> change_status Status.Stopped | (_, _, _) -> () ; done; @@ -487,37 +485,37 @@ let run ?(daemon=true) ezfio_filename = Zmq.Socket.close socket ) in - + (** {3 Log thread} *) - let start_log_thread = + let start_log_thread = let t0 = Unix.gettimeofday () in - Thread.create (fun () -> + Thread.create (fun () -> send_log "status" 0 t0 "Starting log thread"; let socket = Zmq.Socket.create zmq_context Zmq.Socket.xsub - and address = + and address = Printf.sprintf "tcp://*:%d" (port+3) in bind_socket "XSUB" socket address; let pollitem = - Zmq.Poll.mask_of + Zmq.Poll.mask_of [| (socket , Zmq.Poll.In) ; - (debug_socket , Zmq.Poll.In) + (debug_socket , Zmq.Poll.In) |] in while (!status <> Status.Stopped) do let polling = - Zmq.Poll.poll ~timeout:1000 pollitem + Zmq.Poll.poll ~timeout:1000 pollitem in if (polling.(0) = Some Zmq.Poll.In) then begin - let message = + let message = Zmq.Socket.recv_all ~block:false socket |> String.concat " " in @@ -530,7 +528,7 @@ let run ?(daemon=true) ezfio_filename = begin (* Forward subscription from XPUB to XSUB *) Zmq.Socket.recv_all ~block:false debug_socket - |> Zmq.Socket.send_all socket + |> Zmq.Socket.send_all socket end done; Zmq.Socket.set_linger_period socket 1000 ; @@ -539,17 +537,17 @@ let run ?(daemon=true) ezfio_filename = in (** {3 Main thread} *) - let random_walkers n_walks = + let random_walkers n_walks = let rec walkers accu = function - | 0 -> accu - | n -> + | 0 -> accu + | n -> let random_int = Random.int (Strictly_positive_int.to_int n_walks) in - let new_accu = + let new_accu = walkers_array.(random_int) :: accu in - walkers new_accu (n-1) + walkers new_accu (n-1) in walkers [] (Strictly_positive_int.to_int n_walks) |> Array.concat @@ -560,21 +558,21 @@ let run ?(daemon=true) ezfio_filename = let wall0 = Unix.gettimeofday () in - let f () = - + let f () = + change_status Status.Queued; send_log "status" 0 wall0 "Starting main thread"; (** Reply socket *) let rep_socket = Zmq.Socket.create zmq_context Zmq.Socket.rep - and address = + and address = Printf.sprintf "tcp://*:%d" port in bind_socket "REP" rep_socket address; Zmq.Socket.set_receive_high_water_mark rep_socket 100_000; Zmq.Socket.set_send_high_water_mark rep_socket 100_000; - Zmq.Socket.set_immediate rep_socket true; + Zmq.Socket.set_immediate rep_socket true; Zmq.Socket.set_linger_period rep_socket 600_000 ; (** EZFIO Cache *) @@ -595,14 +593,14 @@ let run ?(daemon=true) ezfio_filename = in List.iter (fun x -> if handle_ezfio ("has_"^x) = "T" then - try ignore @@ handle_ezfio ("get_"^x) + try ignore @@ handle_ezfio ("get_"^x) with Failure _ -> ()) Qptypes.all_ezfio_messages; (** Pull socket for computed data *) let pull_socket = Zmq.Socket.create zmq_context Zmq.Socket.pull - and address = + and address = Printf.sprintf "tcp://*:%d" (port+2) in bind_socket "PULL" pull_socket address; @@ -611,7 +609,7 @@ let run ?(daemon=true) ezfio_filename = (** Address of the dataserver *) let server_address = let ip = - Lazy.force Qmcchem_config.ip_address + Lazy.force Qmcchem_config.ip_address in Printf.sprintf "tcp://%s:%d" ip port in @@ -621,7 +619,7 @@ let run ?(daemon=true) ezfio_filename = (** Polling item to poll REP and PULL sockets. *) let pollitem = - Zmq.Poll.mask_of + Zmq.Poll.mask_of [| ( rep_socket, Zmq.Poll.In) ; ( pull_socket, Zmq.Poll.In) ; |] @@ -629,27 +627,27 @@ let run ?(daemon=true) ezfio_filename = (** Handles messages coming into the REP socket. *) - let handle_rep () = + let handle_rep () = let raw_msg = Zmq.Socket.recv_all ~block:false rep_socket in - let t0 = + let t0 = Unix.gettimeofday () in let msg = List.rev_map String.trim raw_msg |> List.rev - |> Message.create + |> Message.create and msg_size = List.fold_left (fun accu x -> accu + (String.length x)) 0 raw_msg in let handle = function | Message.Error _ -> () - | Message.Ezfio ezfio_msg -> - let result = + | Message.Ezfio ezfio_msg -> + let result = handle_ezfio ezfio_msg in - Zmq.Socket.send_all rep_socket + Zmq.Socket.send_all rep_socket [ String.length result |> Printf.sprintf "%d " ; result ] ; @@ -657,53 +655,53 @@ let run ?(daemon=true) ezfio_filename = | Message.GetWalkers n_walks -> begin send_log "req" msg_size t0 "get_walkers"; - let result = + let result = random_walkers n_walks in Zmq.Socket.send_all rep_socket result; send_log "rep" walkers_size t0 "get_walkers" - end - | Message.Register (w,pid) -> + end + | Message.Register (w,pid) -> begin match !status with - | Status.Queued + | Status.Queued | Status.Running -> begin - String.concat " " [ "Register :" ; - Compute_node.to_string w ; + String.concat " " [ "Register :" ; + Compute_node.to_string w ; string_of_int pid ] |> send_log "req" msg_size t0; add_worker w pid; if (!status = Status.Queued) then change_status Status.Running ; Zmq.Socket.send rep_socket "OK"; - send_log "rep" 2 t0 "Register : OK" + send_log "rep" 2 t0 "Register : OK" end - | Status.Stopping + | Status.Stopping | Status.Stopped -> Zmq.Socket.send rep_socket "Failed"; end - | Message.Unregister (w,pid) -> + | Message.Unregister (w,pid) -> begin String.concat " " [ "Unregister :" ; - (Compute_node.to_string w) ; + (Compute_node.to_string w) ; (string_of_int pid) ] |> send_log "req" msg_size t0; Zmq.Socket.send rep_socket "OK"; del_worker w pid; - String.concat " " [ "Unregister :"; + String.concat " " [ "Unregister :"; (Hashtbl.length workers_hash) |> string_of_int ; - "remaining" ] + "remaining" ] |> send_log "rep" 2 t0 ; - let n_connect = + let n_connect = n_connected workers_hash t0 in - match (daemon,n_connect) with + match (daemon,n_connect) with | (false,0) -> change_status Status.Stopped - | (true ,0) -> change_status Status.Queued + | (true ,0) -> change_status Status.Queued | _ -> () end - | Message.Test -> + | Message.Test -> begin Zmq.Socket.send rep_socket "OK"; send_log "rep" 2 t0 "Test" @@ -719,17 +717,17 @@ let run ?(daemon=true) ezfio_filename = let raw_msg = Zmq.Socket.recv_all ~block:false pull_socket in - let t0 = + let t0 = Unix.gettimeofday () in let msg = List.rev_map String.trim raw_msg |> List.rev - |> Message.create + |> Message.create and msg_size = List.fold_left (fun accu x -> accu + (String.length x)) 0 raw_msg in - let recv_log = + let recv_log = send_log "pull" msg_size t0 in @@ -739,7 +737,7 @@ let run ?(daemon=true) ezfio_filename = begin if (status = Status.Running) then touch_worker h pid ; - let log_msg = + let log_msg = Printf.sprintf "Walkers from %s : %d / %d / %d" (key h pid) (Array.length w) (!last_walker) walk_num_tot in @@ -754,13 +752,14 @@ let run ?(daemon=true) ezfio_filename = (Unix.gettimeofday () -. wall0) 1. (Property.to_string Property.Wall) hostname (string_of_int dataserver_pid) 1 - |> Block.of_string + |> Block.of_string in match wall with - | Some wall -> + | Some wall -> begin - output_string !block_channel (Block.to_string wall); - output_char !block_channel '\n'; + output_string !block_channel (Block.to_string_or_bytes wall); + if not Qmcchem_config.binary_io then + output_char !block_channel '\n'; end | _ -> () end @@ -768,16 +767,17 @@ let run ?(daemon=true) ezfio_filename = begin if (status = Status.Running) then touch_worker b.Block.compute_node b.Block.pid ; - output_string !block_channel (Block.to_string b); - output_char !block_channel '\n'; - recv_log (Block.to_string b) + output_string !block_channel (Block.to_string_or_bytes b); + if not Qmcchem_config.binary_io then + output_char !block_channel '\n'; + recv_log (Block.to_short_string b) end | Message.Test | Message.GetWalkers _ | Message.Ezfio _ | Message.Register (_, _) | Message.Unregister (_, _) - -> failwith "Bad message" + -> failwith "Bad message" in handle msg in @@ -785,18 +785,18 @@ let run ?(daemon=true) ezfio_filename = while (!status <> Status.Stopped) do let polling = - Zmq.Poll.poll ~timeout:1000 pollitem + Zmq.Poll.poll ~timeout:1000 pollitem in match polling.(1) with | Some Zmq.Poll.In -> handle_pull !status - | _ -> + | _ -> begin match polling.(0) with | Some Zmq.Poll.In -> handle_rep () - | _ -> + | _ -> begin flush !block_channel ; - let file_size = + let file_size = (Unix.stat block_channel_filename_locked).Unix.st_size in if (file_size > !max_file_size) then @@ -811,13 +811,13 @@ let run ?(daemon=true) ezfio_filename = List.iter (fun socket -> Zmq.Socket.set_linger_period socket 1000 ; Zmq.Socket.close socket) - [ rep_socket ; pull_socket ] + [ rep_socket ; pull_socket ] in Thread.create f in - + (** {2 Finalization} *) (** Cleans all the open files, sockets, etc. @@ -842,16 +842,16 @@ let run ?(daemon=true) ezfio_filename = (** {3 Main function} *) - let t0 = + let t0 = Unix.gettimeofday () in (* Handle signals *) - let handler s = + let handler s = Printf.printf "Dataserver received signal %d... killing\n%!" s; Watchdog.kill (); in - List.iter (fun s -> ignore @@ Sys.signal s (Sys.Signal_handle handler)) + List.iter (fun s -> ignore @@ Sys.signal s (Sys.Signal_handle handler)) [ Sys.sigint ; Sys.sigterm ; @@ -864,7 +864,7 @@ let run ?(daemon=true) ezfio_filename = begin try - (List.iter Thread.join + (List.iter Thread.join [ start_status_thread () ; start_log_thread () ; start_main_thread () ; diff --git a/ocaml/Random_variable.ml b/ocaml/Random_variable.ml index 0e9fafb..13482ab 100644 --- a/ocaml/Random_variable.ml +++ b/ocaml/Random_variable.ml @@ -1,13 +1,13 @@ open Qptypes -type t = +type t = { property : Property.t ; data : Block.t list; } module Average = struct - include Sample + include Sample end module Error = struct @@ -19,42 +19,42 @@ module Variance = struct end module Skewness: sig - type t + type t val to_float : t -> float val of_float : float -> t val to_string : t -> string end = struct - type t = float + type t = float let to_string = string_of_float let to_float x = x let of_float x = x end module Kurtosis: sig - type t + type t val to_float : t -> float val of_float : float -> t val to_string : t -> string end = struct - type t = float + type t = float let to_string = string_of_float let to_float x = x let of_float x = x end module GaussianDist: sig - type t + type t val create : mu:Average.t -> sigma2:Variance.t -> t val eval : g:t -> x:float -> float end = struct - type t = { mu: Average.t ; sigma2: Variance.t } + type t = { mu: Average.t ; sigma2: Variance.t } let create ~mu ~sigma2 = { mu ; sigma2 } let eval ~g ~x = - let { mu ; sigma2 } = + let { mu ; sigma2 } = g in - let mu = + let mu = Average.to_float mu and sigma2 = Variance.to_float sigma2 @@ -62,10 +62,10 @@ end = struct let x2 = (x -. mu) *. ( x -. mu) /. sigma2 in - let pi = + let pi = acos (-1.) in - let c = + let c = 1. /. (sqrt (sigma2 *. (pi +. pi))) in c *. exp ( -0.5 *. x2) @@ -78,7 +78,7 @@ let hashtbl_to_alist table = let hashtbl_change table key f = let elt = - try + try Some (Hashtbl.find table key) with | Not_found -> None @@ -91,7 +91,7 @@ let hashtbl_change table key f = (** Build from raw data. Range values are given in percent. *) let of_raw_data ?(locked=true) ~range property = - let data = + let data = Block.raw_data ~locked () |> List.filter (fun x -> x.Block.property = property) |> List.sort (fun x y -> @@ -109,7 +109,7 @@ let of_raw_data ?(locked=true) ~range property = (Weight.to_float x.Block.weight) +. accu ) 0. data in - + let wmin, wmax = rmin *. total_weight *. 0.01, rmax *. total_weight *. 0.01 @@ -128,13 +128,13 @@ let of_raw_data ?(locked=true) ~range property = (wsum_new, x::l) else (wsum_new, l) - end + end ) (0.,[]) data in List.rev new_data in - let result = + let result = match range with | (0.,100.) -> { property ; data } | (rmin,rmax) -> { property ; data=data_in_range rmin rmax } @@ -146,7 +146,7 @@ let of_raw_data ?(locked=true) ~range property = (** Compute average *) let average { property ; data } = if Property.is_scalar property then - let (num,denom) = + let (num,denom) = List.fold_left (fun (an, ad) x -> let num = (Weight.to_float x.Block.weight) *. (Sample.to_float x.Block.value) @@ -154,7 +154,7 @@ let average { property ; data } = (Weight.to_float x.Block.weight) in (an +. num, ad +. den) ) (0., 0.) data - in + in num /. denom |> Average.of_float else @@ -163,15 +163,15 @@ let average { property ; data } = | [] -> 1 | x :: tl -> Sample.dimension x.Block.value in - let (num,denom) = + let (num,denom) = List.fold_left (fun (an, ad) x -> let num = Array.map (fun y -> (Weight.to_float x.Block.weight) *. y) - (Sample.to_float_array x.Block.value) + (Sample.to_float_array x.Block.value) and den = (Weight.to_float x.Block.weight) - in ( Array.mapi (fun i y -> y +. num.(i)) an , ad +. den) + in ( Array.mapi (fun i y -> y +. num.(i)) an , ad +. den) ) (Array.make dim 0. , 0.) data - in + in let denom_inv = 1. /. denom in @@ -180,22 +180,28 @@ let average { property ; data } = - + (** Compute sum (for CPU/Wall time) *) let sum { property ; data } = List.fold_left (fun accu x -> let num = (Weight.to_float x.Block.weight) *. (Sample.to_float x.Block.value) in accu +. num - ) 0. data + ) 0. data (** Calculation of the average and error bar *) let ave_error { property ; data } = - let rec loop ~sum ~avsq ~ansum ~avsum ~n ?idx = function - | [] -> + (* sum: \sum_k x_k *. w_k + ansum: \sum_k w_k + avsum: \sum_k x_k *. w_k + avcu0: avsum / ansum + avsq: \sum_k (1. -. (w_k /. ansum_k)) *. (x_k -. avcu0)^2 *. w_k) + *) + let rec loop ~sum ~avsq ~ansum ~avsum ~n ?idx = function + | [] -> begin if (n > 0.) then ( Average.of_float (sum /. ansum), @@ -205,12 +211,8 @@ let ave_error { property ; data } = end | (x,w) :: tail -> begin - let avcu0 = - avsum /. ansum - in - let xw = - x *. w - in + let avcu0 = avsum /. ansum in + let xw = x *. w in let ansum, avsum, sum = ansum +. w , avsum +. xw , @@ -220,9 +222,9 @@ let ave_error { property ; data } = ~sum:sum ~avsq:(avsq +. (1. -. (w /. ansum)) *. (x -. avcu0) *. (x -. avcu0) *. w) - ~avsum:avsum + ~avsum:avsum ~ansum:ansum - ~n:(n +. 1.) + ~n:(n +. 1.) end in @@ -242,7 +244,7 @@ let ave_error { property ; data } = (Sample.to_float x.Block.value, Weight.to_float x.Block.weight) ) data - |> ave_error_scalar + |> ave_error_scalar else match data with | [] -> (Average.of_float 0., None) @@ -251,7 +253,7 @@ let ave_error { property ; data } = head.Block.value |> Sample.dimension in - let result = + let result = Array.init dim (fun idx -> List.rev_map (fun x -> (Sample.to_float ~idx x.Block.value, @@ -260,16 +262,16 @@ let ave_error { property ; data } = |> ave_error_scalar ) in - ( Array.map (fun (x,_) -> Average.to_float x) result - |> Average.of_float_array ~dim , + ( Array.map (fun (x,_) -> Average.to_float x) result + |> Average.of_float_array ~dim , if (Array.length result < 2) then None else - Some (Array.map (function + Some (Array.map (function | (_,Some y) -> Error.to_float y - | (_,None) -> 0.) result + | (_,None) -> 0.) result |> Average.of_float_array ~dim) - ) + ) @@ -286,14 +288,14 @@ let fold_blocks ~f { property ; data } = List.fold_left (fun accu block -> let x = Sample.to_float block.Block.value in f accu x - ) init data + ) init data (** Convergence plot *) let convergence { property ; data } = - let rec loop ~sum ~avsq ~ansum ~avsum ~n ~accu = function + let rec loop ~sum ~avsq ~ansum ~avsum ~n ~accu = function | [] -> List.rev accu | head :: tail -> begin @@ -307,7 +309,7 @@ let convergence { property ; data } = and avsum = avsum +. xw and sum = sum +. xw in - let accu = + let accu = if (n > 0.) then (sum /. ansum, sqrt ( abs_float ( avsq /.( ansum *. n))))::accu else @@ -317,9 +319,9 @@ let convergence { property ; data } = ~sum:sum ~avsq:(avsq +. (1. -. (w /. ansum)) *. (x -. avcu0) *. (x -. avcu0) *. w) - ~avsum:avsum + ~avsum:avsum ~ansum:ansum - ~n:(n +. 1.) + ~n:(n +. 1.) ~accu:accu end in @@ -365,7 +367,7 @@ let max_block = (** Create a hash table for merging *) -let create_hash ~create_key ?(update_block_id=(fun x->x)) +let create_hash ~create_key ?(update_block_id=(fun x->x)) ?(update_value=(fun wc vc wb vb sw -> (wc *. vc +. wb *. vb) /. sw) ) ?(update_weight=(fun wc wb -> wc +. wb) ) t = let table = Hashtbl.create 63 @@ -374,7 +376,7 @@ let create_hash ~create_key ?(update_block_id=(fun x->x)) let key = create_key block in let open Block in - hashtbl_change table key (function + hashtbl_change table key (function | Some current -> let wc, wb = Weight.to_float current.weight, @@ -384,7 +386,7 @@ let create_hash ~create_key ?(update_block_id=(fun x->x)) update_weight wc wb in if (Property.is_scalar current.property) then - let vc, vb = + let vc, vb = Sample.to_float current.value, Sample.to_float block.value in Some @@ -396,15 +398,15 @@ let create_hash ~create_key ?(update_block_id=(fun x->x)) compute_node = block.compute_node; } else - let vc, vb = + let vc, vb = Sample.to_float_array current.value, Sample.to_float_array block.value - and dim = + and dim = Sample.dimension current.value in Some { property = current.property ; weight = Weight.of_float sw ; - value = + value = Array.init dim (fun i -> update_value wc vc.(i) wb vb.(i) sw) |> Sample.of_float_array ~dim ; block_id = update_block_id block.block_id; @@ -443,15 +445,15 @@ let merge ~create_key ?update_block_id ?update_value ?update_weight t = (** Merge per block id *) let merge_per_block_id = - merge + merge ~create_key:(fun block -> Block_id.to_string block.Block.block_id) (** Merge per compute_node *) let merge_per_compute_node = merge - ~create_key:(fun block -> - Printf.sprintf "%s" + ~create_key:(fun block -> + Printf.sprintf "%s" (Compute_node.to_string block.Block.compute_node) ) @@ -459,8 +461,8 @@ let merge_per_compute_node = (** Merge per Compute_node and PID *) let merge_per_compute_node_and_pid = merge - ~create_key:(fun block -> - Printf.sprintf "%s %10.10d" + ~create_key:(fun block -> + Printf.sprintf "%s %10.10d" (Compute_node.to_string block.Block.compute_node) (block.Block.pid) ) @@ -469,8 +471,8 @@ let merge_per_compute_node_and_pid = (** Merge per Compute_node and BlockId *) let merge_per_compute_node_and_block_id = merge - ~create_key:(fun block -> - Printf.sprintf "%s %10.10d" + ~create_key:(fun block -> + Printf.sprintf "%s %10.10d" (Compute_node.to_string block.Block.compute_node) (Block_id.to_int block.Block.block_id) ) @@ -510,48 +512,48 @@ let error_x_over_y = function (** Create float, variable operators *) -let one_variable_operator ~update_value p f = - { p with - data = List.rev @@ List.rev_map (fun b -> { b with +let one_variable_operator ~update_value p f = + { p with + data = List.rev @@ List.rev_map (fun b -> { b with Block.value = Sample.of_float (update_value (Sample.to_float b.Block.value) ) } ) p.data } let ( +@ ) p f = one_variable_operator p f ~update_value: (fun x -> x +. f ) -let ( *@ ) p f = one_variable_operator p f +let ( *@ ) p f = one_variable_operator p f ~update_value: (fun x -> x *. f ) -let ( -@ ) p f = one_variable_operator p f +let ( -@ ) p f = one_variable_operator p f ~update_value: (fun x -> x -. f ) -let ( /@ ) p f = one_variable_operator p f +let ( /@ ) p f = one_variable_operator p f ~update_value: (fun x -> x /. f ) (** Create two variable operators *) -let two_variable_operator ~update_value p1 p2 = +let two_variable_operator ~update_value p1 p2 = merge ~update_value ~create_key:(fun block -> - Printf.sprintf "%s %10.10d %10.10d" + Printf.sprintf "%s %10.10d %10.10d" (Compute_node.to_string block.Block.compute_node) - (Block_id.to_int block.Block.block_id) + (Block_id.to_int block.Block.block_id) (block.Block.pid) ) ~update_weight:(fun wc wb -> wc ) { property = p1.property ; - data = List.concat [ p1.data ; p2.data ] } + data = List.concat [ p1.data ; p2.data ] } -let ( +! ) = two_variable_operator +let ( +! ) = two_variable_operator ~update_value: (fun wc vc wb vb sw -> (vc +. vb) ) -let ( *! ) = two_variable_operator +let ( *! ) = two_variable_operator ~update_value: (fun wc vc wb vb sw -> (vc *. vb) ) -let ( -! ) = two_variable_operator +let ( -! ) = two_variable_operator ~update_value: (fun wc vc wb vb sw -> (vc -. vb) ) -let ( /! ) = two_variable_operator +let ( /! ) = two_variable_operator ~update_value: (fun wc vc wb vb sw -> (vc /. vb) ) @@ -560,11 +562,11 @@ let ( /! ) = two_variable_operator (** Merge two consecutive blocks *) let compress = merge - ~create_key:(fun block -> + ~create_key:(fun block -> Printf.sprintf "%s %10.10d %10.10d" (Compute_node.to_string block.Block.compute_node) block.Block.pid (((Block_id.to_int block.Block.block_id)+1)/2)) - ~update_block_id:(fun block_id -> + ~update_block_id:(fun block_id -> ((Block_id.to_int block_id)+1)/2 |> Block_id.of_int ) @@ -576,15 +578,15 @@ let max_value_per_compute_node t = let table = Hashtbl.create 63 in let create_key block = - Printf.sprintf "%s %10.10d" + Printf.sprintf "%s %10.10d" (Compute_node.to_string block.Block.compute_node) - (block.Block.pid) + (block.Block.pid) in List.iter (fun block -> let key = create_key block in let open Block in - hashtbl_change table key (function + hashtbl_change table key (function | Some current -> let vc = Sample.to_float current.value and vb = Sample.to_float block.value @@ -610,36 +612,36 @@ let max_value_per_compute_node t = (** String representation *) -let to_string p = +let to_string p = match p.property with | Property.Cpu -> Printf.sprintf "%s" (Time.string_of_sec (sum p)) | Property.Wall -> Printf.sprintf "%s" (Time.string_of_sec (sum (max_value_per_compute_node p))) | Property.Accep -> Printf.sprintf "%16.10f" (average p |> Average.to_float) - | _ -> + | _ -> begin if Property.is_scalar p.property then match ave_error p with - | (ave, Some error) -> - let (ave, error) = - Average.to_float ave, + | (ave, Some error) -> + let (ave, error) = + Average.to_float ave, Error.to_float error in Printf.sprintf "%16.10f +/- %16.10f" ave error - | (ave, None) -> - let ave = + | (ave, None) -> + let ave = Average.to_float ave in Printf.sprintf "%16.10f" ave else match ave_error p with - | (ave, Some error) -> + | (ave, Some error) -> let idxmax = Average.dimension ave in let rec f accu idx = if (idx < idxmax) then - let (ave, error) = - Average.to_float ~idx ave, + let (ave, error) = + Average.to_float ~idx ave, Error.to_float ~idx error in let s = @@ -650,9 +652,9 @@ let to_string p = accu in (f "[ \n" 0) ^ " ]" - | (ave, None) -> + | (ave, None) -> Average.to_float ave - |> Printf.sprintf "%16.10f" + |> Printf.sprintf "%16.10f" end @@ -666,19 +668,19 @@ let compress_files () = let properties = Lazy.force Block.properties in - + (* Create temporary file *) let dir_name = Block.dir_name in - let dir_name = + let dir_name = Lazy.force dir_name in - let files = + let files = Sys.readdir dir_name |> Array.to_list - |> List.filter (fun x -> + |> List.filter (fun x -> try Str.search_backward (Str.regexp "locked") x (String.length x) >= 0 with @@ -688,10 +690,10 @@ let compress_files () = |> List.rev in - let out_channel_dir = + let out_channel_dir = let rand_num = Random.int 1000000 |> string_of_int in - let dirname = - Filename.concat !Ezfio.ezfio_filename "blocks" + let dirname = + Filename.concat !Ezfio.ezfio_filename "blocks" in if not ( Sys.file_exists dirname ) then Unix.mkdir dirname 0o755; @@ -706,31 +708,31 @@ let compress_files () = raise (Sys_error message) in - let out_channel_name = + let out_channel_name = let hostname = - Lazy.force Qmcchem_config.hostname + Lazy.force Qmcchem_config.hostname and suffix = Unix.getpid () - |> string_of_int + |> string_of_int in String.concat "." [ hostname ; suffix ] in - let block_channel = + let block_channel = Filename.concat out_channel_dir out_channel_name |> open_out in List.iter (fun p -> - let l = + let l = match p with - | Property.Cpu + | Property.Cpu | Property.Accep -> of_raw_data ~locked:false ~range:(0.,100.) p |> merge_per_compute_node - | Property.Wall -> + | Property.Wall -> of_raw_data ~locked:false ~range:(0.,100.) p - |> max_value_per_compute_node + |> max_value_per_compute_node | _ -> of_raw_data ~locked:false ~range:(0.,100.) p (* @@ -738,10 +740,11 @@ let compress_files () = *) in List.iter (fun x -> - output_string block_channel (Block.to_string x); - output_char block_channel '\n'; - ) l.data - ) properties ; + output_string block_channel (Block.to_string_or_bytes x); + if not Qmcchem_config.binary_io then + output_char block_channel '\n'; + ) l.data + ) properties ; close_out block_channel; List.iter Unix.unlink files ; @@ -760,17 +763,17 @@ let autocovariance { property ; data } = match (merge_per_block_id { property ; data }) with { property ; data } -> Array.of_list data in - let x_t = + let x_t = Array.map (fun x -> (Sample.to_float x.Block.value) -. ave) data in - let f i = - let denom = + let f i = + let denom = if (i > 1) then (float_of_int i) else 1. in - let r = + let r = Array.sub x_t 0 i |> Array.fold_left (fun accu x -> - accu +. x *. x_t.(i)) 0. + accu +. x *. x_t.(i)) 0. in r /. denom in @@ -786,14 +789,14 @@ let centered_cumulants { property ; data } = |> Average.to_float in let centered_data = - List.rev_map (fun x -> - ( (Weight.to_float x.Block.weight), + List.rev_map (fun x -> + ( (Weight.to_float x.Block.weight), (Sample.to_float x.Block.value) -. ave ) ) data |> List.rev in - let var = - let (num, denom) = + let var = + let (num, denom) = List.fold_left (fun (a2, ad) (w,x) -> let x2 = x *. x in @@ -802,18 +805,18 @@ let centered_cumulants { property ; data } = in (a2 +. var, ad +. den) ) (0., 0.) centered_data in num /. denom - in + in let centered_data = let sigma_inv = 1. /. (sqrt var) in - List.rev_map (fun x -> - ( (Weight.to_float x.Block.weight), + List.rev_map (fun x -> + ( (Weight.to_float x.Block.weight), ( (Sample.to_float x.Block.value) -. ave ) *. sigma_inv ) ) data |> List.rev in - let (cum3,cum4) = + let (cum3,cum4) = let (cum3, cum4, denom) = List.fold_left (fun (a3, a4, ad) (w,x) -> let x2 = x *. x @@ -823,9 +826,9 @@ let centered_cumulants { property ; data } = and den = w in (a3 +. cum3, a4 +. cum4, ad +. den) ) (0., 0., 0.) centered_data - in + in ( cum3 /. denom, cum4 /. denom -. 3. ) - in + in [| ave ; var ; cum3 ; cum4 |] @@ -833,26 +836,26 @@ let centered_cumulants { property ; data } = (** Computes a histogram *) let histogram { property ; data } = - let min, max = + let min, max = (min_block { property ; data }), (max_block { property ; data }) in - let length = + let length = max -. min and n = List.length data |> float_of_int |> sqrt in - let delta_x = + let delta_x = length /. (n-.1.) and result = - Array.init (int_of_float n + 1) (fun _ -> 0.) + Array.init (int_of_float n + 1) (fun _ -> 0.) in List.iter (fun x -> let w = (Weight.to_float x.Block.weight) - and x = + and x = (Sample.to_float x.Block.value) in let i = @@ -862,7 +865,7 @@ let histogram { property ; data } = result.(i) <- result.(i) +. w ) data ; - let norm = + let norm = 1. /. ( delta_x *. ( Array.fold_left (fun accu x -> accu +. x) 0. result ) ) diff --git a/ocaml/Sample.ml b/ocaml/Sample.ml index 540cb8b..f1f6c51 100644 --- a/ocaml/Sample.ml +++ b/ocaml/Sample.ml @@ -1,6 +1,6 @@ open Sexplib.Std -type t = +type t = | One_dimensional of float | Multidimensional of (float array * int) [@@deriving sexp] @@ -9,27 +9,27 @@ let dimension = function | One_dimensional _ -> 1 | Multidimensional (_,d) -> d -let to_float ?idx x = +let to_float ?idx x = match (idx,x) with | None , One_dimensional x - | Some 0, One_dimensional x -> x + | Some 0, One_dimensional x -> x | Some i, One_dimensional x -> failwith "Index should not be specified in One_dimensional" | None , Multidimensional (x,_) -> x.(0) | Some i, Multidimensional (x,s) when i < s -> x.(i) - | Some i, Multidimensional (x,s) -> - Printf.sprintf "Index out of bounds in Multidimensional + | Some i, Multidimensional (x,s) -> + Printf.sprintf "Index out of bounds in Multidimensional %d not in [0,%d[ " i s |> failwith -let to_float_array = function +let to_float_array = function | One_dimensional _ -> failwith "Should be Multidimensional" | Multidimensional (x,_) -> x -let of_float x = +let of_float x = One_dimensional x -let of_float_array ~dim x = +let of_float_array ~dim x = if (Array.length x) <> dim then failwith "Inconsistent array size in of_float_array" else @@ -39,9 +39,26 @@ let of_float_array ~dim x = let to_string = function | One_dimensional x -> string_of_float x - | Multidimensional (x,_) -> + | Multidimensional (x,_) -> Array.map string_of_float x |> Array.to_list - |> String.concat " " - |> Printf.sprintf "%s" + |> String.concat " " +let to_bytes = function + | One_dimensional x -> Qptypes.bytes_of_float x + | Multidimensional (x,_) -> + let b = Bytes.create (8 * Array.length x) in + Array.iteri (fun i x -> + Int64.bits_of_float x + |> Bytes.set_int64_le b (i*8) ) x; + b + +let of_bytes b = + match Bytes.length b with + | 8 -> let x = Qptypes.float_of_bytes b in + One_dimensional x + | l -> let len = l/8 in + Multidimensional ( Array.init len (fun i -> + Bytes.get_int64_le b (i*8) + |> Int64.float_of_bits ), + len ) diff --git a/ocaml/Sample.mli b/ocaml/Sample.mli index 58a1f64..0a2a545 100644 --- a/ocaml/Sample.mli +++ b/ocaml/Sample.mli @@ -2,7 +2,9 @@ type t [@@deriving sexp] val to_float : ?idx:int -> t -> float val to_float_array : t -> float array val of_float : float -> t -val of_float_array : dim:int -> float array -> t +val of_float_array : dim:int -> float array -> t val to_string : t -> string +val to_bytes : t -> bytes +val of_bytes : bytes -> t val dimension : t -> int diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index e57b9f4..3e93818 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -4,7 +4,15 @@ let global_replace x = |> Str.global_replace (Str.regexp "Float.of_string") "float_of_string" |> Str.global_replace (Str.regexp "Int.to_string") "string_of_int" |> Str.global_replace (Str.regexp "Int.of_string") "int_of_string" - |> Str.global_replace (Str.regexp "String.\\(to\\|of\\)_string") "" + |> Str.global_replace (Str.regexp "Int.to_bytes") "bytes_of_int" + |> Str.global_replace (Str.regexp "Int64.to_bytes") "bytes_of_int64" + |> Str.global_replace (Str.regexp "Float.to_bytes") "bytes_of_float" + |> Str.global_replace (Str.regexp "Float.of_bytes") "float_of_bytes" + |> Str.global_replace (Str.regexp "Int.of_bytes") "int_of_bytes" + |> Str.global_replace (Str.regexp "Int64.of_bytes") "int64_of_bytes" + |> Str.global_replace (Str.regexp "String.\\(to\\|of\\)_string") "" + |> Str.global_replace (Str.regexp "String.to_bytes") "Bytes.of_string" + |> Str.global_replace (Str.regexp "String.of_bytes") "Bytes.to_string" let input_data = " * Positive_float : float @@ -38,7 +46,7 @@ let input_data = " * Negative_int : int if not (x <= 0) then raise (Invalid_argument (Printf.sprintf \"Negative_int : (x <= 0) : x=%d\" x)); - assert (x <= 0) ; + assert (x <= 0) ; * Det_coef : float if (x < -1.) || (x > 1.) then @@ -167,6 +175,34 @@ let input_ezfio = " let untouched = " +let bytes_of_int64 i = + let result = Bytes.create 8 in + Bytes.set_int64_le result 0 i; + result + +let bytes_of_int i = + Int64.of_int i + |> bytes_of_int64 + + +let int64_of_bytes b = + Bytes.get_int64_le b 0 + + +let int_of_bytes b = + int64_of_bytes b + |> Int64.to_int + + +let float_of_bytes b = + int64_of_bytes b + |> Int64.float_of_bits + + +let bytes_of_float f = + Int64.bits_of_float f + |> bytes_of_int64 + " let template = format_of_string " @@ -175,11 +211,15 @@ module %s : sig val to_%s : t -> %s val of_%s : %s %s -> t val to_string : t -> string + val to_bytes : t -> bytes + val of_bytes : bytes -> t end = struct type t = %s [@@deriving sexp] let to_%s x = x let of_%s %s x = ( %s x ) let to_string x = %s.to_string x + let to_bytes x = %s.to_bytes x + let of_bytes b = %s.of_bytes b end " @@ -203,7 +243,7 @@ let parse_input input= and name = String_ext.strip name in let typ_cap = String.capitalize_ascii typ in let newstring = Printf.sprintf template name typ typ typ params_val typ typ - typ typ params ( String_ext.strip text ) typ_cap + typ typ params ( String_ext.strip text ) typ_cap typ_cap typ_cap in List.rev (parse (newstring::result) tail ) in @@ -223,9 +263,11 @@ module %s : sig val get_max : unit -> %s val of_%s : ?min:%s -> ?max:%s -> %s -> t val to_string : t -> string + val to_bytes : t -> bytes end = struct type t = %s [@@deriving sexp] let to_string x = %s.to_string x + let to_bytes x = %s.to_bytes x let get_max () = if (Ezfio.has_%s ()) then Ezfio.get_%s () @@ -251,6 +293,10 @@ end = struct end " +(* + val of_bytes : bytes -> t + let of_bytes x = %s.of_bytes x +*) let parse_input_ezfio input= let parse s = @@ -271,9 +317,10 @@ let parse_input_ezfio input= | [ name ; typ ; ezfio_func ; min ; max ; msg ] -> (name, typ, ezfio_func, min, max, msg) | _ -> assert false in + let typ_cap = String.capitalize_ascii typ in Printf.sprintf ezfio_template - name typ typ typ typ typ typ typ typ (String.capitalize_ascii typ) - ezfio_func ezfio_func max min typ typ max msg min name (String.capitalize_ascii typ) + name typ typ typ typ typ typ typ typ typ_cap typ_cap + ezfio_func ezfio_func max min typ typ max msg min name typ_cap end | _ -> failwith "Error in input_ezfio" in @@ -294,29 +341,30 @@ let input_lines filename = let create_ezfio_handler () = - let lines = + let lines = input_lines "ezfio.ml" - |> List.mapi (fun i l -> if i > 417 then Some l else None) + (* /!\ Change when ezfio.ml changes *) + |> List.mapi (fun i l -> if i > 442 then Some l else None) |> List.filter (fun x -> x <> None) |> List.map (fun x -> match x with | Some x -> x | None -> assert false) in - let functions = + let functions = List.map (fun x -> match String.split_on_char ' ' x with | _ :: x :: "()" :: "=" :: f :: dir :: item :: _-> (x, f, dir, item) | _ :: x :: "=" :: f :: dir :: item :: _-> (x, f, dir, item) | _ -> ("","","","") - ) lines + ) lines in - let has_functions = - List.filter (fun (x,_,_,_) -> String.sub x 0 4 = "has_") functions - and get_functions = + let has_functions = + List.filter (fun (x,_,_,_) -> String.sub x 0 4 = "has_") functions + and get_functions = List.filter (fun (x,_,_,_) -> String.sub x 0 4 = "get_") functions in - let chop s = + let chop s = match (Str.split_delim (Str.regexp ";;") s) with | x :: _ -> x | _ -> assert false @@ -329,17 +377,17 @@ match msg with " ] @ List.map (fun (x,f,d,i) -> let i = chop i in if (String.sub f ((String.length f)-6) 6 = "_array") then - Printf.sprintf " | \"%s\" -> + Printf.sprintf " | \"%s\" -> Ezfio.read_string_array %s %s - |> Ezfio.flattened_ezfio + |> Ezfio.flattened_ezfio |> Array.to_list |> String.concat \" \"" x d i else Printf.sprintf " | \"%s\" -> Ezfio.read_string %s %s" x d i - ) get_functions + ) get_functions ) @ ( List.map (fun (x,_,_,_) -> - Printf.sprintf " | \"%s\" -> if (Ezfio.%s ()) then \"T\" else \"F\"" x x + Printf.sprintf " | \"%s\" -> if (Ezfio.%s ()) then \"T\" else \"F\"" x x ) has_functions ) @ [" | x -> failwith (x^\" : Unknown EZFIO function\")\n;;" ; @@ -347,17 +395,17 @@ match msg with " ] @ ( List.rev_map (fun (x,_,_,_) -> Printf.sprintf " \"%s\" ; " (String.sub x 4 ((String.length x)-4)) - ) has_functions + ) has_functions ) @ ["]"] in - String.concat "\n" result + String.concat "\n" result |> print_endline (** Main *) let () = + print_endline untouched; parse_input input_data ; parse_input_ezfio input_ezfio; - print_endline untouched; create_ezfio_handler () diff --git a/scripts/compile_src.sh b/scripts/compile_src.sh index f9f9cab..e195426 100755 --- a/scripts/compile_src.sh +++ b/scripts/compile_src.sh @@ -18,12 +18,6 @@ OBJ="${OBJ} IRPF90_temp/ZMQ/f77_zmq_module.o" INCLUDES="${INCLUDES} -I AO -I SAMPLING -I TOOLS -I JASTROW -I MAIN -I PROPERTIES -I ZMQ" IRPF90_FLAGS="${IRPF90_FLAGS} ${INCLUDES}" -# Check IRPF90 version -if [[ $( ${IRPF90} -v | python2 -c "import sys ; print float(sys.stdin.read().rsplit('.',1)[0]) >= 1.6") == False ]] -then - echo "IRPF90 version >= 1.6 required" - exit -1 -fi export IRPF90 IRPF90_FLAGS INCLUDES LIB SRC OBJ diff --git a/scripts/create_properties_ezfio.py b/scripts/create_properties_ezfio.py index ed82d8a..e7fb388 100755 --- a/scripts/create_properties_ezfio.py +++ b/scripts/create_properties_ezfio.py @@ -145,6 +145,14 @@ let to_string = function for p in properties_qmcvar: print >>file, """| %(P)s -> "%(P)s" """%{'P':p[1].capitalize(), 'p':p[1]} print >>file, """;; + +let of_bytes x = + Bytes.to_string x + |> of_string + +let to_bytes x = + to_string x + |> Bytes.of_string """ # is_scalar diff --git a/src/MAIN/admc.org b/src/MAIN/admc.org new file mode 100644 index 0000000..ba4432c --- /dev/null +++ b/src/MAIN/admc.org @@ -0,0 +1,277 @@ +#+TITLE: Asynchronous DMC +#+AUTHOR: Anthony Scemama +#+EMAIL: scemama@irsamc.ups-tlse.fr + + +#+PROPERTY: header-args :tangle no :noweb yes + + +* Main program + +** Declarations + + #+NAME: declarations + #+begin_src f90 + include '../types.F' + + integer :: iter + integer :: k_full ! Index of walkers in elec_coord_full + integer :: k_new ! Index of walkers in elec_coord_new + integer :: iw ! Number of copies in branching + integer :: l + + real, allocatable :: elec_coord_new(:,:,:) + + double precision :: w + double precision, allocatable :: E_out(:), w_sum(:) + + double precision, external :: qmc_ranf + + allocate(elec_coord_new(elec_num+1,3,walk_num)) + allocate(E_out(walk_num), w_sum(walk_num)) + #+end_src + +** Main flow + + - Fetch ~walk_num~ electron coordinates in ~elec_coord_full~ + - For each set of coordinates, + - Make a PDMC trajectory, and output the weight ~w~ + - Perform branching depending on the value of the weight + - Store the new sets of coordinates in ~elec_coord_new~ + - When ~elec_coord_new~ is full, send it to the server + + #+begin_src f90 :tangle "admc.irp.f" +program admc + call run + call ezfio_finish +end program admc + +subroutine run + implicit none + + <> + + ! Initialization + if (vmc_algo /= t_Brownian) then + call abrt(irp_here,'DMC should run with Brownian algorithm') + endif + + do iter=1,1000 +! call read_coords() + k_new = 1 + do k_full=1,walk_num + call pdmc_trajectory(k_full, w, E_out(k_full), w_sum(k_full)) + + <> + + if (k_new >= walk_num) then + w_sum(k_full+1:) = 0.d0 + exit + end if + end do + + k_new = k_new-1 + elec_coord_full(1:elec_num+1,1:3,1:k_new) = & + elec_coord_new(1:elec_num+1,1:3,1:k_new) + +! call write_coords(k_new) + call write_energy(walk_num, E_out, w_sum) + end do + +end subroutine run + +<> +<> +<> +<> + #+end_src + +** Branching + + #+NAME: branching + #+begin_src f90 + ! Find number of copies + iw = int(w) + w = w - int(w) + if (qmc_ranf() < w) then + iw = iw+1 + end if + + ! Duplicate walker + do l=1,iw + elec_coord_new(1:elec_num+1,1:3,k_new) = & + elec_coord(1:elec_num+1,1:3) + k_new = k_new+1 + if (k_new >= walk_num) exit + end do + + #+end_src + +* Read/write + +** Read coordinates + + Fetch a new set of coordinates for ~walk_num~ walkers from the pool of coordinates. + + #+NAME: read_coords + #+begin_src f90 +subroutine read_coords() + implicit none + + integer :: i, k + + do k=1,walk_num + do i=1,elec_num + read(*,*) elec_coord_full(i,1:3,k) + end do + end do + + SOFT_TOUCH elec_coord_full + +end subroutine read_coords + #+end_src + +** Write coordinates + + Send the current set of coordinates for ~walk_num~ walkers to the pool of coordinates. + + #+NAME: write_coords + #+begin_src f90 +subroutine write_coords() + implicit none + + integer :: i, k + + do k=1,walk_num + do i=1,elec_num + write(*,*) 'C', elec_coord_full(i,1:3,k) + end do + end do + +end subroutine write_coords + #+end_src + +** Write energy + + Compute the weighted average over the computed energies. + \[ + E = \frac{\sum_i w_i E_i}{\sum_i w_i} + \] + + #+NAME: write_energy + #+begin_src f90 +subroutine write_energy(walk_num_, E_out, w_sum) + implicit none + integer, intent(in) :: walk_num_ + double precision, intent(in) :: E_out(walk_num_) + double precision, intent(in) :: w_sum(walk_num_) + + integer :: i, k + double precision :: E, S + + E = 0.d0 + S = 0.d0 + do k=1,walk_num + S = S + w_sum(k) + E = E + w_sum(k) * E_out(k) + end do + write(*,*) 'E', E/S, S + +end subroutine write_energy + #+end_src + +* PDMC trajectory + + Computes a PDMC trajectory until the weight ~w~ is $1/2 < w < 3/2$. + The energy of the trajectory is computed as + \[ + E = \frac{\sum_i w_i E(R_i)}{\sum_i w_i} + \] + + The function returns: + - ~w~: the last of all $w_i$ + - ~E_out~: The average energy $E$ of the trajectory + - ~w_sum~: The sum of the weights + + #+NAME: declarations_pdmc + #+begin_src f90 + integer :: i,j,l + double precision :: delta + + ! If true, continue to make more steps + logical :: loop + + ! Max number of steps + integer :: imax + integer, parameter :: nmax=10000 + + ! Brownian step variables + double precision :: p,q + real :: delta_x + logical :: accepted + + ! Local energies from the past + double precision :: E_loc_save(4) + double precision :: w + + #+end_src + + #+NAME: pdmc_trajectory + #+begin_src f90 +subroutine pdmc_trajectory(k_full, pdmc_weight, E_out, w_sum) + implicit none + + integer, intent(in) :: k_full + double precision, intent(out) :: pdmc_weight, E_out, w_sum + + <> + + elec_coord(1:elec_num+1,1:3) = elec_coord_full(1:elec_num+1,1:3,k_full) + TOUCH elec_coord + + E_out = 0.d0 + w_sum = 0.d0 + + E_loc_save(1:4) = E_loc + + pdmc_weight = 1.d0 + loop = .True. + + do imax = 1, nmax + + call brownian_step(p,q,accepted,delta_x) + +! delta = (9.d0*E_loc+19.d0*E_loc_save(1)-5.d0*E_loc_save(2)+E_loc_save(3))/24.d0 + delta = E_loc + delta = (delta - E_ref)*p + + if (delta >= 0.d0) then + w = dexp(-dtime_step*delta) + else + w = 2.d0-dexp(dtime_step*delta) + endif + elec_coord(elec_num+1,1) += p*time_step + elec_coord(elec_num+1,2) = E_loc + elec_coord(elec_num+1,3) = pdmc_weight + if (accepted) then + E_loc_save(4) = E_loc_save(3) + E_loc_save(3) = E_loc_save(2) + E_loc_save(2) = E_loc_save(1) + E_loc_save(1) = E_loc + endif + + w_sum = w_sum + pdmc_weight + E_out = E_out + pdmc_weight * E_loc + + pdmc_weight = pdmc_weight * w + + loop = pdmc_weight > 0.5d0 .and. pdmc_weight < 2.0d0 + if (.not.loop) exit + + end do + + E_out = E_out / w_sum + +end subroutine pdmc_trajectory + #+end_src + diff --git a/src/MAIN/admc.py b/src/MAIN/admc.py new file mode 100755 index 0000000..f9ecf2b --- /dev/null +++ b/src/MAIN/admc.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python3 + +from mpi4py import MPI +import sys +import gzip +import random +import math +import subprocess + +admc_exec = "/home/scemama/qmcchem/src/MAIN/admc" +n_walk_per_proc = 10 + +def start(): + return subprocess.Popen( + [ admc_exec, sys.argv[1] ], + stdin=subprocess.PIPE, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE) + + +def read(process,len_walk): + line = process.stdout.readline().decode("utf-8").strip() + walk_num = int(line) + walkers = [] + print(walk_num) + for k in range(walk_num): + w = [] + for i in range(len_walk): + line = process.stdout.readline().decode("utf-8").strip() + w.append( line ) + w = '\n'.join(w) + walkers.append(w) + + _, E, W = process.stdout.readline().decode("utf-8").split() + return walkers, float(E), float(W) + + +def write(process, message): + process.stdin.write(f"{message}\n".encode("utf-8")) + process.stdin.flush() + + +def terminate(process): + process.stdin.close() + process.terminate() + process.wait(timeout=0.2) + +def print_energy(EnergyWeight, Energy2Weight, Weight, N): + e = EnergyWeight / Weight + e2 = Energy2Weight / Weight + err = math.sqrt(abs(e*e - e2) / max(1,(N-1)) ) + print("%f +/- %f"%(e, err)) + return err + +def main(): + try: + input_dir = sys.argv[1] + except: + print("syntax: argv[0] [FILE]") + sys.exit(-1) + + # Pool of electron coordinates + with gzip.open(input_dir+"/electrons/elec_coord_pool.gz","r") as f: + data = f.read().decode("utf-8").split() + + len_walk = int(data[1])*int(data[2]) + icount = 0 + buffer = [] + walkers = [] + for d in data[4:]: + buffer.append(d) + icount += 1 + if (icount == len_walk): + walkers.append(buffer) + buffer = [] + icount = 0 + + walkers = [ '\n'.join(x) for x in walkers ] + do_loop = True + + EnergyWeight = 0. + Energy2Weight = 0. + Weight = 0. + NSamples = 0. + + # Start processes + proc = start() + while do_loop: + + # Once every 1000, shuffle the list of walkers + if random.random() < 0.01: + print("SHUFFLE") + random.shuffle(walkers) + + # Pick new walkers + new_coords = walkers[:n_walk_per_proc] + walkers = walkers[n_walk_per_proc:] + + # Send new walkers to the process + write(proc, '\n'.join(new_coords)) + + # Fetch new walkers from the process + new_coords, e_new, w_new = read(proc, len_walk) + walkers += new_coords + + # Print energy + ew = e_new * w_new + EnergyWeight += ew + Energy2Weight += e_new * ew + Weight += w_new + NSamples += 1. + print (len(walkers)) + err = print_energy(EnergyWeight, Energy2Weight, Weight, NSamples) + + if err < 1.e-3: + do_loop = False + + terminate(proc) + return + + + + +if __name__ == "__main__": + main() diff --git a/src/MAIN/vmc_test.irp.f b/src/MAIN/vmc_test.irp.f index 4876f6d..fb54a9b 100644 --- a/src/MAIN/vmc_test.irp.f +++ b/src/MAIN/vmc_test.irp.f @@ -1,4 +1,4 @@ -program vmc_test +program vmc_test real :: t1,t2 print *, 'Ndet=',det_num print *, 'Ndet alpha beta =',det_alpha_num, det_beta_num diff --git a/src/Makefile b/src/Makefile index 5964b67..68e3a23 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,3 +1,4 @@ +IRPF90=irpf90 IRPF90+= $(IRPF90_FLAGS) include irpf90.make export