[WIP] C API

CMakeLists.txt

@@ -4,27 +4,41 @@ cmake_minimum_required( VERSION 3.17 ) # Require CMake 3.17+
 project( IntegratorXX VERSION 1.2.0 LANGUAGES CXX )
 
 
-add_library( integratorxx INTERFACE )
-add_library( IntegratorXX::IntegratorXX ALIAS integratorxx )
 
 #option( INTEGRATORXX_ENABLE_TESTS "Enable Tets" ON )
+option( INTEGRATORXX_ENABLE_C_API "Enable C Bindings" ON )
 
 
+if(INTEGRATORXX_ENABLE_C_API)
+  set( INTEGRATORXX_HEADER_ONLY FALSE )
+  enable_language(C)
+  add_subdirectory(src/c_api)
+else()
+  # Header only if no C-API
+  set( INTEGRATORXX_HEADER_ONLY TRUE )
+  add_library( integratorxx INTERFACE )
+endif()
+add_library( IntegratorXX::IntegratorXX ALIAS integratorxx )
 
 
 
 # Target features
-target_compile_features( integratorxx INTERFACE cxx_std_17 )
+if(INTEGRATORXX_HEADER_ONLY)
+  set(INTEGRATORXX_PROPERTY_TYPE INTERFACE)
+else()
+  set(INTEGRATORXX_PROPERTY_TYPE PUBLIC)
+endif()
+target_compile_features( integratorxx ${INTEGRATORXX_PROPERTY_TYPE} cxx_std_17 )
 target_include_directories( integratorxx
-  INTERFACE
+  ${INTEGRATORXX_PROPERTY_TYPE}
     $<INSTALL_INTERFACE:include>
     $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include>
 )
 
 include(CheckCXXCompilerFlag)
 check_cxx_compiler_flag("-Wno-missing-braces" INTEGRATORXX_HAS_NO_MISSING_BRACES )
 if( INTEGRATORXX_HAS_NO_MISSING_BRACES )
-	target_compile_options( integratorxx INTERFACE $<$<COMPILE_LANGUAGE:CXX>: -Wno-missing-braces> )
+	target_compile_options( integratorxx ${INTEGRATORXX_PROPERTY_TYPE} $<$<COMPILE_LANGUAGE:CXX>: -Wno-missing-braces> )
 endif()
 
 

include/integratorxx/c_api.h

