From 6ee4518bb17a3d3dd3aba0b793321a6a474ef374 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 20 Mar 2019 20:46:47 +0200 Subject: [PATCH] Complex Bessel functions using amos. Not working yet, probably because of unset machine constants in i2mach.f, d1mach.f. Former-commit-id: 94a15b3bbbcb26b21533866c984076aeefcecb5e --- CMakeLists.txt | 5 +++ README.rst | 4 ++ amos/CMakeLists.txt | 8 ++-- amos/amos.h | 42 +++++++++++++++++++ qpms/CMakeLists.txt | 3 +- qpms/bessel.c | 97 ++++++++++++++++++++++++++++++++++++++++---- qpms/error.c | 12 ++++++ qpms/qpms_error.h | 10 +++++ qpms/qpms_specfunc.h | 4 +- qpms/translations.c | 3 ++ setup.py | 57 ++++++++++++++++++++++++-- 11 files changed, 227 insertions(+), 18 deletions(-) create mode 100644 amos/amos.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 502b77b..22cc8e1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,6 +12,11 @@ macro(use_c99) endif () endmacro(use_c99) +# We need the amos libraries to be compiled with -fPIC +# but at the same time, we do not want to make a .so file, +# so distutils package it with the rest of qpms c/cython lib. +set(CMAKE_POSITION_INDEPENDENT_CODE ON) + set (QPMS_VERSION_MAJOR 0) #set (QPMS_VERSION_MINOR 3) cmake_add_fortran_subdirectory (amos diff --git a/README.rst b/README.rst index 34a7e37..4ce7002 100644 --- a/README.rst +++ b/README.rst @@ -13,8 +13,12 @@ just run ``apt-get install libgsl-dev`` under root. Alternatively, you can `get the source `_ get the source and compile it yourself. +You also need a fresh enough version of ``cmake``. + After GSL is installed, you can install qpms to your local python library using:: + cmake . + make amos python3 setup.py install --user If GSL is not installed the standard library path on your system, you might diff --git a/amos/CMakeLists.txt b/amos/CMakeLists.txt index d3717dc..73d9107 100644 --- a/amos/CMakeLists.txt +++ b/amos/CMakeLists.txt @@ -1,10 +1,10 @@ enable_language (Fortran) include(FortranCInterface) -FortranCInterface_HEADER(amos.h +FortranCInterface_HEADER(amos_mangling.h MACRO_NAMESPACE "AMOS_" SYMBOL_NAMESPACE "amos_" - SYMBOLS zbesj zbesy zbesi zbesk zsqrt + SYMBOLS zbesj zbesy zbesh zbesi zbesk zsqrt ) add_library(amos @@ -13,4 +13,6 @@ add_library(amos i1mach.f zacon.f zbesj.f zbuni.f zkscl.f zrati.f zs1s2.f zuni1.f zuoik.f xerror.f zairy.f zbesy.f zbunk.f zlog.f zseri.f zuchk.f zuni2.f zwrsk.f ) - + +target_include_directories (amos PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) + diff --git a/amos/amos.h b/amos/amos.h new file mode 100644 index 0000000..854a28f --- /dev/null +++ b/amos/amos.h @@ -0,0 +1,42 @@ +#ifndef AMOS_H +#define AMOS_H +#include "amos_mangling.h" + +#define INTEGER_t int +#define DOUBLE_PRECISION_t double + +void amos_zbesj(const DOUBLE_PRECISION_t *zr, + const DOUBLE_PRECISION_t *zi, + const DOUBLE_PRECISION_t *fnu, + const INTEGER_t *kode, + const INTEGER_t *n, + DOUBLE_PRECISION_t *cyr, + DOUBLE_PRECISION_t *cyi, + INTEGER_t *nz, + INTEGER_t *ierr); + +void amos_zbesy(const DOUBLE_PRECISION_t *zr, + const DOUBLE_PRECISION_t *zi, + const DOUBLE_PRECISION_t *fnu, + const INTEGER_t *kode, + const INTEGER_t *n, + DOUBLE_PRECISION_t *cyr, + DOUBLE_PRECISION_t *cyi, + INTEGER_t *nz, + DOUBLE_PRECISION_t *cwrkr, + DOUBLE_PRECISION_t *cwrki, + INTEGER_t *ierr); + +void amos_zbesh(const DOUBLE_PRECISION_t *zr, + const DOUBLE_PRECISION_t *zi, + const DOUBLE_PRECISION_t *fnu, + const INTEGER_t *kode, + const INTEGER_t *m, + const INTEGER_t *n, + DOUBLE_PRECISION_t *cyr, + DOUBLE_PRECISION_t *cyi, + INTEGER_t *nz, + INTEGER_t *ierr); + + +#endif diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index b4c53c2..f2ca92d 100644 --- a/qpms/CMakeLists.txt +++ b/qpms/CMakeLists.txt @@ -8,7 +8,8 @@ set (DIRS ${GSL_INCLUDE_DIRS} ${GSLCBLAS_INCLUDE_DIRS}) include_directories(${DIRS}) add_library (qpms translations.c tmatrices.c vecprint.c vswf.c wigner.c - lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c) + lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c + bessel.c) use_c99() set(LIBS ${LIBS} ${GSL_LIBRARIES} ${GSLCBLAS_LIBRARIES}) diff --git a/qpms/bessel.c b/qpms/bessel.c index 825e30d..84a9336 100644 --- a/qpms/bessel.c +++ b/qpms/bessel.c @@ -7,6 +7,7 @@ #include #include #include "qpms_error.h" +#include #ifndef M_LN2 #define M_LN2 0.69314718055994530942 /* log_e 2 */ @@ -16,6 +17,36 @@ static inline complex double ipow(int x) { return cpow(I,x); } +#if 0 +// Inspired by scipy/special/_spherical_bessel.pxd +static inline complex double spherical_jn(qpms_l_t l, complex double z) { + if (isnan(creal(z)) || isnan(cimag(z))) return NAN+I*NAN; + if (l < 0) QPMS_WTF; + if (fpclassify(creal(z)) == FP_INFINITE) + if(0 == cimag(z)) + return 0; + else + return INFINITY + I * INFINITY; + if (z == 0) + if (l == 0) return 1; + else return 0; + return csqrt(M_PI_2/z) * zbesj(l + .5, z); +} + +static inline complex double spherical_yn(qpms_l_t l, complex double z) { + if (isnan(creal(z)) || isnan(cimag(z))) return NAN+I*NAN; + if (l < 0) QPMS_WTF; + if (fpclassify(creal(z)) == FP_INFINITE) + if(0 == cimag(z)) + return 0; + else + return INFINITY + I * INFINITY; + if (z == 0) + return NAN; + return csqrt(M_PI_2/z) * zbesy(l + .5, z); +} +#endif + // There is a big issue with gsl's precision of spherical bessel function; these have to be implemented differently qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array) { int retval; @@ -50,12 +81,60 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double assert(0); } -qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex double x, complex double *result_array) { - int retval = qpms_sph_bessel_realx_fill(typ, lmax, creal(x), result_array); +// TODO DOC +qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex double x, complex double *res) { if(!cimag(x)) - return retval; - else - QPMS_NOT_IMPLEMENTED("complex argument of Bessel function."); + return qpms_sph_bessel_realx_fill(typ, lmax, creal(x), res); + else if (isnan(creal(x)) || isnan(cimag(x))) + for(qpms_l_t l = 0; l <= lmax; ++l) res[l] = NAN + I*NAN; + else if (lmax < 0) QPMS_WTF; + else if (fpclassify(creal(x)) == FP_INFINITE) + for(qpms_l_t l = 0; l <= lmax; ++l) res[l] = INFINITY + I * INFINITY; + else { + const DOUBLE_PRECISION_t zr = creal(x), zi = cimag(x), fnu = 0.5; + const INTEGER_t n = lmax + 1, kode = 1 /* No exponential scaling */; + DOUBLE_PRECISION_t cyr[n], cyi[n]; + INTEGER_t ierr, nz; + unsigned int kindchar; // Only for error output + const complex double prefac = csqrt(M_PI_2/x); + switch(typ) { + case QPMS_BESSEL_REGULAR: + kindchar = 'j'; + amos_zbesj(&zr, &zi, &fnu, &kode, &n, cyr, cyi, &nz, &ierr); + break; + case QPMS_BESSEL_SINGULAR: + kindchar = 'y'; + { + DOUBLE_PRECISION_t cwrkr[lmax + 1], cwrki[lmax + 1]; + amos_zbesy(&zr, &zi, &fnu, &kode, &n, cyr, cyi, &nz, cwrkr, cwrki, &ierr); + } + break; + case QPMS_HANKEL_PLUS: + case QPMS_HANKEL_MINUS: + kindchar = 'h'; + { + const INTEGER_t m = (typ == QPMS_HANKEL_PLUS) ? 1 : 2; + amos_zbesh(&zr, &zi, &fnu, &kode, &m, &n, cyr, cyi, &nz, &ierr); + } + break; + default: + QPMS_WTF; + } + // TODO check for underflows? (nz != 0) + if (ierr == 0 || ierr == 3) { + for (qpms_l_t l = 0; l <= lmax; ++l) + res[l] = prefac * (cyr[l] + I * cyi[l]); + if (ierr == 3) + QPMS_WARN("Amos's zbes%c computation done but losses of significance " + "by argument reduction produce less than half of machine accuracy.", + kindchar); + return QPMS_SUCCESS; //TODO maybe something else if ierr == 3 + } + else + QPMS_PR_ERROR("Amos's zbes%c failed with ierr == %d.", + kindchar, (int) ierr); + } + return QPMS_SUCCESS; } static inline ptrdiff_t akn_index(qpms_l_t n, qpms_l_t k) { @@ -84,10 +163,10 @@ static inline qpms_errno_t qpms_sbessel_calculator_ensure_lMax(qpms_sbessel_calc } } -complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, double x) { +complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, complex double x) { if(QPMS_SUCCESS != qpms_sbessel_calculator_ensure_lMax(c, n)) abort(); - complex double z = I/x; // FIXME this should be imaginary double, but gcc is broken? + complex double z = I/x; complex double result = 0; for (qpms_l_t k = n; k >= 0; --k) // can we use fma for complex? @@ -98,14 +177,14 @@ complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, do } qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t * const c, - const qpms_l_t lMax, const double x, complex double * const target) { + const qpms_l_t lMax, const complex double x, complex double * const target) { if(QPMS_SUCCESS != qpms_sbessel_calculator_ensure_lMax(c, lMax)) abort(); memset(target, 0, sizeof(complex double) * lMax); complex double kahancomp[lMax]; memset(kahancomp, 0, sizeof(complex double) * lMax); for(qpms_l_t k = 0; k <= lMax; ++k){ - double xp = pow(x, -k-1); + double xp = cpow(x, -k-1); for(qpms_l_t l = k; l <= lMax; ++l) ckahanadd(target + l, kahancomp + l, c->akn[akn_index(l,k)] * xp * ipow(k-l-1)); } diff --git a/qpms/error.c b/qpms/error.c index 305cb11..2a8ec18 100644 --- a/qpms/error.c +++ b/qpms/error.c @@ -44,6 +44,18 @@ void qpms_pr_error_at_flf(const char *filename, unsigned int linenum, abort(); } +void qpms_warn_at_flf(const char *filename, unsigned int linenum, + const char *func, + const char *fmt, ...) { + fprintf(stderr, "%s:%u, %s: ", filename, linenum, func); + va_list ap; + va_start(ap, fmt); + vfprintf(stderr, fmt, ap); + va_end(ap); + fputc('\n', stderr); + fflush(stderr); +} + void qpms_pr_debug_at_flf(const char *filename, unsigned int linenum, const char *func, const char *fmt, ...) { diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h index 325a28a..c6591fc 100644 --- a/qpms/qpms_error.h +++ b/qpms/qpms_error.h @@ -13,18 +13,28 @@ void qpms_pr_error_at_flf(const char *filename, unsigned int linenum, const char *func, const char *fmt, ...); +/// Print a warning message to stderr and flush the buffer. Don't call this directly, use QPMS_WARN(). +void qpms_warn_at_flf(const char *filename, unsigned int linenum, + const char *func, + const char *fmt, ...); + void qpms_pr_debug_at_flf(const char *filename, unsigned int linenum, const char *func, const char *fmt, ...); //void qpms_error_at_line(const char *filename, unsigned int linenum, // const char *fmt, ...); + +#define QPMS_WARN(msg, ...) qpms_warn_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__) + #define QPMS_CRASHING_MALLOC(pointer, size) {(pointer) = malloc(size); if(!pointer && (size)) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Allocation of %zd bytes for " #pointer " failed.", (size_t) (size));} #define QPMS_CRASHING_REALLOC(pointer, size) {(pointer) = realloc(pointer, size); if(!pointer && (size)) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Rellocation of %zd bytes for " #pointer " failed.", (size_t) (size));} #define QPMS_WTF qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error.") +#define QPMS_PR_ERROR(msg, ...) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__) + #define QPMS_ENSURE_SUCCESS(x) {if(x) QPMS_WTF;} #define QPMS_ENSURE(x, msg, ...) {if(!(x)) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); } diff --git a/qpms/qpms_specfunc.h b/qpms/qpms_specfunc.h index f95d819..fe2fd55 100644 --- a/qpms/qpms_specfunc.h +++ b/qpms/qpms_specfunc.h @@ -24,9 +24,9 @@ void qpms_sbessel_calculator_pfree(qpms_sbessel_calculator_t *c); qpms_errno_t qpms_sbessel_calc_fill(qpms_sbessel_calculator_t *c, qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array); -complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, double x); +complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, complex double x); qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t *c, qpms_l_t lmax, - double x, complex double *result_array); + complex double x, complex double *result_array); /****************************************************************************** diff --git a/qpms/translations.c b/qpms/translations.c index 4dbdbea..ddbb618 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -1,5 +1,6 @@ #include #include "qpms_types.h" +#include "qpms_specfunc.h" #include "gaunt.h" #include "translations.h" #include "indexing.h" // TODO replace size_t and int with own index types here @@ -89,6 +90,7 @@ static inline int gauntB_Q_max(int M, int n, int mu, int nu) { return imin(n, imin(nu, (n+nu+1-abs(M+mu))/2)); } +#if 0 //Apparently, this is duplicit to that in bessel.c (which also support complex x) int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) { int retval; double tmparr[lmax+1]; @@ -121,6 +123,7 @@ int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double * } assert(0); } +#endif static inline double qpms_trans_normlogfac(qpms_normalisation_t norm, int m, int n, int mu, int nu) { diff --git a/setup.py b/setup.py index cf067dd..486d0c3 100755 --- a/setup.py +++ b/setup.py @@ -17,6 +17,50 @@ print("You might want to add additional library path to LD_LIBRARY_PATH (especia if("LD_LIBRARY_PATH" in os.environ): print(os.environ['LD_LIBRARY_PATH'].split(':')) + +amos_sources = [ + 'amos/d1mach.f', + 'amos/dgamln.f', + 'amos/i1mach.f', + 'amos/xerror.f', + 'amos/zabs.f', + 'amos/zacai.f', + 'amos/zacon.f', + 'amos/zairy.f', + 'amos/zasyi.f', + 'amos/zbesh.f', + 'amos/zbesi.f', + 'amos/zbesj.f', + 'amos/zbesk.f', + 'amos/zbesy.f', + 'amos/zbinu.f', + 'amos/zbiry.f', + 'amos/zbknu.f', + 'amos/zbuni.f', + 'amos/zbunk.f', + 'amos/zdiv.f', + 'amos/zexp.f', + 'amos/zkscl.f', + 'amos/zlog.f', + 'amos/zmlri.f', + 'amos/zmlt.f', + 'amos/zrati.f', + 'amos/zs1s2.f', + 'amos/zseri.f', + 'amos/zshch.f', + 'amos/zsqrt.f', + 'amos/zuchk.f', + 'amos/zunhj.f', + 'amos/zuni1.f', + 'amos/zuni2.f', + 'amos/zunik.f', + 'amos/zunk1.f', + 'amos/zunk2.f', + 'amos/zuoik.f', + 'amos/zwrsk.f', + ] + + qpms_c = Extension('qpms_c', sources = ['qpms/qpms_c.pyx', #'qpms/hexpoints_c.pyx', 'qpms/gaunt.c',#'qpms/gaunt.h','qpms/vectors.h','qpms/translations.h', @@ -28,7 +72,9 @@ qpms_c = Extension('qpms_c', 'qpms/vswf.c', # FIXME many things from vswf.c are not required by this module, but they have many dependencies (following in this list); maybe I want to move all the "basespec stuff" 'qpms/legendre.c', 'qpms/tmatrices.c', - 'qpms/error.c'], + 'qpms/error.c', + 'qpms/bessel.c', + ], extra_compile_args=['-std=c99','-ggdb', '-O3', '-DQPMS_COMPILE_PYTHON_EXTENSIONS', # this is required #'-DQPMS_USE_OMP', @@ -37,17 +83,22 @@ qpms_c = Extension('qpms_c', ], libraries=['gsl', 'lapacke', 'blas', 'gslcblas', 'pthread', #'omp' # TODO resolve the problem with openblas (missing gotoblas symbol) and preferable use other blas library + #('amos', dict(sources=amos_sources) ), ], - runtime_library_dirs=os.environ['LD_LIBRARY_PATH'].split(':') if 'LD_LIBRARY_PATH' in os.environ else [] + include_dirs=['amos'], + extra_link_args=['amos/libamos.a'], + runtime_library_dirs=os.environ['LD_LIBRARY_PATH'].split(':') if 'LD_LIBRARY_PATH' in os.environ else [], + #extra_objects = ['amos/libamos.a'], # FIXME apparently, I would like to eliminate the need to cmake/make first ) setup(name='qpms', version = "0.2.995", packages=['qpms'], +# libraries = [('amos', {'sources': amos_sources} )], setup_requires=['cython>=0.28',], install_requires=['cython>=0.28','quaternion','spherical_functions','scipy>=0.18.0', 'sympy>=1.2'], dependency_links=['https://github.com/moble/quaternion/archive/v2.0.tar.gz','https://github.com/moble/spherical_functions/archive/master.zip'], - ext_modules=cythonize([qpms_c], include_path=['qpms'], gdb_debug=True), + ext_modules=cythonize([qpms_c], include_path=['qpms', 'amos'], gdb_debug=True), cmdclass = {'build_ext': build_ext}, )