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

@@ -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

README.rst

@@ -13,8 +13,12 @@ just run ``apt-get install libgsl-dev`` under root. Alternatively,
 you can `get the source 
 <https://www.gnu.org/software/gsl/>`_ 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 

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

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

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

qpms/bessel.c

@@ -7,6 +7,7 @@
 #include <gsl/gsl_sf_bessel.h>
 #include <complex.h>
 #include "qpms_error.h"
+#include <amos.h>
 
 #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) {
+// TODO DOC
-  int retval = qpms_sph_bessel_realx_fill(typ, lmax, creal(x), result_array);
+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;
+    return qpms_sph_bessel_realx_fill(typ, lmax, creal(x), res);
-  else
+  else if (isnan(creal(x)) || isnan(cimag(x))) 
-    QPMS_NOT_IMPLEMENTED("complex argument of Bessel function.");
+    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));
   }

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, ...) {

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__); }

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);
 
 
 /******************************************************************************

qpms/translations.c

@@ -1,5 +1,6 @@
 #include <math.h>
 #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) {

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