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},
)