@@ -0,0 +1,158 @@
+#pragma once
+
+/*** Error codes ***/
+#define INTXX_SUCCESS       0
+#define INTXX_INVALID_QUAD -1
+#define INTXX_NULL_QUADPTR -2
+#define INTXX_NULL_INFOPTR -3
+#define INTXX_INVALID_OPT  -4
+#define INTXX_INVALID_ARG  -5
+#define INTXX_INVALID_OUT  -6
+
+/*** Quadrature Classes ***/
+#define INTXX_PRM_QUAD 1
+#define INTXX_RAD_QUAD 2
+#define INTXX_ANG_QUAD 3
+#define INDXX_SPH_QUAD 4
+
+/*** Primitive Quadratures ***/
+#define INTXX_PRMQ_MASK           0x0000FF
+#define INTXX_PRMQ_UNIFORM        0x000001
+#define INTXX_PRMQ_GAUSSLEG       0x000002
+#define INTXX_PRMQ_GAUSSCHEB_1    0x000003
+#define INTXX_PRMQ_GAUSSCHEB_2    0x000004
+#define INTXX_PRMQ_GAUSSCHEB_2MOD 0x000005
+#define INTXX_PRMQ_GAUSSCHEB_3    0x000006
+#define INTXX_PRMQ_GAUSSLOB       0x000007
+
+/*** Radial Quadratures ***/
+#define INTXX_RADQ_MASK  0x00FF00
+#define INTXX_RADQ_BECKE 0x000100 // Becke
+#define INTXX_RADQ_MHL   0x000200 // Murray-Handy-Laming
+#define INTXX_RADQ_TA    0x000300 // Treuter-Ahlrichs
+#define INTXX_RADQ_MK    0x000400 // Mura-Knowles
+
+/*** Angular (S2) Quadratures ***/
+#define INTXX_ANGQ_MASK 0xFF0000
+#define INTXX_ANGQ_LEB  0x010000 // Lebedev-Laikov
+#define INTXX_ANGQ_DEL  0x020000 // Delley
+#define INTXX_ANGQ_AB   0x030000 // Ahrens-Beylkin
+#define INTXX_ANGQ_WOM  0x040000 // Womersley
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// Forward defs of types
+struct intxx_quad_type;
+
+typedef struct {
+  int n; ///< Number of parameters
+
+  const char** names;        ///< Names of the params
+  const char** descriptions; ///< Long descriptions of params
+
+  int (*set_name)(struct intxx_quad_type* p, const char* name, double v);
+  ///< Set parameter function
+  int (*get_name)(struct intxx_quad_type* p, const char* name, double* v);
+    ///< Get parameter function
+  int (*generate)(struct intxx_quad_type* p);
+  int (*destroy) (struct intxx_quad_type* p);
+} intxx_quad_params_type;
+
+typedef struct {
+  int number; ///< Quadrature identifier
+  int kind;   ///< Type of quadrature (PRM, RAD, ANG, SPH)
+  int dim;    ///< Dimension of the grid (1D or 3D)
+
+  const char* name; ///< Name of the functional, e.g. "Becke"
+  // TODO References
+
+  intxx_quad_params_type ext_params; ///< External params
+
+  int (*generate)(struct intxx_quad_type* p);
+  int (*destroy) (struct intxx_quad_type* p);
+  int (*set_npts)(struct intxx_quad_type* p, int);
+} intxx_quad_info_type;
+
+struct intxx_quad_type {
+  int npoints;
+  const intxx_quad_info_type* info;
+  void* _state_quad; ///< Internal state, NOT FOR EXTERNAL USE
+  void* _state_parm; ///< Internal state, NOT FOR EXTERNAL USE
+};
+
+typedef struct intxx_quad_type intxx_quad_type;
+
+
+
+/**
+ *  @brief Initialize a specified quadrature
+ *
+ *  See error code returns for how to interpret failures.
+ *
+ *  @param[out] p Quadrature instance of the specified type
+ *  @param[in]  quad Type of the specified quadrature
+ *
+ *  @returns INTXX_SUCCESS:      no errors were encountered
+ *           INTXX_NULL_QUADPTR: `p == NULL`
+ *           INTXX_INVALID_QUAD: `quad` is invalid
+ */
+int  intxx_quad_init(intxx_quad_type* p, int quad);
+
+/// Frees a quadrature instance. No throw gurantee
+void intxx_quad_end (intxx_quad_type* p);
+
+/**
+ *  @brief Set the number of quadrature points
+ *
+ *  Sets the number of nodes for the passed quadrature 
+ *  instance. Only sensible for quadratures in which
+ *  the size is a free parameter (most).
+ *
+ *  See error code returns for how to interpret failures.
+ *
+ *  @param[out] p The quadrature for which to set npts.
+ *  @param[in]  npts The number of points
+ *
+ *  @returns INTXX_SUCCESS:      no errors were encountered
+ *           INTXX_NULL_QUADPTR: `p == NULL`
+ *           INTXX_NULL_INFOPTR: `p->info == NULL` (uninit)
+ *           INTXX_INVALID_ARG:  invalid npts (e.g. npts < 0)
+ */
+int intxx_quad_set_npts(intxx_quad_type* p, int npts);
+
+/**
+ *  @brief Retrieve the number of points for a quadrature
+ *
+ *  Retreieves the number of grid points for a specified 
+ *  quadrature. In necessicary, this will be initialized 
+ *  from the default parameters, otherwise it will be
+ *  read from p->npoints.
+ *
+ *  See error code returns for how to interpret failures.
+ *
+ *  @param[in]   p    The quadrature for which to get npts.
+ *  @param[out]  npts The number of points
+ *
+ *  @returns INTXX_SUCCESS:      no errors were encountered
+ *           INTXX_NULL_QUADPTR: `p == NULL`
+ *           INTXX_NULL_INFOPTR: `p->info == NULL` (uninit)
+ *           INTXX_INVALID_ARG:  `npts == NULL` 
+ *           INTXX_INVALID_OUT:  `npts` is not valid (npts < 0)
+ */ 
+int intxx_quad_get_npts(intxx_quad_type* p, int* npts);
+
+
+int intxx_generate_quad(intxx_quad_type* p);
+int intxx_destroy_quad(intxx_quad_type* p);
+
+int intxx_get_ext_param_name(intxx_quad_type* p, const char* name, double* v);
+int intxx_set_ext_param_name(intxx_quad_type* p, const char* name, double v);
+
+int intxx_get_rad_scal(intxx_quad_type* p, double *R);
+int intxx_set_rad_scal(intxx_quad_type* p, double R);
+
+#ifdef __cplusplus
+}
+#endif

include/integratorxx/quadratures/radial/becke.hpp

@@ -47,6 +47,9 @@ class BeckeRadialTraits {
   inline auto radial_jacobian(PointType x) const noexcept {
     return R_ * 2.0 / std::pow(1.0 - x, 2);
   }
+
+  /// Return radial scaling factor
+  auto R() const { return R_; }
 };
 
 /**

include/integratorxx/quadratures/radial/mhl.hpp

@@ -47,6 +47,8 @@ class MurrayHandyLamingRadialTraits {
     return R_ * M * std::pow(x, M-1) / std::pow(1.0 - x, M+1);
   }
 
+  /// Return radial scaling factor
+  auto R() const { return R_; }
 };
 
 

include/integratorxx/quadratures/radial/muraknowles.hpp

@@ -155,6 +155,8 @@ class MuraKnowlesRadialTraits {
     return R_ * 3.0 * x2 / (1.0 - x2 * x);
   }
 
+  /// Return radial scaling factor
+  auto R() const { return R_; }
 }; 
 
 template <typename PointType, typename WeightType>

include/integratorxx/quadratures/radial/radial_transform.hpp

@@ -12,6 +12,9 @@ struct RadialTransformQuadrature :
 
 public:
 
+  using quad_type   = BaseQuad;
+  using traits_type = RadialTraits;
+
   using point_type       = typename base_type::point_type;
   using weight_type      = typename base_type::weight_type;
   using point_container  = typename base_type::point_container;

include/integratorxx/quadratures/radial/treutlerahlrichs.hpp

@@ -64,6 +64,8 @@ class TreutlerAhlrichsRadialTraits {
     return R_ * pow_term / ln_2 * ( alpha_ * log_term / (a+x) + (1./(1.-x)) );
   }
 
+  /// Return radial scaling factor
+  auto R() const { return R_; }
 };
 
 

src/c_api/CMakeLists.txt

@@ -0,0 +1 @@
+add_library( integratorxx c_api.c c_primitive.cxx c_radial.cxx c_misc.cxx c_angular.cxx )

src/c_api/c_angular.cxx

@@ -0,0 +1,117 @@
+#include <integratorxx/quadratures/s2.hpp>
+
+#include "c_internal.h"
+#include "c_api_util.hpp"
+
+
+extern "C" {
+
+using leb_quad_type    = IntegratorXX::LebedevLaikov<double>;
+using delley_quad_type = IntegratorXX::Delley<double>;
+using wom_quad_type    = IntegratorXX::Womersley<double>;
+using ab_quad_type     = IntegratorXX::AhrensBeylkin<double>;
+
+
+/*******************************************************/
+/****** Forward declare C-API for ANG Quadratures ******/
+/*******************************************************/
+
+#define FWD_ANG(cname) \
+int intxx_get_##cname##_info(intxx_quad_info_type* p); \
+int intxx_generate_##cname(intxx_quad_type* p); \
+int intxx_destroy_##cname(intxx_quad_type* p); \
+
+FWD_ANG(leb)
+FWD_ANG(delley)
+FWD_ANG(wom)
+FWD_ANG(ab)
+
+#undef FWD_ANG
+
+/*******************************************************/
+/****** Runtime Generator for ANG Quadrature Info ******/
+/*******************************************************/
+
+int intxx_get_angq_info(intxx_quad_info_type* p, int quad) {
+  if(p == NULL) return INTXX_NULL_INFOPTR;
+
+  quad = quad & INTXX_ANGQ_MASK;
+  p->number = quad;
+  p->kind   = INTXX_ANG_QUAD;
+  p->dim    = 3;
+  switch(quad) {
+    case INTXX_ANGQ_LEB:
+      p->name = "LEBEDEV_LAIKOV";
+      return intxx_get_leb_info(p);
+    case INTXX_ANGQ_DEL:
+      p->name = "DELLEY";
+      return intxx_get_delley_info(p);
+    case INTXX_ANGQ_AB:
+      p->name = "AHRENS_BEYLKIN";
+      return intxx_get_ab_info(p);
+    case INTXX_ANGQ_WOM:
+      p->name = "WOMERSLEY";
+      return intxx_get_wom_info(p);
+    default:
+      return INTXX_INVALID_QUAD;
+  }
+  return INTXX_SUCCESS;
+}
+
+/**********************************************/
+/****** NPTS Setters for ANG Quadratures ******/
+/**********************************************/
+
+#define INTXX_ANG_SET_NPTS(cname, nsp) \
+int intxx_##cname##_set_npts(intxx_quad_type* p, int npts) {\
+  using namespace IntegratorXX::detail::nsp;\
+  if(algebraic_order_by_npts(npts) > 0) {\
+    p->npoints = npts;\
+    return INTXX_SUCCESS;\
+  } else {\
+    p->npoints = -1;\
+    return INTXX_INVALID_ARG;\
+  }\
+}
+
+INTXX_ANG_SET_NPTS(leb, lebedev)
+INTXX_ANG_SET_NPTS(delley, delley)
+INTXX_ANG_SET_NPTS(ab, ahrensbeylkin)
+INTXX_ANG_SET_NPTS(wom, womersley)
+
+/*************************************************/
+/****** Info generation for ANG Quadratures ******/
+/*************************************************/
+
+#define INTXX_NOPARAM_GET_INFO_IMPL(cname) \
+int intxx_get_##cname##_info(intxx_quad_info_type* p) { \
+  return intxx_noparam_info(p, &intxx_generate_##cname, \
+    &intxx_destroy_##cname, intxx_##cname##_set_npts);  \
+}
+
+INTXX_NOPARAM_GET_INFO_IMPL(leb);
+INTXX_NOPARAM_GET_INFO_IMPL(delley);
+INTXX_NOPARAM_GET_INFO_IMPL(ab);
+INTXX_NOPARAM_GET_INFO_IMPL(wom);
+
+#undef INTXX_NOPARAM_GET_INFO_IMPL
+
+/*************************************************************/
+/****** Allocate/Deallocate ANG Quadrature States ******/
+/*************************************************************/
+
+#define INTXX_GD_BASIC_IMPL(cname, qname, ...)   \
+int intxx_generate_##cname(intxx_quad_type* p) { \
+  return intxx_generate_noparam_impl<qname>(p, ##__VA_ARGS__);  \
+}                                                \
+int intxx_destroy_##cname(intxx_quad_type* p) {  \
+  return intxx_destroy_impl<qname>(p);   \
+}
+
+INTXX_GD_BASIC_IMPL(leb, leb_quad_type);
+INTXX_GD_BASIC_IMPL(delley, delley_quad_type);
+INTXX_GD_BASIC_IMPL(ab, ab_quad_type);
+INTXX_GD_BASIC_IMPL(wom, wom_quad_type);
+
+#undef INTXX_GD_BASIC_IMPL
+}