From 625ccb1d8decf294b67b735be6a63cf26782759c Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 23 Sep 2023 16:49:42 -0700 Subject: [PATCH 01/22] Start C API, able to generate metadata and set/get npts from primitive quadratures --- CMakeLists.txt | 24 ++++++++-- include/integratorxx/c_api.h | 89 ++++++++++++++++++++++++++++++++++++ src/c_api/CMakeLists.txt | 1 + src/c_api/c_api.c | 87 +++++++++++++++++++++++++++++++++++ src/c_api/c_internal.h | 16 +++++++ src/c_api/c_radial.cxx | 84 ++++++++++++++++++++++++++++++++++ test/CMakeLists.txt | 3 ++ test/c_api.cxx | 79 ++++++++++++++++++++++++++++++++ 8 files changed, 378 insertions(+), 5 deletions(-) create mode 100644 include/integratorxx/c_api.h create mode 100644 src/c_api/CMakeLists.txt create mode 100644 src/c_api/c_api.c create mode 100644 src/c_api/c_internal.h create mode 100644 src/c_api/c_radial.cxx create mode 100644 test/c_api.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index 45337ce..1188a67 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,19 +4,33 @@ 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} $ $ ) @@ -24,7 +38,7 @@ target_include_directories( integratorxx 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 $<$: -Wno-missing-braces> ) + target_compile_options( integratorxx ${INTEGRATORXX_PROPERTY_TYPE} $<$: -Wno-missing-braces> ) endif() diff --git a/include/integratorxx/c_api.h b/include/integratorxx/c_api.h new file mode 100644 index 0000000..281587f --- /dev/null +++ b/include/integratorxx/c_api.h @@ -0,0 +1,89 @@ +#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 + +/*** 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 + const double *values; ///< Default values of params + + void (*set)(struct intxx_quad_type* p, const double** v); + ///< Set function +} intxx_quad_params_type; + +typedef struct { + int number; ///< Quadrature identifier + int kind; ///< Type of quadrature (PRM, RAD, ANG, SPH) + + const char* name; ///< Name of the functional, e.g. "Becke" + // TODO References + + intxx_quad_params_type ext_params; ///< External params + + void (*init)(struct intxx_quad_type* p); +} intxx_quad_info_type; + +struct intxx_quad_type { + int npoints; + const intxx_quad_info_type* info; +}; + +typedef struct intxx_quad_type intxx_quad_type; + + + + +int intxx_quad_init(intxx_quad_type* p, int quad); +void intxx_quad_end (intxx_quad_type* p); + +int intxx_quad_set_npts(intxx_quad_type* p, int npts); +int intxx_quad_get_npts(intxx_quad_type* p, int* npts); + +#ifdef __cplusplus +} +#endif diff --git a/src/c_api/CMakeLists.txt b/src/c_api/CMakeLists.txt new file mode 100644 index 0000000..12616f8 --- /dev/null +++ b/src/c_api/CMakeLists.txt @@ -0,0 +1 @@ +add_library( integratorxx c_api.c c_radial.cxx ) diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c new file mode 100644 index 0000000..548a0ec --- /dev/null +++ b/src/c_api/c_api.c @@ -0,0 +1,87 @@ +#include "c_internal.h" +#include +#include + + + +int intxx_quad_init(intxx_quad_type* p, int quad) { + // Local type defs + typedef intxx_quad_info_type info_type; + + // Sanity check + if(p == NULL) return INTXX_NULL_QUADPTR; + p->info = NULL; + p->npoints = -1; + + if(quad < 0) return INTXX_INVALID_QUAD; + + // Determine quadrature class + int is_prmq = quad & INTXX_PRMQ_MASK; + int is_radq = quad & INTXX_RADQ_MASK; + int is_angq = quad & INTXX_ANGQ_MASK; + int is_sphq = is_radq && is_angq; + + // Passed quadrature had to be something sane + if(!is_prmq && !is_radq && !is_angq) + return INTXX_INVALID_QUAD; + + // Primitive quadratures cannot be mixed with angular ones + if(is_prmq && is_radq) return INTXX_INVALID_QUAD; + if(is_prmq && is_angq) return INTXX_INVALID_QUAD; + + // Get info + info_type* finfo = (info_type*)malloc(sizeof(info_type)); + intxx_default_quad_info(finfo); // set (invalid) state + p->info = finfo; + + int error; + if(is_prmq) { + // Get primitive quadrature info + error = intxx_get_prmq_info(finfo, quad); + } else if(is_sphq) { + // TODO: Get spherical quadrature info + } else if(is_radq) { + // TODO: Get radial quadrature info + } else { + // Angular by exclusion + // TODO: Get angular quadrature info + } + + return error; +} + +void intxx_quad_end(intxx_quad_type* p) { + // Stateless - just bail + if(p == NULL || p->info == NULL) return; + + // Free info + free((void*)p->info); + p->info = NULL; +} + + + +int intxx_quad_set_npts(intxx_quad_type* p, int npts) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + // NPTS must be > 0 + if(npts <= 0) return INTXX_INVALID_ARG; + + // TODO: Handle the case when NPTS is derived from params + p->npoints = npts; + return INTXX_SUCCESS; +} + +int intxx_quad_get_npts(intxx_quad_type* p, int* npts) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + // Return memory must be valid + if(npts == NULL) return INTXX_INVALID_ARG; + + // TODO: Handle the case when NPTS is derived from params + *npts = p->npoints; + return INTXX_SUCCESS; +} + diff --git a/src/c_api/c_internal.h b/src/c_api/c_internal.h new file mode 100644 index 0000000..1ead231 --- /dev/null +++ b/src/c_api/c_internal.h @@ -0,0 +1,16 @@ +#pragma once +#include + +#ifdef __cplusplus +extern "C" { +#endif + +void intxx_default_quad_info(intxx_quad_info_type* p); +int intxx_get_radq_info(intxx_quad_info_type* p, int quad); +int intxx_get_angq_info(intxx_quad_info_type* p, int quad); +int intxx_get_prmq_info(intxx_quad_info_type* p, int quad); + +#ifdef __cplusplus +} +#endif + diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx new file mode 100644 index 0000000..6571643 --- /dev/null +++ b/src/c_api/c_radial.cxx @@ -0,0 +1,84 @@ +#include "c_internal.h" +#include + +extern "C" { + +int intxx_get_uniform_info(intxx_quad_info_type* p); +int intxx_get_gaussleg_info(intxx_quad_info_type* p); +int intxx_get_gausslob_info(intxx_quad_info_type* p); +int intxx_get_gausscheb1_info(intxx_quad_info_type* p); +int intxx_get_gausscheb2_info(intxx_quad_info_type* p); +int intxx_get_gausscheb2m_info(intxx_quad_info_type* p); +int intxx_get_gausscheb3_info(intxx_quad_info_type* p); + +void intxx_default_quad_info(intxx_quad_info_type* p) { + p->number = 0; // Invalid + p->kind = 0; // Invalid + p->name = "Default"; // Default state + p->init = NULL; +} + +int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + quad = quad & INTXX_PRMQ_MASK; + p->number = quad; + p->kind = INTXX_PRM_QUAD; + switch(quad) { + case INTXX_PRMQ_UNIFORM: + p->name = "UNIFORM"; + return intxx_get_uniform_info(p); + case INTXX_PRMQ_GAUSSLEG: + p->name = "GAUSS_LEGENDRE"; + return intxx_get_gaussleg_info(p); + case INTXX_PRMQ_GAUSSCHEB_1: + p->name = "GAUSS_CHEBYSHEV_1"; + return intxx_get_gausscheb1_info(p); + case INTXX_PRMQ_GAUSSCHEB_2: + p->name = "GAUSS_CHEBYSHEV_2"; + return intxx_get_gausscheb2_info(p); + case INTXX_PRMQ_GAUSSCHEB_2MOD: + p->name = "GAUSS_CHEBYSHEV_2MOD"; + return intxx_get_gausscheb2m_info(p); + case INTXX_PRMQ_GAUSSCHEB_3: + p->name = "GAUSS_CHEBYSHEV_3"; + return intxx_get_gausscheb3_info(p); + case INTXX_PRMQ_GAUSSLOB: + p->name = "GAUSS_LOBATTO"; + return intxx_get_gausslob_info(p); + default: + return INTXX_INVALID_QUAD; + } +} + +int intxx_noparam_info(intxx_quad_info_type* p) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + p->ext_params.n = 0; /// No External Parameters + + return INTXX_SUCCESS; +} + +int intxx_get_uniform_info(intxx_quad_info_type* p) { + return intxx_noparam_info(p); +} +int intxx_get_gausslob_info(intxx_quad_info_type* p) { + return intxx_noparam_info(p); +} +int intxx_get_gaussleg_info(intxx_quad_info_type* p) { + return intxx_noparam_info(p); +} +int intxx_get_gausscheb1_info(intxx_quad_info_type* p) { + return intxx_noparam_info(p); +} +int intxx_get_gausscheb2_info(intxx_quad_info_type* p) { + return intxx_noparam_info(p); +} +int intxx_get_gausscheb2m_info(intxx_quad_info_type* p) { + return intxx_noparam_info(p); +} +int intxx_get_gausscheb3_info(intxx_quad_info_type* p) { + return intxx_noparam_info(p); +} + +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 865e52e..7860108 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -24,6 +24,9 @@ target_link_libraries( integratorxx_common_ut PUBLIC Catch2::Catch2WithMain inte add_executable( quadrature_manipulation quadrature_manipulation.cxx ) target_link_libraries( quadrature_manipulation PUBLIC integratorxx_common_ut ) +add_executable( c_api c_api.cxx ) +target_link_libraries( c_api PUBLIC integratorxx_common_ut ) + add_executable( 1d_quadratures 1d_quadratures.cxx ) target_link_libraries( 1d_quadratures PUBLIC integratorxx_common_ut ) diff --git a/test/c_api.cxx b/test/c_api.cxx new file mode 100644 index 0000000..79c100a --- /dev/null +++ b/test/c_api.cxx @@ -0,0 +1,79 @@ +#include "catch2/catch_all.hpp" +#include +#include + +TEST_CASE("C API") { + + intxx_quad_type quad; + int error; + + SECTION("Invalid") { + error = intxx_quad_init(&quad, -1); + REQUIRE(error == INTXX_INVALID_QUAD); + REQUIRE(quad.info == NULL); + REQUIRE(quad.npoints == -1); + } + + SECTION("Primitive") { + + const char* name; + const int base_npts = 100; + + SECTION("Uniform") { + error = intxx_quad_init(&quad, INTXX_PRMQ_UNIFORM); + name = "UNIFORM"; + } + + SECTION("Gauss-Legendre") { + error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSLEG); + name = "GAUSS_LEGENDRE"; + } + + SECTION("Gauss-Lobatto") { + error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSLOB); + name = "GAUSS_LOBATTO"; + } + + SECTION("Gauss-Chebyshev 1") { + error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_1); + name = "GAUSS_CHEBYSHEV_1"; + } + + SECTION("Gauss-Chebyshev 2") { + error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_2); + name = "GAUSS_CHEBYSHEV_2"; + } + + SECTION("Gauss-Chebyshev 2MOD") { + error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_2MOD); + name = "GAUSS_CHEBYSHEV_2MOD"; + } + + SECTION("Gauss-Chebyshev 3") { + error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_3); + name = "GAUSS_CHEBYSHEV_3"; + } + + REQUIRE(error == INTXX_SUCCESS); + + // Meta data + REQUIRE(quad.info != NULL); + REQUIRE(quad.npoints == -1); + REQUIRE(quad.info->ext_params.n == 0); + REQUIRE(!strcmp(quad.info->name, name)); + + // Set NPTS + error = intxx_quad_set_npts(&quad, base_npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(quad.npoints == base_npts); + + // Get NPTS + int npts; + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(npts == base_npts); + } + + intxx_quad_end(&quad); + +} From 3cfb79a2aee278767b9813eafb5af8f72b4d53ec Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 23 Sep 2023 17:20:14 -0700 Subject: [PATCH 02/22] Dox++ --- include/integratorxx/c_api.h | 57 ++++++++++++++++++++++++++++++++++-- src/c_api/c_api.c | 2 +- src/c_api/c_radial.cxx | 2 +- test/c_api.cxx | 7 ++++- 4 files changed, 62 insertions(+), 6 deletions(-) diff --git a/include/integratorxx/c_api.h b/include/integratorxx/c_api.h index 281587f..766af76 100644 --- a/include/integratorxx/c_api.h +++ b/include/integratorxx/c_api.h @@ -7,6 +7,7 @@ #define INTXX_NULL_INFOPTR -3 #define INTXX_INVALID_OPT -4 #define INTXX_INVALID_ARG -5 +#define INTXX_INVALID_OUT -5 /*** Quadrature Classes ***/ #define INTXX_PRM_QUAD 1 @@ -52,7 +53,7 @@ typedef struct { const char** descriptions; ///< Long descriptions of params const double *values; ///< Default values of params - void (*set)(struct intxx_quad_type* p, const double** v); + //void (*set)(struct intxx_quad_type* p, const double** v); ///< Set function } intxx_quad_params_type; @@ -65,7 +66,7 @@ typedef struct { intxx_quad_params_type ext_params; ///< External params - void (*init)(struct intxx_quad_type* p); + //void (*init)(struct intxx_quad_type* p); } intxx_quad_info_type; struct intxx_quad_type { @@ -77,11 +78,61 @@ 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); #ifdef __cplusplus diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c index 548a0ec..3293311 100644 --- a/src/c_api/c_api.c +++ b/src/c_api/c_api.c @@ -82,6 +82,6 @@ int intxx_quad_get_npts(intxx_quad_type* p, int* npts) { // TODO: Handle the case when NPTS is derived from params *npts = p->npoints; - return INTXX_SUCCESS; + return *npts > 0 ? INTXX_SUCCESS : INTXX_INVALID_OUT; } diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index 6571643..7fb56a1 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -15,7 +15,7 @@ void intxx_default_quad_info(intxx_quad_info_type* p) { p->number = 0; // Invalid p->kind = 0; // Invalid p->name = "Default"; // Default state - p->init = NULL; + //p->init = NULL; } int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { diff --git a/test/c_api.cxx b/test/c_api.cxx index 79c100a..498fc01 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -62,13 +62,18 @@ TEST_CASE("C API") { REQUIRE(quad.info->ext_params.n == 0); REQUIRE(!strcmp(quad.info->name, name)); + // Get before set + int npts; + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_INVALID_OUT); + REQUIRE(npts == -1); + // Set NPTS error = intxx_quad_set_npts(&quad, base_npts); REQUIRE(error == INTXX_SUCCESS); REQUIRE(quad.npoints == base_npts); // Get NPTS - int npts; error = intxx_quad_get_npts(&quad, &npts); REQUIRE(error == INTXX_SUCCESS); REQUIRE(npts == base_npts); From 055d8226451ffd4683ecc69e107c646dda26f500 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 23 Sep 2023 18:48:33 -0700 Subject: [PATCH 03/22] Added ability to generate and destroy quadratures, tested with UniformTrapezoid --- include/integratorxx/c_api.h | 8 +++- src/c_api/c_api.c | 14 ++++++ src/c_api/c_radial.cxx | 88 ++++++++++++++++++++++++++++++++---- test/c_api.cxx | 27 +++++++++++ 4 files changed, 127 insertions(+), 10 deletions(-) diff --git a/include/integratorxx/c_api.h b/include/integratorxx/c_api.h index 766af76..1ad100e 100644 --- a/include/integratorxx/c_api.h +++ b/include/integratorxx/c_api.h @@ -66,12 +66,14 @@ typedef struct { intxx_quad_params_type ext_params; ///< External params - //void (*init)(struct intxx_quad_type* p); + int (*generate)(struct intxx_quad_type* p); + int (*destroy) (struct intxx_quad_type* p); } intxx_quad_info_type; struct intxx_quad_type { int npoints; const intxx_quad_info_type* info; + void* _state; ///< Internal state, NOT FOR PUBLIC USE }; typedef struct intxx_quad_type intxx_quad_type; @@ -135,6 +137,10 @@ int intxx_quad_set_npts(intxx_quad_type* p, int npts); */ 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); + #ifdef __cplusplus } #endif diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c index 3293311..c7f08ff 100644 --- a/src/c_api/c_api.c +++ b/src/c_api/c_api.c @@ -11,6 +11,7 @@ int intxx_quad_init(intxx_quad_type* p, int quad) { // Sanity check if(p == NULL) return INTXX_NULL_QUADPTR; p->info = NULL; + p->_state = NULL; p->npoints = -1; if(quad < 0) return INTXX_INVALID_QUAD; @@ -84,4 +85,17 @@ int intxx_quad_get_npts(intxx_quad_type* p, int* npts) { *npts = p->npoints; return *npts > 0 ? INTXX_SUCCESS : INTXX_INVALID_OUT; } + +int intxx_generate_quad(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + return p->info->generate(p); +} +int intxx_destroy_quad(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + return p->info->destroy(p); +} diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index 7fb56a1..5528c4f 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -1,3 +1,5 @@ +#include + #include "c_internal.h" #include @@ -11,11 +13,15 @@ int intxx_get_gausscheb2_info(intxx_quad_info_type* p); int intxx_get_gausscheb2m_info(intxx_quad_info_type* p); int intxx_get_gausscheb3_info(intxx_quad_info_type* p); +int intxx_generate_uniform(intxx_quad_type* p); +int intxx_destroy_uniform(intxx_quad_type* p); + void intxx_default_quad_info(intxx_quad_info_type* p) { p->number = 0; // Invalid p->kind = 0; // Invalid p->name = "Default"; // Default state - //p->init = NULL; + p->generate = NULL; + p->destroy = NULL; } int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { @@ -51,34 +57,98 @@ int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { } } -int intxx_noparam_info(intxx_quad_info_type* p) { +int intxx_noparam_info(intxx_quad_info_type* p, + int (*g)(intxx_quad_type*), + int (*d)(intxx_quad_type*)) { if(p == NULL) return INTXX_NULL_INFOPTR; p->ext_params.n = 0; /// No External Parameters + p->generate = g; + p->destroy = d; return INTXX_SUCCESS; } int intxx_get_uniform_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p); + return intxx_noparam_info(p, &intxx_generate_uniform, + &intxx_destroy_uniform); } int intxx_get_gausslob_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p); + return intxx_noparam_info(p, NULL, NULL); } int intxx_get_gaussleg_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p); + return intxx_noparam_info(p, NULL, NULL); } int intxx_get_gausscheb1_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p); + return intxx_noparam_info(p, NULL, NULL); } int intxx_get_gausscheb2_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p); + return intxx_noparam_info(p, NULL, NULL); } int intxx_get_gausscheb2m_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p); + return intxx_noparam_info(p, NULL, NULL); } int intxx_get_gausscheb3_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p); + return intxx_noparam_info(p, NULL, NULL); +} + + + +int intxx_generate_uniform(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + using namespace IntegratorXX; + using quad_type = UniformTrapezoid; + using alloc_type = std::allocator; + using alloc_traits = std::allocator_traits; + + int npts = p->npoints; + if( npts <= 0 ) { + // Return error code + } + + if(p->_state != NULL) { + // Check if we need to regenerate + } + + + // Allocate and construct the quadrature instance + alloc_type alloc; + auto ptr = alloc_traits::allocate(alloc, 1); + alloc_traits::construct(alloc, ptr, npts, 0.0, 1.0); + + // Store the pointer + p->_state = (void*)ptr; + + return INTXX_SUCCESS; + +} + +int intxx_destroy_uniform(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + using namespace IntegratorXX; + using quad_type = UniformTrapezoid; + using alloc_type = std::allocator; + using alloc_traits = std::allocator_traits; + + if(p->_state != NULL) { + alloc_type alloc; + + // Cross your fingers... + auto ptr = p->_state; + auto quad_ptr = reinterpret_cast(ptr); + + // Destroy and deallocate + alloc_traits::destroy(alloc, quad_ptr); + alloc_traits::deallocate(alloc, quad_ptr, 1); + + // Null out state + p->_state = NULL; + } + return INTXX_SUCCESS; } } diff --git a/test/c_api.cxx b/test/c_api.cxx index 498fc01..23fd297 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -1,7 +1,11 @@ #include "catch2/catch_all.hpp" +#include "quad_matcher.hpp" #include #include + +using namespace IntegratorXX; + TEST_CASE("C API") { intxx_quad_type quad; @@ -18,10 +22,15 @@ TEST_CASE("C API") { const char* name; const int base_npts = 100; + using base_quad_type = QuadratureBase, std::vector>; + + std::unique_ptr base_quad = nullptr; SECTION("Uniform") { + using quad_type = UniformTrapezoid; error = intxx_quad_init(&quad, INTXX_PRMQ_UNIFORM); name = "UNIFORM"; + base_quad = std::make_unique(base_npts, 0.0, 1.0); } SECTION("Gauss-Legendre") { @@ -59,6 +68,7 @@ TEST_CASE("C API") { // Meta data REQUIRE(quad.info != NULL); REQUIRE(quad.npoints == -1); + REQUIRE(quad._state == NULL); REQUIRE(quad.info->ext_params.n == 0); REQUIRE(!strcmp(quad.info->name, name)); @@ -77,6 +87,23 @@ TEST_CASE("C API") { error = intxx_quad_get_npts(&quad, &npts); REQUIRE(error == INTXX_SUCCESS); REQUIRE(npts == base_npts); + + // Check Quadrature Generation and Destruction + if(base_quad) { + intxx_generate_quad(&quad); + REQUIRE(quad._state != NULL); + + // Check validity of the state + auto state_as_quad = reinterpret_cast(quad._state); + REQUIRE(state_as_quad->npts() == npts); + for(auto i = 0; i < npts; ++ i) { + REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad->points()[i], 1e-15)); + REQUIRE_THAT(state_as_quad->weights()[i], Matchers::WithinAbs(name, base_quad->weights()[i], 1e-15)); + } + + intxx_destroy_quad(&quad); + REQUIRE(quad._state == NULL); + } } intxx_quad_end(&quad); From 504312e21cc0a32d91157cd07d0013fe5d96f92f Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sat, 23 Sep 2023 19:07:03 -0700 Subject: [PATCH 04/22] Minor refactor of C impls --- src/c_api/c_radial.cxx | 98 ++++++++++++++++++++++++------------------ 1 file changed, 57 insertions(+), 41 deletions(-) diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index 5528c4f..c596155 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -3,6 +3,59 @@ #include "c_internal.h" #include +template +auto generate_quadrature(Args&&... args) { + using alloc_type = std::allocator; + using alloc_traits = std::allocator_traits; + + alloc_type alloc; + auto ptr = alloc_traits::allocate(alloc, 1); + alloc_traits::construct(alloc, ptr, std::forward(args)...); + + return ptr; +} + +template +int intxx_generate_impl(intxx_quad_type* p, Args&&... args) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + using quad_type = QuadType; + + int npts = p->npoints; + if( npts <= 0 ) { + // Return error code + } + + if(p->_state != NULL) { + // Check if we need to regenerate + } + + // Allocate and construct the quadrature instance + auto ptr = generate_quadrature(npts, std::forward(args)... ); + + // Store the pointer + p->_state = (void*)ptr; + + return INTXX_SUCCESS; + +} + +template +void destroy_quadrature(void* ptr) { + using alloc_type = std::allocator; + using alloc_traits = std::allocator_traits; + + alloc_type alloc; + + // Cross your fingers... + auto quad_ptr = reinterpret_cast(ptr); + + // Destroy and deallocate + alloc_traits::destroy(alloc, quad_ptr); + alloc_traits::deallocate(alloc, quad_ptr, 1); +} + extern "C" { int intxx_get_uniform_info(intxx_quad_info_type* p); @@ -93,37 +146,9 @@ int intxx_get_gausscheb3_info(intxx_quad_info_type* p) { } - int intxx_generate_uniform(intxx_quad_type* p) { - if(p == NULL) return INTXX_NULL_QUADPTR; - if(p->info == NULL) return INTXX_NULL_INFOPTR; - - using namespace IntegratorXX; - using quad_type = UniformTrapezoid; - using alloc_type = std::allocator; - using alloc_traits = std::allocator_traits; - - int npts = p->npoints; - if( npts <= 0 ) { - // Return error code - } - - if(p->_state != NULL) { - // Check if we need to regenerate - } - - - // Allocate and construct the quadrature instance - alloc_type alloc; - auto ptr = alloc_traits::allocate(alloc, 1); - alloc_traits::construct(alloc, ptr, npts, 0.0, 1.0); - - // Store the pointer - p->_state = (void*)ptr; - - return INTXX_SUCCESS; - -} + return intxx_generate_impl>(p, 0.0, 1.0); +} int intxx_destroy_uniform(intxx_quad_type* p) { if(p == NULL) return INTXX_NULL_QUADPTR; @@ -131,19 +156,10 @@ int intxx_destroy_uniform(intxx_quad_type* p) { using namespace IntegratorXX; using quad_type = UniformTrapezoid; - using alloc_type = std::allocator; - using alloc_traits = std::allocator_traits; if(p->_state != NULL) { - alloc_type alloc; - - // Cross your fingers... - auto ptr = p->_state; - auto quad_ptr = reinterpret_cast(ptr); - - // Destroy and deallocate - alloc_traits::destroy(alloc, quad_ptr); - alloc_traits::deallocate(alloc, quad_ptr, 1); + // Destroy quadrature + destroy_quadrature(p->_state); // Null out state p->_state = NULL; From 10a0da38741c6ae15b2db7237095b976f88411bc Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 10:23:08 -0700 Subject: [PATCH 05/22] Cleanup state on intxx_quad_end if required --- src/c_api/c_api.c | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c index c7f08ff..00eb1e8 100644 --- a/src/c_api/c_api.c +++ b/src/c_api/c_api.c @@ -26,7 +26,8 @@ int intxx_quad_init(intxx_quad_type* p, int quad) { if(!is_prmq && !is_radq && !is_angq) return INTXX_INVALID_QUAD; - // Primitive quadratures cannot be mixed with angular ones + // Primitive quadratures cannot be mixed with angular or + // radial quadratures if(is_prmq && is_radq) return INTXX_INVALID_QUAD; if(is_prmq && is_angq) return INTXX_INVALID_QUAD; @@ -55,6 +56,9 @@ void intxx_quad_end(intxx_quad_type* p) { // Stateless - just bail if(p == NULL || p->info == NULL) return; + // Destroy quadrature state if populated + intxx_destroy_quad(p); + // Free info free((void*)p->info); p->info = NULL; @@ -90,12 +94,16 @@ int intxx_generate_quad(intxx_quad_type* p) { if(p == NULL) return INTXX_NULL_QUADPTR; if(p->info == NULL) return INTXX_NULL_INFOPTR; - return p->info->generate(p); + if(p->info->generate) + return p->info->generate(p); + else return INTXX_SUCCESS; } int intxx_destroy_quad(intxx_quad_type* p) { if(p == NULL) return INTXX_NULL_QUADPTR; if(p->info == NULL) return INTXX_NULL_INFOPTR; - return p->info->destroy(p); + if(p->info->destroy) + return p->info->destroy(p); + else return INTXX_SUCCESS; } From cdd4d0ea62d9c4e8cb29b7ab73603d9d8acb5f98 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 10:37:29 -0700 Subject: [PATCH 06/22] Added C API for Gauss{Legendre,Lobatto} --- src/c_api/c_radial.cxx | 64 +++++++++++++++++++++++++++++------------- test/c_api.cxx | 4 +++ 2 files changed, 49 insertions(+), 19 deletions(-) diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index c596155..516fd96 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -1,4 +1,4 @@ -#include +#include #include "c_internal.h" #include @@ -56,6 +56,24 @@ void destroy_quadrature(void* ptr) { alloc_traits::deallocate(alloc, quad_ptr, 1); } +template +int intxx_destroy_impl(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + using namespace IntegratorXX; + using quad_type = QuadType; + + if(p->_state != NULL) { + // Destroy quadrature + destroy_quadrature(p->_state); + + // Null out state + p->_state = NULL; + } + return INTXX_SUCCESS; +} + extern "C" { int intxx_get_uniform_info(intxx_quad_info_type* p); @@ -68,6 +86,10 @@ int intxx_get_gausscheb3_info(intxx_quad_info_type* p); int intxx_generate_uniform(intxx_quad_type* p); int intxx_destroy_uniform(intxx_quad_type* p); +int intxx_generate_gaussleg(intxx_quad_type* p); +int intxx_destroy_gaussleg(intxx_quad_type* p); +int intxx_generate_gausslob(intxx_quad_type* p); +int intxx_destroy_gausslob(intxx_quad_type* p); void intxx_default_quad_info(intxx_quad_info_type* p) { p->number = 0; // Invalid @@ -127,10 +149,12 @@ int intxx_get_uniform_info(intxx_quad_info_type* p) { &intxx_destroy_uniform); } int intxx_get_gausslob_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, NULL, NULL); + return intxx_noparam_info(p, &intxx_generate_gausslob, + &intxx_destroy_gausslob); } int intxx_get_gaussleg_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, NULL, NULL); + return intxx_noparam_info(p, &intxx_generate_gaussleg, + &intxx_destroy_gaussleg); } int intxx_get_gausscheb1_info(intxx_quad_info_type* p) { return intxx_noparam_info(p, NULL, NULL); @@ -145,26 +169,28 @@ int intxx_get_gausscheb3_info(intxx_quad_info_type* p) { return intxx_noparam_info(p, NULL, NULL); } - +using uniform_quad_type = IntegratorXX::UniformTrapezoid; int intxx_generate_uniform(intxx_quad_type* p) { - return intxx_generate_impl>(p, 0.0, 1.0); + return intxx_generate_impl(p, 0.0, 1.0); } - int intxx_destroy_uniform(intxx_quad_type* p) { - if(p == NULL) return INTXX_NULL_QUADPTR; - if(p->info == NULL) return INTXX_NULL_INFOPTR; - - using namespace IntegratorXX; - using quad_type = UniformTrapezoid; + return intxx_destroy_impl(p); +} - if(p->_state != NULL) { - // Destroy quadrature - destroy_quadrature(p->_state); - - // Null out state - p->_state = NULL; - } - return INTXX_SUCCESS; + +#define INTXX_GD_BASIC_IMPL(cname, qname) \ +int intxx_generate_##cname(intxx_quad_type* p) { \ + return intxx_generate_impl(p); \ +} \ +int intxx_destroy_##cname(intxx_quad_type* p) { \ + return intxx_destroy_impl(p); \ } +using gaussleg_quad_type = IntegratorXX::GaussLegendre; +using gausslob_quad_type = IntegratorXX::GaussLobatto; +INTXX_GD_BASIC_IMPL(gaussleg, gaussleg_quad_type); +INTXX_GD_BASIC_IMPL(gausslob, gausslob_quad_type); + +#undef INTXX_GD_BASIC_IMPL + } diff --git a/test/c_api.cxx b/test/c_api.cxx index 23fd297..82a5d5a 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -34,13 +34,17 @@ TEST_CASE("C API") { } SECTION("Gauss-Legendre") { + using quad_type = GaussLegendre; error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSLEG); name = "GAUSS_LEGENDRE"; + base_quad = std::make_unique(base_npts); } SECTION("Gauss-Lobatto") { + using quad_type = GaussLobatto; error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSLOB); name = "GAUSS_LOBATTO"; + base_quad = std::make_unique(base_npts); } SECTION("Gauss-Chebyshev 1") { From 43d9e1b1c7d095083d62ad4e81f36cf94bbe05dc Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 10:45:04 -0700 Subject: [PATCH 07/22] Added C API for GaussChebyshev{1,2,2Modified,3} --- src/c_api/c_radial.cxx | 41 +++++++++++++++++++++++++++++++---------- test/c_api.cxx | 8 ++++++++ 2 files changed, 39 insertions(+), 10 deletions(-) diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index 516fd96..1b5b94f 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -90,13 +90,21 @@ int intxx_generate_gaussleg(intxx_quad_type* p); int intxx_destroy_gaussleg(intxx_quad_type* p); int intxx_generate_gausslob(intxx_quad_type* p); int intxx_destroy_gausslob(intxx_quad_type* p); +int intxx_generate_gausscheb1(intxx_quad_type* p); +int intxx_destroy_gausscheb1(intxx_quad_type* p); +int intxx_generate_gausscheb2(intxx_quad_type* p); +int intxx_destroy_gausscheb2(intxx_quad_type* p); +int intxx_generate_gausscheb3(intxx_quad_type* p); +int intxx_destroy_gausscheb3(intxx_quad_type* p); +int intxx_generate_gausscheb2mod(intxx_quad_type* p); +int intxx_destroy_gausscheb2mod(intxx_quad_type* p); void intxx_default_quad_info(intxx_quad_info_type* p) { p->number = 0; // Invalid p->kind = 0; // Invalid p->name = "Default"; // Default state - p->generate = NULL; - p->destroy = NULL; + p->generate = NULL; // Invalid + p->destroy = NULL; // Invalid } int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { @@ -157,16 +165,20 @@ int intxx_get_gaussleg_info(intxx_quad_info_type* p) { &intxx_destroy_gaussleg); } int intxx_get_gausscheb1_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, NULL, NULL); + return intxx_noparam_info(p, &intxx_generate_gausscheb1, + &intxx_destroy_gausscheb1); } int intxx_get_gausscheb2_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, NULL, NULL); + return intxx_noparam_info(p, &intxx_generate_gausscheb2, + &intxx_destroy_gausscheb2); } int intxx_get_gausscheb2m_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, NULL, NULL); + return intxx_noparam_info(p, &intxx_generate_gausscheb2mod, + &intxx_destroy_gausscheb2mod); } int intxx_get_gausscheb3_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, NULL, NULL); + return intxx_noparam_info(p, &intxx_generate_gausscheb3, + &intxx_destroy_gausscheb3); } using uniform_quad_type = IntegratorXX::UniformTrapezoid; @@ -186,10 +198,19 @@ int intxx_destroy_##cname(intxx_quad_type* p) { \ return intxx_destroy_impl(p); \ } -using gaussleg_quad_type = IntegratorXX::GaussLegendre; -using gausslob_quad_type = IntegratorXX::GaussLobatto; -INTXX_GD_BASIC_IMPL(gaussleg, gaussleg_quad_type); -INTXX_GD_BASIC_IMPL(gausslob, gausslob_quad_type); +using gaussleg_quad_type = IntegratorXX::GaussLegendre; +using gausslob_quad_type = IntegratorXX::GaussLobatto; +using gausscheb1_quad_type = IntegratorXX::GaussChebyshev1; +using gausscheb2_quad_type = IntegratorXX::GaussChebyshev2; +using gausscheb3_quad_type = IntegratorXX::GaussChebyshev3; +using gausscheb2mod_quad_type = IntegratorXX::GaussChebyshev2Modified; + +INTXX_GD_BASIC_IMPL(gaussleg, gaussleg_quad_type ); +INTXX_GD_BASIC_IMPL(gausslob, gausslob_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb1, gausscheb1_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb2, gausscheb2_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb3, gausscheb3_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb2mod, gausscheb2mod_quad_type); #undef INTXX_GD_BASIC_IMPL diff --git a/test/c_api.cxx b/test/c_api.cxx index 82a5d5a..c287af3 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -48,23 +48,31 @@ TEST_CASE("C API") { } SECTION("Gauss-Chebyshev 1") { + using quad_type = GaussChebyshev1; error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_1); name = "GAUSS_CHEBYSHEV_1"; + base_quad = std::make_unique(base_npts); } SECTION("Gauss-Chebyshev 2") { + using quad_type = GaussChebyshev2; error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_2); name = "GAUSS_CHEBYSHEV_2"; + base_quad = std::make_unique(base_npts); } SECTION("Gauss-Chebyshev 2MOD") { + using quad_type = GaussChebyshev2Modified; error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_2MOD); name = "GAUSS_CHEBYSHEV_2MOD"; + base_quad = std::make_unique(base_npts); } SECTION("Gauss-Chebyshev 3") { + using quad_type = GaussChebyshev3; error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_3); name = "GAUSS_CHEBYSHEV_3"; + base_quad = std::make_unique(base_npts); } REQUIRE(error == INTXX_SUCCESS); From 0af91e537dee29acfa7f2df8d4b24f1bbf39fb1c Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 11:05:43 -0700 Subject: [PATCH 08/22] Radial name added to C API and tested --- src/c_api/c_api.c | 1 + src/c_api/c_radial.cxx | 26 ++++++++++++++++++ test/c_api.cxx | 62 +++++++++++++++++++++++++++++++++++++++++- 3 files changed, 88 insertions(+), 1 deletion(-) diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c index 00eb1e8..94bd73f 100644 --- a/src/c_api/c_api.c +++ b/src/c_api/c_api.c @@ -44,6 +44,7 @@ int intxx_quad_init(intxx_quad_type* p, int quad) { // TODO: Get spherical quadrature info } else if(is_radq) { // TODO: Get radial quadrature info + error = intxx_get_radq_info(finfo, quad); } else { // Angular by exclusion // TODO: Get angular quadrature info diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index 1b5b94f..624e6d3 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -140,6 +140,32 @@ int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { } } + +int intxx_get_radq_info(intxx_quad_info_type* p, int quad) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + quad = quad & INTXX_RADQ_MASK; + p->number = quad; + p->kind = INTXX_RAD_QUAD; + switch(quad) { + case INTXX_RADQ_BECKE: + p->name = "BECKE"; + break; + case INTXX_RADQ_MHL: + p->name = "MURRAY_HANDY_LAMING"; + break; + case INTXX_RADQ_TA: + p->name = "TREUTLER_AHLRICHS"; + break; + case INTXX_RADQ_MK: + p->name = "MURA_KNOWLES"; + break; + default: + return INTXX_INVALID_QUAD; + } + return INTXX_SUCCESS; +} + int intxx_noparam_info(intxx_quad_info_type* p, int (*g)(intxx_quad_type*), int (*d)(intxx_quad_type*)) { diff --git a/test/c_api.cxx b/test/c_api.cxx index c287af3..971d617 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -18,12 +18,12 @@ TEST_CASE("C API") { REQUIRE(quad.npoints == -1); } + SECTION("Primitive") { const char* name; const int base_npts = 100; using base_quad_type = QuadratureBase, std::vector>; - std::unique_ptr base_quad = nullptr; SECTION("Uniform") { @@ -118,6 +118,66 @@ TEST_CASE("C API") { } } + + SECTION("Radial") { + const char* name; + const int base_npts = 100; + using base_quad_type = QuadratureBase, std::vector>; + std::unique_ptr base_quad_default = nullptr; + + SECTION("Becke") { + //using quad_type = UniformTrapezoid; + error = intxx_quad_init(&quad, INTXX_RADQ_BECKE); + name = "BECKE"; + //base_quad = std::make_unique(base_npts); + } + + SECTION("MHL") { + //using quad_type = UniformTrapezoid; + error = intxx_quad_init(&quad, INTXX_RADQ_MHL); + name = "MURRAY_HANDY_LAMING"; + //base_quad = std::make_unique(base_npts); + } + + SECTION("TA") { + //using quad_type = UniformTrapezoid; + error = intxx_quad_init(&quad, INTXX_RADQ_TA); + name = "TREUTLER_AHLRICHS"; + //base_quad = std::make_unique(base_npts); + } + + SECTION("MK") { + //using quad_type = UniformTrapezoid; + error = intxx_quad_init(&quad, INTXX_RADQ_MK); + name = "MURA_KNOWLES"; + //base_quad = std::make_unique(base_npts); + } + + REQUIRE(error == INTXX_SUCCESS); + + // Meta data + REQUIRE(quad.info != NULL); + REQUIRE(quad.npoints == -1); + REQUIRE(quad._state == NULL); + REQUIRE(!strcmp(quad.info->name, name)); + + // Get before set + int npts; + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_INVALID_OUT); + REQUIRE(npts == -1); + + // Set NPTS + error = intxx_quad_set_npts(&quad, base_npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(quad.npoints == base_npts); + + // Get NPTS + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(npts == base_npts); + } + intxx_quad_end(&quad); } From ff0f04c0b4f612aaf78178e803fdab6c7f77c133 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 13:38:36 -0700 Subject: [PATCH 09/22] Added generation of Radial quadratures to C API, includes storage of traits classes --- include/integratorxx/c_api.h | 12 +- .../quadratures/radial/radial_transform.hpp | 3 + src/c_api/c_api.c | 16 +- src/c_api/c_radial.cxx | 297 +++++++++++++----- test/c_api.cxx | 47 ++- 5 files changed, 274 insertions(+), 101 deletions(-) diff --git a/include/integratorxx/c_api.h b/include/integratorxx/c_api.h index 1ad100e..a32b3a4 100644 --- a/include/integratorxx/c_api.h +++ b/include/integratorxx/c_api.h @@ -51,10 +51,13 @@ typedef struct { const char** names; ///< Names of the params const char** descriptions; ///< Long descriptions of params - const double *values; ///< Default values of params - //void (*set)(struct intxx_quad_type* p, const double** v); - ///< Set function + void (*set_name)(struct intxx_quad_type* p, const char* name, double v); + ///< Set parameter function + void (*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 { @@ -73,7 +76,8 @@ typedef struct { struct intxx_quad_type { int npoints; const intxx_quad_info_type* info; - void* _state; ///< Internal state, NOT FOR PUBLIC USE + 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; diff --git a/include/integratorxx/quadratures/radial/radial_transform.hpp b/include/integratorxx/quadratures/radial/radial_transform.hpp index 93f8096..418d00b 100644 --- a/include/integratorxx/quadratures/radial/radial_transform.hpp +++ b/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; diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c index 94bd73f..00e607d 100644 --- a/src/c_api/c_api.c +++ b/src/c_api/c_api.c @@ -10,8 +10,9 @@ int intxx_quad_init(intxx_quad_type* p, int quad) { // Sanity check if(p == NULL) return INTXX_NULL_QUADPTR; - p->info = NULL; - p->_state = NULL; + p->info = NULL; + p->_state_quad = NULL; + p->_state_parm = NULL; p->npoints = -1; if(quad < 0) return INTXX_INVALID_QUAD; @@ -50,6 +51,12 @@ int intxx_quad_init(intxx_quad_type* p, int quad) { // TODO: Get angular quadrature info } + if(error) return error; + + if(finfo->ext_params.n && finfo->ext_params.generate) { + error = finfo->ext_params.generate(p); + } + return error; } @@ -57,6 +64,11 @@ void intxx_quad_end(intxx_quad_type* p) { // Stateless - just bail if(p == NULL || p->info == NULL) return; + // Destroy traits if required + if(p->info->ext_params.n && p->info->ext_params.destroy) { + p->info->ext_params.destroy(p); + } + // Destroy quadrature state if populated intxx_destroy_quad(p); diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index 624e6d3..3df70b3 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -1,11 +1,12 @@ #include +#include #include "c_internal.h" #include -template +template auto generate_quadrature(Args&&... args) { - using alloc_type = std::allocator; + using alloc_type = std::allocator; using alloc_traits = std::allocator_traits; alloc_type alloc; @@ -16,7 +17,7 @@ auto generate_quadrature(Args&&... args) { } template -int intxx_generate_impl(intxx_quad_type* p, Args&&... args) { +int intxx_generate_noparam_impl(intxx_quad_type* p, Args&&... args) { if(p == NULL) return INTXX_NULL_QUADPTR; if(p->info == NULL) return INTXX_NULL_INFOPTR; @@ -27,7 +28,7 @@ int intxx_generate_impl(intxx_quad_type* p, Args&&... args) { // Return error code } - if(p->_state != NULL) { + if(p->_state_quad != NULL) { // Check if we need to regenerate } @@ -35,21 +36,36 @@ int intxx_generate_impl(intxx_quad_type* p, Args&&... args) { auto ptr = generate_quadrature(npts, std::forward(args)... ); // Store the pointer - p->_state = (void*)ptr; + p->_state_quad = (void*)ptr; return INTXX_SUCCESS; } -template + +template +int intxx_generate_param_impl(intxx_quad_type* p, Args&&... args) { + if(p == NULL) return INTXX_NULL_QUADPTR; + + if(p->_state_parm == NULL) { + // Return error code + } + + // Cross your fingers... + const auto _params = reinterpret_cast(p->_state_parm); + return intxx_generate_noparam_impl(p, *_params, std::forward(args)... ); + +} + +template void destroy_quadrature(void* ptr) { - using alloc_type = std::allocator; + using alloc_type = std::allocator; using alloc_traits = std::allocator_traits; alloc_type alloc; // Cross your fingers... - auto quad_ptr = reinterpret_cast(ptr); + auto quad_ptr = reinterpret_cast(ptr); // Destroy and deallocate alloc_traits::destroy(alloc, quad_ptr); @@ -59,45 +75,87 @@ void destroy_quadrature(void* ptr) { template int intxx_destroy_impl(intxx_quad_type* p) { if(p == NULL) return INTXX_NULL_QUADPTR; - if(p->info == NULL) return INTXX_NULL_INFOPTR; - using namespace IntegratorXX; using quad_type = QuadType; - if(p->_state != NULL) { + if(p->_state_quad != NULL) { // Destroy quadrature - destroy_quadrature(p->_state); + destroy_quadrature(p->_state_quad); + + // Null out state + p->_state_quad = NULL; + } + return INTXX_SUCCESS; +} + +template +int intxx_generate_params_impl(intxx_quad_type* p, Args&&... args) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(p->_state_parm != NULL) { + return INTXX_SUCCESS; + } + + // Allocate and construct the quadrature instance + auto ptr = generate_quadrature(std::forward(args)... ); + + // Store the pointer + p->_state_parm = (void*)ptr; + + return INTXX_SUCCESS; + +} + +template +int intxx_destroy_params_impl(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + + if(p->_state_parm != NULL) { + // Destroy traits object + destroy_quadrature(p->_state_parm); // Null out state - p->_state = NULL; + p->_state_parm = NULL; } return INTXX_SUCCESS; } + extern "C" { -int intxx_get_uniform_info(intxx_quad_info_type* p); -int intxx_get_gaussleg_info(intxx_quad_info_type* p); -int intxx_get_gausslob_info(intxx_quad_info_type* p); -int intxx_get_gausscheb1_info(intxx_quad_info_type* p); -int intxx_get_gausscheb2_info(intxx_quad_info_type* p); -int intxx_get_gausscheb2m_info(intxx_quad_info_type* p); -int intxx_get_gausscheb3_info(intxx_quad_info_type* p); - -int intxx_generate_uniform(intxx_quad_type* p); -int intxx_destroy_uniform(intxx_quad_type* p); -int intxx_generate_gaussleg(intxx_quad_type* p); -int intxx_destroy_gaussleg(intxx_quad_type* p); -int intxx_generate_gausslob(intxx_quad_type* p); -int intxx_destroy_gausslob(intxx_quad_type* p); -int intxx_generate_gausscheb1(intxx_quad_type* p); -int intxx_destroy_gausscheb1(intxx_quad_type* p); -int intxx_generate_gausscheb2(intxx_quad_type* p); -int intxx_destroy_gausscheb2(intxx_quad_type* p); -int intxx_generate_gausscheb3(intxx_quad_type* p); -int intxx_destroy_gausscheb3(intxx_quad_type* p); -int intxx_generate_gausscheb2mod(intxx_quad_type* p); -int intxx_destroy_gausscheb2mod(intxx_quad_type* p); +#define FWD_PRM(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_PRM(uniform) +FWD_PRM(gaussleg) +FWD_PRM(gausslob) +FWD_PRM(gausscheb1) +FWD_PRM(gausscheb2) +FWD_PRM(gausscheb2m) +FWD_PRM(gausscheb3) + +#undef FWD_PRM + +#define FWD_RAD(cname) \ +int intxx_get_##cname##_quad_info(intxx_quad_info_type* p); \ +int intxx_generate_##cname##_quad(intxx_quad_type* p); \ +int intxx_destroy_##cname##_quad(intxx_quad_type* p); \ +int intxx_generate_##cname##_quad_params(intxx_quad_type* p); \ +int intxx_destroy_##cname##_quad_params(intxx_quad_type* p); + +FWD_RAD(becke) +FWD_RAD(mhl) +FWD_RAD(ta) +FWD_RAD(mk) + +#undef FWD_RAD + + + + void intxx_default_quad_info(intxx_quad_info_type* p) { p->number = 0; // Invalid @@ -150,22 +208,26 @@ int intxx_get_radq_info(intxx_quad_info_type* p, int quad) { switch(quad) { case INTXX_RADQ_BECKE: p->name = "BECKE"; - break; + return intxx_get_becke_quad_info(p); case INTXX_RADQ_MHL: p->name = "MURRAY_HANDY_LAMING"; - break; + return intxx_get_mhl_quad_info(p); case INTXX_RADQ_TA: p->name = "TREUTLER_AHLRICHS"; - break; + return intxx_get_ta_quad_info(p); case INTXX_RADQ_MK: p->name = "MURA_KNOWLES"; - break; + return intxx_get_mk_quad_info(p); default: return INTXX_INVALID_QUAD; } return INTXX_SUCCESS; } +/******************************************************/ +/****** Info generation for quads w/o parameters ******/ +/******************************************************/ + int intxx_noparam_info(intxx_quad_info_type* p, int (*g)(intxx_quad_type*), int (*d)(intxx_quad_type*)) { @@ -178,66 +240,137 @@ int intxx_noparam_info(intxx_quad_info_type* p, return INTXX_SUCCESS; } -int intxx_get_uniform_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, &intxx_generate_uniform, - &intxx_destroy_uniform); -} -int intxx_get_gausslob_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, &intxx_generate_gausslob, - &intxx_destroy_gausslob); -} -int intxx_get_gaussleg_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, &intxx_generate_gaussleg, - &intxx_destroy_gaussleg); -} -int intxx_get_gausscheb1_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, &intxx_generate_gausscheb1, - &intxx_destroy_gausscheb1); -} -int intxx_get_gausscheb2_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, &intxx_generate_gausscheb2, - &intxx_destroy_gausscheb2); -} -int intxx_get_gausscheb2m_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, &intxx_generate_gausscheb2mod, - &intxx_destroy_gausscheb2mod); -} -int intxx_get_gausscheb3_info(intxx_quad_info_type* p) { - return intxx_noparam_info(p, &intxx_generate_gausscheb3, - &intxx_destroy_gausscheb3); +#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); \ } -using uniform_quad_type = IntegratorXX::UniformTrapezoid; -int intxx_generate_uniform(intxx_quad_type* p) { - return intxx_generate_impl(p, 0.0, 1.0); +INTXX_NOPARAM_GET_INFO_IMPL(uniform); +INTXX_NOPARAM_GET_INFO_IMPL(gaussleg); +INTXX_NOPARAM_GET_INFO_IMPL(gausslob); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb1); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb2); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb2m); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb3); + +#undef INTXX_NOPARAM_GET_INFO_IMPL + +/********************************************************/ +/****** Info generation for quads w rad parameters ******/ +/********************************************************/ + + +int intxx_radscal_info(intxx_quad_info_type* p, + int (*gparm)(intxx_quad_type*), + int (*dparm)(intxx_quad_type*), + int (*gquad)(intxx_quad_type*), + int (*dquad)(intxx_quad_type*)) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + static const char* names[] = {"RAD_SCAL"}; + static const char* desc[] = {"Radial Scaling Factor"}; + p->ext_params = { + 1, names, desc, + NULL, /* set_name */ + NULL, /* get_name */ + gparm, /* generate */ + dparm /* destroy */ + }; + + p->generate = gquad; + p->destroy = dquad; + + return INTXX_SUCCESS; } -int intxx_destroy_uniform(intxx_quad_type* p) { - return intxx_destroy_impl(p); + +#define INTXX_RADSCAL_GET_INFO_IMPL(cname) \ +int intxx_get_##cname##_quad_info(intxx_quad_info_type* p) { \ + return intxx_radscal_info(p, \ + intxx_generate_##cname##_quad_params, \ + intxx_destroy_##cname##_quad_params, \ + intxx_generate_##cname##_quad, \ + intxx_destroy_##cname##_quad \ + ); \ } +INTXX_RADSCAL_GET_INFO_IMPL(becke) +INTXX_RADSCAL_GET_INFO_IMPL(mhl) +INTXX_RADSCAL_GET_INFO_IMPL(ta) +INTXX_RADSCAL_GET_INFO_IMPL(mk) + +#undef INTXX_RADSCAL_GET_INFO_IMPL + + +/***********************************************/ +/****** Allocate/Deallocate Radial Traits ******/ +/***********************************************/ + +#define INTXX_RAD_TRAITS_IMPL(cname, traits_type) \ +int intxx_generate_##cname##_quad_params(intxx_quad_type* p) { \ + return intxx_generate_params_impl(p); \ +} \ +int intxx_destroy_##cname##_quad_params(intxx_quad_type* p) { \ + return intxx_destroy_params_impl(p); \ +} + +using becke_traits_type = IntegratorXX::BeckeRadialTraits; +using mhl_traits_type = IntegratorXX::MurrayHandyLamingRadialTraits<2>; +using ta_traits_type = IntegratorXX::TreutlerAhlrichsRadialTraits; +using mk_traits_type = IntegratorXX::MuraKnowlesRadialTraits; -#define INTXX_GD_BASIC_IMPL(cname, qname) \ +INTXX_RAD_TRAITS_IMPL(becke, becke_traits_type); +INTXX_RAD_TRAITS_IMPL(mhl, mhl_traits_type ); +INTXX_RAD_TRAITS_IMPL(ta, ta_traits_type ); +INTXX_RAD_TRAITS_IMPL(mk, mk_traits_type ); + +#undef INTXX_RAD_TRAITS_IMPL + + + + +#define INTXX_GD_BASIC_IMPL(cname, qname, ...) \ int intxx_generate_##cname(intxx_quad_type* p) { \ - return intxx_generate_impl(p); \ -} \ + return intxx_generate_noparam_impl(p, ##__VA_ARGS__); \ +} \ int intxx_destroy_##cname(intxx_quad_type* p) { \ - return intxx_destroy_impl(p); \ + return intxx_destroy_impl(p); \ } -using gaussleg_quad_type = IntegratorXX::GaussLegendre; -using gausslob_quad_type = IntegratorXX::GaussLobatto; -using gausscheb1_quad_type = IntegratorXX::GaussChebyshev1; -using gausscheb2_quad_type = IntegratorXX::GaussChebyshev2; -using gausscheb3_quad_type = IntegratorXX::GaussChebyshev3; -using gausscheb2mod_quad_type = IntegratorXX::GaussChebyshev2Modified; +using uniform_quad_type = IntegratorXX::UniformTrapezoid; +using gaussleg_quad_type = IntegratorXX::GaussLegendre; +using gausslob_quad_type = IntegratorXX::GaussLobatto; +using gausscheb1_quad_type = IntegratorXX::GaussChebyshev1; +using gausscheb2_quad_type = IntegratorXX::GaussChebyshev2; +using gausscheb3_quad_type = IntegratorXX::GaussChebyshev3; +using gausscheb2m_quad_type = IntegratorXX::GaussChebyshev2Modified; +INTXX_GD_BASIC_IMPL(uniform, uniform_quad_type, 0.0, 1.0); INTXX_GD_BASIC_IMPL(gaussleg, gaussleg_quad_type ); INTXX_GD_BASIC_IMPL(gausslob, gausslob_quad_type ); INTXX_GD_BASIC_IMPL(gausscheb1, gausscheb1_quad_type ); INTXX_GD_BASIC_IMPL(gausscheb2, gausscheb2_quad_type ); INTXX_GD_BASIC_IMPL(gausscheb3, gausscheb3_quad_type ); -INTXX_GD_BASIC_IMPL(gausscheb2mod, gausscheb2mod_quad_type); +INTXX_GD_BASIC_IMPL(gausscheb2m, gausscheb2m_quad_type); #undef INTXX_GD_BASIC_IMPL +#define INTXX_GD_RAD_IMPL(cname, qname, tname) \ +int intxx_generate_##cname##_quad(intxx_quad_type* p) { \ + return intxx_generate_param_impl(p); \ +}\ +int intxx_destroy_##cname##_quad(intxx_quad_type* p) { \ + return intxx_destroy_impl(p); \ +} + +using becke_quad_type = IntegratorXX::Becke; +using mhl_quad_type = IntegratorXX::MurrayHandyLaming; +using ta_quad_type = IntegratorXX::TreutlerAhlrichs; +using mk_quad_type = IntegratorXX::MuraKnowles; + +INTXX_GD_RAD_IMPL(becke, becke_quad_type, becke_traits_type); +INTXX_GD_RAD_IMPL(mhl, mhl_quad_type, mhl_traits_type); +INTXX_GD_RAD_IMPL(ta, ta_quad_type, ta_traits_type); +INTXX_GD_RAD_IMPL(mk, mk_quad_type, mk_traits_type); + } diff --git a/test/c_api.cxx b/test/c_api.cxx index 971d617..92d3463 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -80,7 +80,7 @@ TEST_CASE("C API") { // Meta data REQUIRE(quad.info != NULL); REQUIRE(quad.npoints == -1); - REQUIRE(quad._state == NULL); + REQUIRE(quad._state_quad == NULL); REQUIRE(quad.info->ext_params.n == 0); REQUIRE(!strcmp(quad.info->name, name)); @@ -103,10 +103,10 @@ TEST_CASE("C API") { // Check Quadrature Generation and Destruction if(base_quad) { intxx_generate_quad(&quad); - REQUIRE(quad._state != NULL); + REQUIRE(quad._state_quad != NULL); // Check validity of the state - auto state_as_quad = reinterpret_cast(quad._state); + auto state_as_quad = reinterpret_cast(quad._state_quad); REQUIRE(state_as_quad->npts() == npts); for(auto i = 0; i < npts; ++ i) { REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad->points()[i], 1e-15)); @@ -114,7 +114,7 @@ TEST_CASE("C API") { } intxx_destroy_quad(&quad); - REQUIRE(quad._state == NULL); + REQUIRE(quad._state_quad == NULL); } } @@ -126,31 +126,31 @@ TEST_CASE("C API") { std::unique_ptr base_quad_default = nullptr; SECTION("Becke") { - //using quad_type = UniformTrapezoid; + using quad_type = Becke; error = intxx_quad_init(&quad, INTXX_RADQ_BECKE); name = "BECKE"; - //base_quad = std::make_unique(base_npts); + base_quad_default = std::make_unique(base_npts); } SECTION("MHL") { - //using quad_type = UniformTrapezoid; + using quad_type = MurrayHandyLaming; error = intxx_quad_init(&quad, INTXX_RADQ_MHL); name = "MURRAY_HANDY_LAMING"; - //base_quad = std::make_unique(base_npts); + base_quad_default = std::make_unique(base_npts); } SECTION("TA") { - //using quad_type = UniformTrapezoid; + using quad_type = TreutlerAhlrichs; error = intxx_quad_init(&quad, INTXX_RADQ_TA); name = "TREUTLER_AHLRICHS"; - //base_quad = std::make_unique(base_npts); + base_quad_default = std::make_unique(base_npts); } SECTION("MK") { - //using quad_type = UniformTrapezoid; + using quad_type = MuraKnowles; error = intxx_quad_init(&quad, INTXX_RADQ_MK); name = "MURA_KNOWLES"; - //base_quad = std::make_unique(base_npts); + base_quad_default = std::make_unique(base_npts); } REQUIRE(error == INTXX_SUCCESS); @@ -158,9 +158,13 @@ TEST_CASE("C API") { // Meta data REQUIRE(quad.info != NULL); REQUIRE(quad.npoints == -1); - REQUIRE(quad._state == NULL); + REQUIRE(quad._state_quad == NULL); REQUIRE(!strcmp(quad.info->name, name)); + REQUIRE(quad.info->ext_params.n == 1); + REQUIRE(!strcmp(quad.info->ext_params.names[0], "RAD_SCAL")); + REQUIRE(!strcmp(quad.info->ext_params.descriptions[0], "Radial Scaling Factor")); + // Get before set int npts; error = intxx_quad_get_npts(&quad, &npts); @@ -176,6 +180,23 @@ TEST_CASE("C API") { error = intxx_quad_get_npts(&quad, &npts); REQUIRE(error == INTXX_SUCCESS); REQUIRE(npts == base_npts); + + // Check Quadrature Generation and Destruction + if(base_quad_default) { + intxx_generate_quad(&quad); + REQUIRE(quad._state_quad != NULL); + + // Check validity of the state + auto state_as_quad = reinterpret_cast(quad._state_quad); + REQUIRE(state_as_quad->npts() == npts); + for(auto i = 0; i < npts; ++ i) { + REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad_default->points()[i], 1e-15)); + REQUIRE_THAT(state_as_quad->weights()[i], Matchers::WithinAbs(name, base_quad_default->weights()[i], 1e-15)); + } + + intxx_destroy_quad(&quad); + REQUIRE(quad._state_quad == NULL); + } } intxx_quad_end(&quad); From 253ae5c5c8b034c5af05abb3b3e64a85cd771f42 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 13:54:00 -0700 Subject: [PATCH 10/22] Added C_API tests to CTest --- test/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7860108..e34b50e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -54,3 +54,5 @@ add_test( NAME QUADRATURES_LEGENDRE COMMAND gausslegendre ) add_test( NAME QUADRATURES_LOBATTO COMMAND gausslobatto ) add_test( NAME QUADRATURES_CHEBYSHEV1 COMMAND gausschebyshev1 ) add_test( NAME QUADRATURES_CHEBYSHEV2 COMMAND gausschebyshev2 ) +add_test( NAME C_API COMMAND c_api ) +add_test( NAME C_API_MEMORY COMMAND valgrind --leak-check=full --error-exitcode=100 $ ) From 8da6d79016b0a27e2d861e017cf51a5dc2fde3f7 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 15:03:38 -0700 Subject: [PATCH 11/22] [CI] Install Valgrind in CI container to enable mem tests for C API --- .github/workflows/build_and_test_compiler_zoo.yml | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/.github/workflows/build_and_test_compiler_zoo.yml b/.github/workflows/build_and_test_compiler_zoo.yml index 4192e5a..a15de2a 100644 --- a/.github/workflows/build_and_test_compiler_zoo.yml +++ b/.github/workflows/build_and_test_compiler_zoo.yml @@ -40,6 +40,10 @@ jobs: shell: bash run: cmake --build ${{runner.workspace}}/build -j2 + - name: Install Valgrind + shell: bash + run: apt install valgrind + - name: Test shell: bash run: cmake --build ${{runner.workspace}}/build --target test @@ -72,6 +76,10 @@ jobs: shell: bash run: cmake --build ${{runner.workspace}}/build -j2 + - name: Install Valgrind + shell: bash + run: apt install valgrind + - name: Test shell: bash run: cmake --build ${{runner.workspace}}/build --target test From 092e19d18ce18341ab11ab7d6ec8243cf0f46ff3 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 15:57:17 -0700 Subject: [PATCH 12/22] Added a C API getter for named parameters, added path for radial quadrature RSCAL getter + UT --- include/integratorxx/c_api.h | 14 +++-- .../integratorxx/quadratures/radial/becke.hpp | 3 + .../integratorxx/quadratures/radial/mhl.hpp | 2 + .../quadratures/radial/muraknowles.hpp | 2 + .../quadratures/radial/treutlerahlrichs.hpp | 2 + src/c_api/c_api.c | 15 +++++ src/c_api/c_radial.cxx | 58 +++++++++++++++++-- test/c_api.cxx | 15 +++++ 8 files changed, 103 insertions(+), 8 deletions(-) diff --git a/include/integratorxx/c_api.h b/include/integratorxx/c_api.h index a32b3a4..df1cab9 100644 --- a/include/integratorxx/c_api.h +++ b/include/integratorxx/c_api.h @@ -7,7 +7,7 @@ #define INTXX_NULL_INFOPTR -3 #define INTXX_INVALID_OPT -4 #define INTXX_INVALID_ARG -5 -#define INTXX_INVALID_OUT -5 +#define INTXX_INVALID_OUT -6 /*** Quadrature Classes ***/ #define INTXX_PRM_QUAD 1 @@ -52,9 +52,9 @@ typedef struct { const char** names; ///< Names of the params const char** descriptions; ///< Long descriptions of params - void (*set_name)(struct intxx_quad_type* p, const char* name, double v); - ///< Set parameter function - void (*get_name)(struct intxx_quad_type* p, const char* name, double* v); + 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); @@ -145,6 +145,12 @@ 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 diff --git a/include/integratorxx/quadratures/radial/becke.hpp b/include/integratorxx/quadratures/radial/becke.hpp index 1c711ae..7c56fa0 100644 --- a/include/integratorxx/quadratures/radial/becke.hpp +++ b/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_; } }; /** diff --git a/include/integratorxx/quadratures/radial/mhl.hpp b/include/integratorxx/quadratures/radial/mhl.hpp index a45fe44..91565d8 100644 --- a/include/integratorxx/quadratures/radial/mhl.hpp +++ b/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_; } }; diff --git a/include/integratorxx/quadratures/radial/muraknowles.hpp b/include/integratorxx/quadratures/radial/muraknowles.hpp index 23a26ff..cf00f0b 100644 --- a/include/integratorxx/quadratures/radial/muraknowles.hpp +++ b/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 diff --git a/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp b/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp index 9aa5e96..46c41a6 100644 --- a/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp +++ b/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_; } }; diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c index 00e607d..e7732cd 100644 --- a/src/c_api/c_api.c +++ b/src/c_api/c_api.c @@ -120,3 +120,18 @@ int intxx_destroy_quad(intxx_quad_type* p) { return p->info->destroy(p); else return INTXX_SUCCESS; } + +int intxx_get_ext_param_name(intxx_quad_type* p, const char* name, double* v) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + if(p->info->ext_params.n <= 0) { + return INTXX_INVALID_OPT; + } + + if(p->info->ext_params.get_name == NULL) { + // Something more descriptive? + return INTXX_INVALID_OPT; + } + + return p->info->ext_params.get_name(p, name, v); +} diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index 3df70b3..6b5c282 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -3,6 +3,7 @@ #include "c_internal.h" #include +#include template auto generate_quadrature(Args&&... args) { @@ -121,6 +122,39 @@ int intxx_destroy_params_impl(intxx_quad_type* p) { return INTXX_SUCCESS; } +template +int intxx_get_rad_scal_impl(intxx_quad_type* p, double* R) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(R == NULL) { + // Return error code + } + + auto ext_param = p->info->ext_params; + + // Check that this quadrature has a radial scaling factor + bool r_found = false; + for(int i = 0; i < ext_param.n; ++i) { + if(r_found) break; + const auto name = ext_param.names[i]; + r_found = strcmp(name, "RAD_SCAL"); + } + + if(not r_found) { + // Return error code + } + + if(p->_state_parm == NULL) { + // Return error code + } + + // Read data + *R = reinterpret_cast(p->_state_parm)->R(); + + return INTXX_SUCCESS; +} + extern "C" { @@ -144,7 +178,8 @@ int intxx_get_##cname##_quad_info(intxx_quad_info_type* p); \ int intxx_generate_##cname##_quad(intxx_quad_type* p); \ int intxx_destroy_##cname##_quad(intxx_quad_type* p); \ int intxx_generate_##cname##_quad_params(intxx_quad_type* p); \ -int intxx_destroy_##cname##_quad_params(intxx_quad_type* p); +int intxx_destroy_##cname##_quad_params(intxx_quad_type* p); \ +int intxx_get_name_##cname##_quad(intxx_quad_type* p, const char* name, double *v); FWD_RAD(becke) FWD_RAD(mhl) @@ -262,6 +297,8 @@ INTXX_NOPARAM_GET_INFO_IMPL(gausscheb3); int intxx_radscal_info(intxx_quad_info_type* p, + int (*sn)(intxx_quad_type*,const char*,double), + int (*gn)(intxx_quad_type*,const char*,double*), int (*gparm)(intxx_quad_type*), int (*dparm)(intxx_quad_type*), int (*gquad)(intxx_quad_type*), @@ -272,8 +309,8 @@ int intxx_radscal_info(intxx_quad_info_type* p, static const char* desc[] = {"Radial Scaling Factor"}; p->ext_params = { 1, names, desc, - NULL, /* set_name */ - NULL, /* get_name */ + sn, /* set_name */ + gn, /* get_name */ gparm, /* generate */ dparm /* destroy */ }; @@ -287,6 +324,8 @@ int intxx_radscal_info(intxx_quad_info_type* p, #define INTXX_RADSCAL_GET_INFO_IMPL(cname) \ int intxx_get_##cname##_quad_info(intxx_quad_info_type* p) { \ return intxx_radscal_info(p, \ + NULL, \ + intxx_get_name_##cname##_quad, \ intxx_generate_##cname##_quad_params, \ intxx_destroy_##cname##_quad_params, \ intxx_generate_##cname##_quad, \ @@ -326,8 +365,19 @@ INTXX_RAD_TRAITS_IMPL(mk, mk_traits_type ); #undef INTXX_RAD_TRAITS_IMPL +/**********************************/ +/****** Radial Quad GET_NAME ******/ +/**********************************/ +#define INTXX_RAD_GETSET_IMPL(cname, traits_type) \ +int intxx_get_name_##cname##_quad(intxx_quad_type* p, const char* name, double *v) { \ + if(strcmp(name, "RAD_SCAL")){ return INTXX_INVALID_OPT; } \ + return intxx_get_rad_scal_impl(p, v);\ +} - +INTXX_RAD_GETSET_IMPL(becke, becke_traits_type); +INTXX_RAD_GETSET_IMPL(mhl, mhl_traits_type ); +INTXX_RAD_GETSET_IMPL(ta, ta_traits_type ); +INTXX_RAD_GETSET_IMPL(mk, mk_traits_type ); #define INTXX_GD_BASIC_IMPL(cname, qname, ...) \ int intxx_generate_##cname(intxx_quad_type* p) { \ diff --git a/test/c_api.cxx b/test/c_api.cxx index 92d3463..557d07b 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -122,14 +122,18 @@ TEST_CASE("C API") { SECTION("Radial") { const char* name; const int base_npts = 100; + const double RSCAL = 2.0; using base_quad_type = QuadratureBase, std::vector>; std::unique_ptr base_quad_default = nullptr; + std::unique_ptr base_quad_scaled = nullptr; SECTION("Becke") { using quad_type = Becke; + using traits_type = typename quad_type::traits_type; error = intxx_quad_init(&quad, INTXX_RADQ_BECKE); name = "BECKE"; base_quad_default = std::make_unique(base_npts); + base_quad_scaled = std::make_unique(base_npts,traits_type(RSCAL)); } SECTION("MHL") { @@ -161,10 +165,16 @@ TEST_CASE("C API") { REQUIRE(quad._state_quad == NULL); REQUIRE(!strcmp(quad.info->name, name)); + + // Check default parameters REQUIRE(quad.info->ext_params.n == 1); REQUIRE(!strcmp(quad.info->ext_params.names[0], "RAD_SCAL")); REQUIRE(!strcmp(quad.info->ext_params.descriptions[0], "Radial Scaling Factor")); + double R; + intxx_get_ext_param_name(&quad, "RAD_SCAL", &R); + REQUIRE(R == 1.0); + // Get before set int npts; error = intxx_quad_get_npts(&quad, &npts); @@ -197,6 +207,11 @@ TEST_CASE("C API") { intxx_destroy_quad(&quad); REQUIRE(quad._state_quad == NULL); } + + // Check Scaling Works + if(base_quad_scaled) { + + } } intxx_quad_end(&quad); From 7c377ee47028ea6d0b4d3e970a685dc9b710c403 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 16:31:29 -0700 Subject: [PATCH 13/22] Added ext_param setters for C API, tested paths for radial quadratures --- src/c_api/c_api.c | 15 +++++++++++++ src/c_api/c_radial.cxx | 48 ++++++++++++++++++++++++++++++++++++++++-- test/c_api.cxx | 26 +++++++++++++++++++++++ 3 files changed, 87 insertions(+), 2 deletions(-) diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c index e7732cd..f6ea373 100644 --- a/src/c_api/c_api.c +++ b/src/c_api/c_api.c @@ -135,3 +135,18 @@ int intxx_get_ext_param_name(intxx_quad_type* p, const char* name, double* v) { return p->info->ext_params.get_name(p, name, v); } + +int intxx_set_ext_param_name(intxx_quad_type* p, const char* name, double v) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + if(p->info->ext_params.n <= 0) { + return INTXX_INVALID_OPT; + } + + if(p->info->ext_params.set_name == NULL) { + // Something more descriptive? + return INTXX_INVALID_OPT; + } + + return p->info->ext_params.set_name(p, name, v); +} diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index 6b5c282..0fce7d0 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -155,6 +155,45 @@ int intxx_get_rad_scal_impl(intxx_quad_type* p, double* R) { return INTXX_SUCCESS; } +template +int intxx_set_rad_scal_impl(intxx_quad_type* p, double R) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(R <= 0.0) { + // Return error code + } + + auto ext_param = p->info->ext_params; + + // Check that this quadrature has a radial scaling factor + bool r_found = false; + for(int i = 0; i < ext_param.n; ++i) { + if(r_found) break; + const auto name = ext_param.names[i]; + r_found = strcmp(name, "RAD_SCAL"); + } + + if(not r_found) { + // Return error code + } + + if(p->_state_parm == NULL) { + // Return error code + } + + // Overwrite + *reinterpret_cast(p->_state_parm) = TraitsType(R); + + // Regenerate if quadrature is populated + if(p->_state_quad) { + intxx_destroy_quad(p); + intxx_generate_quad(p); + } + + return INTXX_SUCCESS; +} + extern "C" { @@ -179,7 +218,8 @@ int intxx_generate_##cname##_quad(intxx_quad_type* p); \ int intxx_destroy_##cname##_quad(intxx_quad_type* p); \ int intxx_generate_##cname##_quad_params(intxx_quad_type* p); \ int intxx_destroy_##cname##_quad_params(intxx_quad_type* p); \ -int intxx_get_name_##cname##_quad(intxx_quad_type* p, const char* name, double *v); +int intxx_get_name_##cname##_quad(intxx_quad_type* p, const char* name, double *v); \ +int intxx_set_name_##cname##_quad(intxx_quad_type* p, const char* name, double v); FWD_RAD(becke) FWD_RAD(mhl) @@ -324,7 +364,7 @@ int intxx_radscal_info(intxx_quad_info_type* p, #define INTXX_RADSCAL_GET_INFO_IMPL(cname) \ int intxx_get_##cname##_quad_info(intxx_quad_info_type* p) { \ return intxx_radscal_info(p, \ - NULL, \ + intxx_set_name_##cname##_quad, \ intxx_get_name_##cname##_quad, \ intxx_generate_##cname##_quad_params, \ intxx_destroy_##cname##_quad_params, \ @@ -372,6 +412,10 @@ INTXX_RAD_TRAITS_IMPL(mk, mk_traits_type ); int intxx_get_name_##cname##_quad(intxx_quad_type* p, const char* name, double *v) { \ if(strcmp(name, "RAD_SCAL")){ return INTXX_INVALID_OPT; } \ return intxx_get_rad_scal_impl(p, v);\ +}\ +int intxx_set_name_##cname##_quad(intxx_quad_type* p, const char* name, double v) { \ + if(strcmp(name, "RAD_SCAL")){ return INTXX_INVALID_OPT; } \ + return intxx_set_rad_scal_impl(p, v);\ } INTXX_RAD_GETSET_IMPL(becke, becke_traits_type); diff --git a/test/c_api.cxx b/test/c_api.cxx index 557d07b..deea65d 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -210,7 +210,33 @@ TEST_CASE("C API") { // Check Scaling Works if(base_quad_scaled) { + intxx_set_ext_param_name(&quad, "RAD_SCAL", RSCAL); + // Check change stuck + intxx_get_ext_param_name(&quad, "RAD_SCAL", &R); + REQUIRE(R == RSCAL); + + // Regenerate and check + intxx_generate_quad(&quad); + + auto state_as_quad = reinterpret_cast(quad._state_quad); + REQUIRE(state_as_quad->npts() == npts); + for(auto i = 0; i < npts; ++ i) { + REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad_scaled->points()[i], 1e-15)); + REQUIRE_THAT(state_as_quad->weights()[i], Matchers::WithinAbs(name, base_quad_scaled->weights()[i], 1e-15)); + } + + + // Set without destroying to make sure quadratures are regerated + intxx_set_ext_param_name(&quad, "RAD_SCAL", 1.0); + state_as_quad = reinterpret_cast(quad._state_quad); + REQUIRE(state_as_quad->npts() == npts); + for(auto i = 0; i < npts; ++ i) { + REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad_default->points()[i], 1e-15)); + REQUIRE_THAT(state_as_quad->weights()[i], Matchers::WithinAbs(name, base_quad_default->weights()[i], 1e-15)); + } + + } } From e1cf8df69fade401b1ff863ed410a08568318d36 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 16:32:44 -0700 Subject: [PATCH 14/22] Ammend 8da6d79016b0a27e2d861e017cf51a5dc2fde3f7 --- .github/workflows/build_and_test_compiler_zoo.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build_and_test_compiler_zoo.yml b/.github/workflows/build_and_test_compiler_zoo.yml index a15de2a..0db7a6e 100644 --- a/.github/workflows/build_and_test_compiler_zoo.yml +++ b/.github/workflows/build_and_test_compiler_zoo.yml @@ -42,7 +42,7 @@ jobs: - name: Install Valgrind shell: bash - run: apt install valgrind + run: apt install -y valgrind - name: Test shell: bash @@ -78,7 +78,7 @@ jobs: - name: Install Valgrind shell: bash - run: apt install valgrind + run: apt install -y valgrind - name: Test shell: bash From faf09f40973f5575f588e7837c1a0ef53ff09f71 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Sun, 24 Sep 2023 16:42:03 -0700 Subject: [PATCH 15/22] Ammend e1cf8df69fade401b1ff863ed410a08568318d36 --- .github/workflows/build_and_test_compiler_zoo.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build_and_test_compiler_zoo.yml b/.github/workflows/build_and_test_compiler_zoo.yml index 0db7a6e..e5e4db8 100644 --- a/.github/workflows/build_and_test_compiler_zoo.yml +++ b/.github/workflows/build_and_test_compiler_zoo.yml @@ -42,7 +42,7 @@ jobs: - name: Install Valgrind shell: bash - run: apt install -y valgrind + run: apt update && apt install -y valgrind - name: Test shell: bash @@ -78,7 +78,7 @@ jobs: - name: Install Valgrind shell: bash - run: apt install -y valgrind + run: apt update && apt install -y valgrind - name: Test shell: bash From 3d31ea4d87389aee78b9ffd83cd432d16fc34105 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 25 Sep 2023 09:06:33 -0700 Subject: [PATCH 16/22] Make C_API_MEMORY test contingent on valgrind discovery, disable in CI --- .github/workflows/build_and_test_compiler_zoo.yml | 8 -------- test/CMakeLists.txt | 9 ++++++++- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/.github/workflows/build_and_test_compiler_zoo.yml b/.github/workflows/build_and_test_compiler_zoo.yml index e5e4db8..4192e5a 100644 --- a/.github/workflows/build_and_test_compiler_zoo.yml +++ b/.github/workflows/build_and_test_compiler_zoo.yml @@ -40,10 +40,6 @@ jobs: shell: bash run: cmake --build ${{runner.workspace}}/build -j2 - - name: Install Valgrind - shell: bash - run: apt update && apt install -y valgrind - - name: Test shell: bash run: cmake --build ${{runner.workspace}}/build --target test @@ -76,10 +72,6 @@ jobs: shell: bash run: cmake --build ${{runner.workspace}}/build -j2 - - name: Install Valgrind - shell: bash - run: apt update && apt install -y valgrind - - name: Test shell: bash run: cmake --build ${{runner.workspace}}/build --target test diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e34b50e..ea59a35 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -55,4 +55,11 @@ add_test( NAME QUADRATURES_LOBATTO COMMAND gausslobatto ) add_test( NAME QUADRATURES_CHEBYSHEV1 COMMAND gausschebyshev1 ) add_test( NAME QUADRATURES_CHEBYSHEV2 COMMAND gausschebyshev2 ) add_test( NAME C_API COMMAND c_api ) -add_test( NAME C_API_MEMORY COMMAND valgrind --leak-check=full --error-exitcode=100 $ ) + +# Memory check +find_program( VALGRIND_EXECUTABLE valgrind ) +if(VALGRIND_EXECUTABLE) + message(STATUS "Found Valgrind: ${VALGRIND_EXECUTABLE}") + message(STATUS "Enabling C_API_MEMORY test") + add_test( NAME C_API_MEMORY COMMAND ${VALGRIND_EXECUTABLE} --leak-check=full --error-exitcode=100 $ ) +endif() From b326aefece8e87b33d9eea9496bb90150e09df98 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 25 Sep 2023 10:06:50 -0700 Subject: [PATCH 17/22] Refactor C API source --- src/c_api/CMakeLists.txt | 2 +- src/c_api/c_api_util.hpp | 195 ++++++++++++++++++++ src/c_api/c_internal.h | 13 ++ src/c_api/c_misc.cxx | 59 ++++++ src/c_api/c_primitive.cxx | 115 ++++++++++++ src/c_api/c_radial.cxx | 372 +++----------------------------------- 6 files changed, 410 insertions(+), 346 deletions(-) create mode 100644 src/c_api/c_api_util.hpp create mode 100644 src/c_api/c_misc.cxx create mode 100644 src/c_api/c_primitive.cxx diff --git a/src/c_api/CMakeLists.txt b/src/c_api/CMakeLists.txt index 12616f8..16401c1 100644 --- a/src/c_api/CMakeLists.txt +++ b/src/c_api/CMakeLists.txt @@ -1 +1 @@ -add_library( integratorxx c_api.c c_radial.cxx ) +add_library( integratorxx c_api.c c_primitive.cxx c_radial.cxx c_misc.cxx ) diff --git a/src/c_api/c_api_util.hpp b/src/c_api/c_api_util.hpp new file mode 100644 index 0000000..e736cfb --- /dev/null +++ b/src/c_api/c_api_util.hpp @@ -0,0 +1,195 @@ +#pragma once +#include "c_internal.h" +#include +#include + +template +auto generate_quadrature(Args&&... args) { + using alloc_type = std::allocator; + using alloc_traits = std::allocator_traits; + + alloc_type alloc; + auto ptr = alloc_traits::allocate(alloc, 1); + alloc_traits::construct(alloc, ptr, std::forward(args)...); + + return ptr; +} + +template +int intxx_generate_noparam_impl(intxx_quad_type* p, Args&&... args) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + using quad_type = QuadType; + + int npts = p->npoints; + if( npts <= 0 ) { + // Return error code + } + + if(p->_state_quad != NULL) { + // Check if we need to regenerate + } + + // Allocate and construct the quadrature instance + auto ptr = generate_quadrature(npts, std::forward(args)... ); + + // Store the pointer + p->_state_quad = (void*)ptr; + + return INTXX_SUCCESS; + +} + + +template +int intxx_generate_param_impl(intxx_quad_type* p, Args&&... args) { + if(p == NULL) return INTXX_NULL_QUADPTR; + + if(p->_state_parm == NULL) { + // Return error code + } + + // Cross your fingers... + const auto _params = reinterpret_cast(p->_state_parm); + return intxx_generate_noparam_impl(p, *_params, std::forward(args)... ); + +} + +template +void destroy_quadrature(void* ptr) { + using alloc_type = std::allocator; + using alloc_traits = std::allocator_traits; + + alloc_type alloc; + + // Cross your fingers... + auto quad_ptr = reinterpret_cast(ptr); + + // Destroy and deallocate + alloc_traits::destroy(alloc, quad_ptr); + alloc_traits::deallocate(alloc, quad_ptr, 1); +} + +template +int intxx_destroy_impl(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + + using quad_type = QuadType; + + if(p->_state_quad != NULL) { + // Destroy quadrature + destroy_quadrature(p->_state_quad); + + // Null out state + p->_state_quad = NULL; + } + return INTXX_SUCCESS; +} + +template +int intxx_generate_params_impl(intxx_quad_type* p, Args&&... args) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(p->_state_parm != NULL) { + return INTXX_SUCCESS; + } + + // Allocate and construct the quadrature instance + auto ptr = generate_quadrature(std::forward(args)... ); + + // Store the pointer + p->_state_parm = (void*)ptr; + + return INTXX_SUCCESS; + +} + +template +int intxx_destroy_params_impl(intxx_quad_type* p) { + if(p == NULL) return INTXX_NULL_QUADPTR; + + if(p->_state_parm != NULL) { + // Destroy traits object + destroy_quadrature(p->_state_parm); + + // Null out state + p->_state_parm = NULL; + } + return INTXX_SUCCESS; +} + +template +int intxx_get_rad_scal_impl(intxx_quad_type* p, double* R) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(R == NULL) { + // Return error code + } + + auto ext_param = p->info->ext_params; + + // Check that this quadrature has a radial scaling factor + bool r_found = false; + for(int i = 0; i < ext_param.n; ++i) { + if(r_found) break; + const auto name = ext_param.names[i]; + r_found = strcmp(name, "RAD_SCAL"); + } + + if(not r_found) { + // Return error code + } + + if(p->_state_parm == NULL) { + // Return error code + } + + // Read data + *R = reinterpret_cast(p->_state_parm)->R(); + + return INTXX_SUCCESS; +} + +template +int intxx_set_rad_scal_impl(intxx_quad_type* p, double R) { + if(p == NULL) return INTXX_NULL_QUADPTR; + if(p->info == NULL) return INTXX_NULL_INFOPTR; + + if(R <= 0.0) { + // Return error code + } + + auto ext_param = p->info->ext_params; + + // Check that this quadrature has a radial scaling factor + bool r_found = false; + for(int i = 0; i < ext_param.n; ++i) { + if(r_found) break; + const auto name = ext_param.names[i]; + r_found = strcmp(name, "RAD_SCAL"); + } + + if(not r_found) { + // Return error code + } + + if(p->_state_parm == NULL) { + // Return error code + } + + // Overwrite + *reinterpret_cast(p->_state_parm) = TraitsType(R); + + // Regenerate if quadrature is populated + if(p->_state_quad) { + intxx_destroy_quad(p); + intxx_generate_quad(p); + } + + return INTXX_SUCCESS; +} + + diff --git a/src/c_api/c_internal.h b/src/c_api/c_internal.h index 1ead231..b7b62a8 100644 --- a/src/c_api/c_internal.h +++ b/src/c_api/c_internal.h @@ -10,6 +10,19 @@ int intxx_get_radq_info(intxx_quad_info_type* p, int quad); int intxx_get_angq_info(intxx_quad_info_type* p, int quad); int intxx_get_prmq_info(intxx_quad_info_type* p, int quad); + +int intxx_noparam_info(intxx_quad_info_type* p, + int (*g)(intxx_quad_type*), + int (*d)(intxx_quad_type*)); + +int intxx_radscal_info(intxx_quad_info_type* p, + int (*sn)(intxx_quad_type*,const char*,double), + int (*gn)(intxx_quad_type*,const char*,double*), + int (*gparm)(intxx_quad_type*), + int (*dparm)(intxx_quad_type*), + int (*gquad)(intxx_quad_type*), + int (*dquad)(intxx_quad_type*)); + #ifdef __cplusplus } #endif diff --git a/src/c_api/c_misc.cxx b/src/c_api/c_misc.cxx new file mode 100644 index 0000000..639b797 --- /dev/null +++ b/src/c_api/c_misc.cxx @@ -0,0 +1,59 @@ +#include "c_internal.h" +#include + +extern "C" { + +void intxx_default_quad_info(intxx_quad_info_type* p) { + p->number = 0; // Invalid + p->kind = 0; // Invalid + p->name = "Default"; // Default state + p->generate = NULL; // Invalid + p->destroy = NULL; // Invalid +} + +/******************************************************/ +/****** Info generation for quads w/o parameters ******/ +/******************************************************/ + +int intxx_noparam_info(intxx_quad_info_type* p, + int (*g)(intxx_quad_type*), + int (*d)(intxx_quad_type*)) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + p->ext_params.n = 0; /// No External Parameters + p->generate = g; + p->destroy = d; + + return INTXX_SUCCESS; +} + +/********************************************************/ +/****** Info generation for quads w rad parameters ******/ +/********************************************************/ + +int intxx_radscal_info(intxx_quad_info_type* p, + int (*sn)(intxx_quad_type*,const char*,double), + int (*gn)(intxx_quad_type*,const char*,double*), + int (*gparm)(intxx_quad_type*), + int (*dparm)(intxx_quad_type*), + int (*gquad)(intxx_quad_type*), + int (*dquad)(intxx_quad_type*)) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + static const char* names[] = {"RAD_SCAL"}; + static const char* desc[] = {"Radial Scaling Factor"}; + p->ext_params = { + 1, names, desc, + sn, /* set_name */ + gn, /* get_name */ + gparm, /* generate */ + dparm /* destroy */ + }; + + p->generate = gquad; + p->destroy = dquad; + + return INTXX_SUCCESS; +} + +} diff --git a/src/c_api/c_primitive.cxx b/src/c_api/c_primitive.cxx new file mode 100644 index 0000000..21fc568 --- /dev/null +++ b/src/c_api/c_primitive.cxx @@ -0,0 +1,115 @@ +#include + +#include "c_internal.h" +#include "c_api_util.hpp" + +extern "C" { + +// C++ types for internal quadrature states +using uniform_quad_type = IntegratorXX::UniformTrapezoid; +using gaussleg_quad_type = IntegratorXX::GaussLegendre; +using gausslob_quad_type = IntegratorXX::GaussLobatto; +using gausscheb1_quad_type = IntegratorXX::GaussChebyshev1; +using gausscheb2_quad_type = IntegratorXX::GaussChebyshev2; +using gausscheb3_quad_type = IntegratorXX::GaussChebyshev3; +using gausscheb2m_quad_type = IntegratorXX::GaussChebyshev2Modified; + +/*************************************************************/ +/****** Forward declare C-API for Primitive Quadratures ******/ +/*************************************************************/ + +#define FWD_PRM(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_PRM(uniform) +FWD_PRM(gaussleg) +FWD_PRM(gausslob) +FWD_PRM(gausscheb1) +FWD_PRM(gausscheb2) +FWD_PRM(gausscheb2m) +FWD_PRM(gausscheb3) + +#undef FWD_PRM + +/*************************************************************/ +/****** Runtime Generator for Primitive Quadrature Info ******/ +/*************************************************************/ + +int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { + if(p == NULL) return INTXX_NULL_INFOPTR; + + quad = quad & INTXX_PRMQ_MASK; + p->number = quad; + p->kind = INTXX_PRM_QUAD; + switch(quad) { + case INTXX_PRMQ_UNIFORM: + p->name = "UNIFORM"; + return intxx_get_uniform_info(p); + case INTXX_PRMQ_GAUSSLEG: + p->name = "GAUSS_LEGENDRE"; + return intxx_get_gaussleg_info(p); + case INTXX_PRMQ_GAUSSCHEB_1: + p->name = "GAUSS_CHEBYSHEV_1"; + return intxx_get_gausscheb1_info(p); + case INTXX_PRMQ_GAUSSCHEB_2: + p->name = "GAUSS_CHEBYSHEV_2"; + return intxx_get_gausscheb2_info(p); + case INTXX_PRMQ_GAUSSCHEB_2MOD: + p->name = "GAUSS_CHEBYSHEV_2MOD"; + return intxx_get_gausscheb2m_info(p); + case INTXX_PRMQ_GAUSSCHEB_3: + p->name = "GAUSS_CHEBYSHEV_3"; + return intxx_get_gausscheb3_info(p); + case INTXX_PRMQ_GAUSSLOB: + p->name = "GAUSS_LOBATTO"; + return intxx_get_gausslob_info(p); + default: + return INTXX_INVALID_QUAD; + } +} + +/*******************************************************/ +/****** Info generation for Primitive 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_NOPARAM_GET_INFO_IMPL(uniform); +INTXX_NOPARAM_GET_INFO_IMPL(gaussleg); +INTXX_NOPARAM_GET_INFO_IMPL(gausslob); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb1); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb2); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb2m); +INTXX_NOPARAM_GET_INFO_IMPL(gausscheb3); + +#undef INTXX_NOPARAM_GET_INFO_IMPL + +/*************************************************************/ +/****** Allocate/Deallocate Primitive Quadrature States ******/ +/*************************************************************/ + +#define INTXX_GD_BASIC_IMPL(cname, qname, ...) \ +int intxx_generate_##cname(intxx_quad_type* p) { \ + return intxx_generate_noparam_impl(p, ##__VA_ARGS__); \ +} \ +int intxx_destroy_##cname(intxx_quad_type* p) { \ + return intxx_destroy_impl(p); \ +} + +INTXX_GD_BASIC_IMPL(uniform, uniform_quad_type, 0.0, 1.0); +INTXX_GD_BASIC_IMPL(gaussleg, gaussleg_quad_type ); +INTXX_GD_BASIC_IMPL(gausslob, gausslob_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb1, gausscheb1_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb2, gausscheb2_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb3, gausscheb3_quad_type ); +INTXX_GD_BASIC_IMPL(gausscheb2m, gausscheb2m_quad_type); + +#undef INTXX_GD_BASIC_IMPL + +} // extern C diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index 0fce7d0..076ff28 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -1,216 +1,26 @@ -#include #include #include "c_internal.h" -#include -#include - -template -auto generate_quadrature(Args&&... args) { - using alloc_type = std::allocator; - using alloc_traits = std::allocator_traits; - - alloc_type alloc; - auto ptr = alloc_traits::allocate(alloc, 1); - alloc_traits::construct(alloc, ptr, std::forward(args)...); - - return ptr; -} - -template -int intxx_generate_noparam_impl(intxx_quad_type* p, Args&&... args) { - if(p == NULL) return INTXX_NULL_QUADPTR; - if(p->info == NULL) return INTXX_NULL_INFOPTR; - - using quad_type = QuadType; - - int npts = p->npoints; - if( npts <= 0 ) { - // Return error code - } - - if(p->_state_quad != NULL) { - // Check if we need to regenerate - } - - // Allocate and construct the quadrature instance - auto ptr = generate_quadrature(npts, std::forward(args)... ); - - // Store the pointer - p->_state_quad = (void*)ptr; - - return INTXX_SUCCESS; - -} - - -template -int intxx_generate_param_impl(intxx_quad_type* p, Args&&... args) { - if(p == NULL) return INTXX_NULL_QUADPTR; - - if(p->_state_parm == NULL) { - // Return error code - } - - // Cross your fingers... - const auto _params = reinterpret_cast(p->_state_parm); - return intxx_generate_noparam_impl(p, *_params, std::forward(args)... ); - -} - -template -void destroy_quadrature(void* ptr) { - using alloc_type = std::allocator; - using alloc_traits = std::allocator_traits; - - alloc_type alloc; - - // Cross your fingers... - auto quad_ptr = reinterpret_cast(ptr); - - // Destroy and deallocate - alloc_traits::destroy(alloc, quad_ptr); - alloc_traits::deallocate(alloc, quad_ptr, 1); -} - -template -int intxx_destroy_impl(intxx_quad_type* p) { - if(p == NULL) return INTXX_NULL_QUADPTR; - - using quad_type = QuadType; - - if(p->_state_quad != NULL) { - // Destroy quadrature - destroy_quadrature(p->_state_quad); - - // Null out state - p->_state_quad = NULL; - } - return INTXX_SUCCESS; -} - -template -int intxx_generate_params_impl(intxx_quad_type* p, Args&&... args) { - if(p == NULL) return INTXX_NULL_QUADPTR; - if(p->info == NULL) return INTXX_NULL_INFOPTR; - - if(p->_state_parm != NULL) { - return INTXX_SUCCESS; - } - - // Allocate and construct the quadrature instance - auto ptr = generate_quadrature(std::forward(args)... ); - - // Store the pointer - p->_state_parm = (void*)ptr; - - return INTXX_SUCCESS; - -} - -template -int intxx_destroy_params_impl(intxx_quad_type* p) { - if(p == NULL) return INTXX_NULL_QUADPTR; - - if(p->_state_parm != NULL) { - // Destroy traits object - destroy_quadrature(p->_state_parm); - - // Null out state - p->_state_parm = NULL; - } - return INTXX_SUCCESS; -} - -template -int intxx_get_rad_scal_impl(intxx_quad_type* p, double* R) { - if(p == NULL) return INTXX_NULL_QUADPTR; - if(p->info == NULL) return INTXX_NULL_INFOPTR; - - if(R == NULL) { - // Return error code - } - - auto ext_param = p->info->ext_params; - - // Check that this quadrature has a radial scaling factor - bool r_found = false; - for(int i = 0; i < ext_param.n; ++i) { - if(r_found) break; - const auto name = ext_param.names[i]; - r_found = strcmp(name, "RAD_SCAL"); - } - - if(not r_found) { - // Return error code - } - - if(p->_state_parm == NULL) { - // Return error code - } - - // Read data - *R = reinterpret_cast(p->_state_parm)->R(); - - return INTXX_SUCCESS; -} - -template -int intxx_set_rad_scal_impl(intxx_quad_type* p, double R) { - if(p == NULL) return INTXX_NULL_QUADPTR; - if(p->info == NULL) return INTXX_NULL_INFOPTR; - - if(R <= 0.0) { - // Return error code - } - - auto ext_param = p->info->ext_params; - - // Check that this quadrature has a radial scaling factor - bool r_found = false; - for(int i = 0; i < ext_param.n; ++i) { - if(r_found) break; - const auto name = ext_param.names[i]; - r_found = strcmp(name, "RAD_SCAL"); - } - - if(not r_found) { - // Return error code - } - - if(p->_state_parm == NULL) { - // Return error code - } - - // Overwrite - *reinterpret_cast(p->_state_parm) = TraitsType(R); - - // Regenerate if quadrature is populated - if(p->_state_quad) { - intxx_destroy_quad(p); - intxx_generate_quad(p); - } - - return INTXX_SUCCESS; -} +#include "c_api_util.hpp" extern "C" { -#define FWD_PRM(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_PRM(uniform) -FWD_PRM(gaussleg) -FWD_PRM(gausslob) -FWD_PRM(gausscheb1) -FWD_PRM(gausscheb2) -FWD_PRM(gausscheb2m) -FWD_PRM(gausscheb3) +// C++ types for internal quadrature states +using becke_traits_type = IntegratorXX::BeckeRadialTraits; +using mhl_traits_type = IntegratorXX::MurrayHandyLamingRadialTraits<2>; +using ta_traits_type = IntegratorXX::TreutlerAhlrichsRadialTraits; +using mk_traits_type = IntegratorXX::MuraKnowlesRadialTraits; + +using becke_quad_type = IntegratorXX::Becke; +using mhl_quad_type = IntegratorXX::MurrayHandyLaming; +using ta_quad_type = IntegratorXX::TreutlerAhlrichs; +using mk_quad_type = IntegratorXX::MuraKnowles; -#undef FWD_PRM +/**********************************************************/ +/****** Forward declare C-API for Radial Quadratures ******/ +/**********************************************************/ #define FWD_RAD(cname) \ int intxx_get_##cname##_quad_info(intxx_quad_info_type* p); \ @@ -228,51 +38,9 @@ FWD_RAD(mk) #undef FWD_RAD - - - - -void intxx_default_quad_info(intxx_quad_info_type* p) { - p->number = 0; // Invalid - p->kind = 0; // Invalid - p->name = "Default"; // Default state - p->generate = NULL; // Invalid - p->destroy = NULL; // Invalid -} - -int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { - if(p == NULL) return INTXX_NULL_INFOPTR; - - quad = quad & INTXX_PRMQ_MASK; - p->number = quad; - p->kind = INTXX_PRM_QUAD; - switch(quad) { - case INTXX_PRMQ_UNIFORM: - p->name = "UNIFORM"; - return intxx_get_uniform_info(p); - case INTXX_PRMQ_GAUSSLEG: - p->name = "GAUSS_LEGENDRE"; - return intxx_get_gaussleg_info(p); - case INTXX_PRMQ_GAUSSCHEB_1: - p->name = "GAUSS_CHEBYSHEV_1"; - return intxx_get_gausscheb1_info(p); - case INTXX_PRMQ_GAUSSCHEB_2: - p->name = "GAUSS_CHEBYSHEV_2"; - return intxx_get_gausscheb2_info(p); - case INTXX_PRMQ_GAUSSCHEB_2MOD: - p->name = "GAUSS_CHEBYSHEV_2MOD"; - return intxx_get_gausscheb2m_info(p); - case INTXX_PRMQ_GAUSSCHEB_3: - p->name = "GAUSS_CHEBYSHEV_3"; - return intxx_get_gausscheb3_info(p); - case INTXX_PRMQ_GAUSSLOB: - p->name = "GAUSS_LOBATTO"; - return intxx_get_gausslob_info(p); - default: - return INTXX_INVALID_QUAD; - } -} - +/**********************************************************/ +/****** Runtime Generator for Radial Quadrature Info ******/ +/**********************************************************/ int intxx_get_radq_info(intxx_quad_info_type* p, int quad) { if(p == NULL) return INTXX_NULL_INFOPTR; @@ -299,67 +67,9 @@ int intxx_get_radq_info(intxx_quad_info_type* p, int quad) { return INTXX_SUCCESS; } -/******************************************************/ -/****** Info generation for quads w/o parameters ******/ -/******************************************************/ - -int intxx_noparam_info(intxx_quad_info_type* p, - int (*g)(intxx_quad_type*), - int (*d)(intxx_quad_type*)) { - if(p == NULL) return INTXX_NULL_INFOPTR; - - p->ext_params.n = 0; /// No External Parameters - p->generate = g; - p->destroy = d; - - return INTXX_SUCCESS; -} - -#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_NOPARAM_GET_INFO_IMPL(uniform); -INTXX_NOPARAM_GET_INFO_IMPL(gaussleg); -INTXX_NOPARAM_GET_INFO_IMPL(gausslob); -INTXX_NOPARAM_GET_INFO_IMPL(gausscheb1); -INTXX_NOPARAM_GET_INFO_IMPL(gausscheb2); -INTXX_NOPARAM_GET_INFO_IMPL(gausscheb2m); -INTXX_NOPARAM_GET_INFO_IMPL(gausscheb3); - -#undef INTXX_NOPARAM_GET_INFO_IMPL - -/********************************************************/ -/****** Info generation for quads w rad parameters ******/ -/********************************************************/ - - -int intxx_radscal_info(intxx_quad_info_type* p, - int (*sn)(intxx_quad_type*,const char*,double), - int (*gn)(intxx_quad_type*,const char*,double*), - int (*gparm)(intxx_quad_type*), - int (*dparm)(intxx_quad_type*), - int (*gquad)(intxx_quad_type*), - int (*dquad)(intxx_quad_type*)) { - if(p == NULL) return INTXX_NULL_INFOPTR; - - static const char* names[] = {"RAD_SCAL"}; - static const char* desc[] = {"Radial Scaling Factor"}; - p->ext_params = { - 1, names, desc, - sn, /* set_name */ - gn, /* get_name */ - gparm, /* generate */ - dparm /* destroy */ - }; - - p->generate = gquad; - p->destroy = dquad; - - return INTXX_SUCCESS; -} +/****************************************************/ +/****** Info generation for Radial Quadratures ******/ +/****************************************************/ #define INTXX_RADSCAL_GET_INFO_IMPL(cname) \ int intxx_get_##cname##_quad_info(intxx_quad_info_type* p) { \ @@ -380,7 +90,6 @@ INTXX_RADSCAL_GET_INFO_IMPL(mk) #undef INTXX_RADSCAL_GET_INFO_IMPL - /***********************************************/ /****** Allocate/Deallocate Radial Traits ******/ /***********************************************/ @@ -393,11 +102,6 @@ int intxx_destroy_##cname##_quad_params(intxx_quad_type* p) { \ return intxx_destroy_params_impl(p); \ } -using becke_traits_type = IntegratorXX::BeckeRadialTraits; -using mhl_traits_type = IntegratorXX::MurrayHandyLamingRadialTraits<2>; -using ta_traits_type = IntegratorXX::TreutlerAhlrichsRadialTraits; -using mk_traits_type = IntegratorXX::MuraKnowlesRadialTraits; - INTXX_RAD_TRAITS_IMPL(becke, becke_traits_type); INTXX_RAD_TRAITS_IMPL(mhl, mhl_traits_type ); INTXX_RAD_TRAITS_IMPL(ta, ta_traits_type ); @@ -408,6 +112,7 @@ INTXX_RAD_TRAITS_IMPL(mk, mk_traits_type ); /**********************************/ /****** Radial Quad GET_NAME ******/ /**********************************/ + #define INTXX_RAD_GETSET_IMPL(cname, traits_type) \ int intxx_get_name_##cname##_quad(intxx_quad_type* p, const char* name, double *v) { \ if(strcmp(name, "RAD_SCAL")){ return INTXX_INVALID_OPT; } \ @@ -423,31 +128,11 @@ INTXX_RAD_GETSET_IMPL(mhl, mhl_traits_type ); INTXX_RAD_GETSET_IMPL(ta, ta_traits_type ); INTXX_RAD_GETSET_IMPL(mk, mk_traits_type ); -#define INTXX_GD_BASIC_IMPL(cname, qname, ...) \ -int intxx_generate_##cname(intxx_quad_type* p) { \ - return intxx_generate_noparam_impl(p, ##__VA_ARGS__); \ -} \ -int intxx_destroy_##cname(intxx_quad_type* p) { \ - return intxx_destroy_impl(p); \ -} - -using uniform_quad_type = IntegratorXX::UniformTrapezoid; -using gaussleg_quad_type = IntegratorXX::GaussLegendre; -using gausslob_quad_type = IntegratorXX::GaussLobatto; -using gausscheb1_quad_type = IntegratorXX::GaussChebyshev1; -using gausscheb2_quad_type = IntegratorXX::GaussChebyshev2; -using gausscheb3_quad_type = IntegratorXX::GaussChebyshev3; -using gausscheb2m_quad_type = IntegratorXX::GaussChebyshev2Modified; +#undef INTXX_RAD_GETSET_IMPL -INTXX_GD_BASIC_IMPL(uniform, uniform_quad_type, 0.0, 1.0); -INTXX_GD_BASIC_IMPL(gaussleg, gaussleg_quad_type ); -INTXX_GD_BASIC_IMPL(gausslob, gausslob_quad_type ); -INTXX_GD_BASIC_IMPL(gausscheb1, gausscheb1_quad_type ); -INTXX_GD_BASIC_IMPL(gausscheb2, gausscheb2_quad_type ); -INTXX_GD_BASIC_IMPL(gausscheb3, gausscheb3_quad_type ); -INTXX_GD_BASIC_IMPL(gausscheb2m, gausscheb2m_quad_type); - -#undef INTXX_GD_BASIC_IMPL +/**********************************************************/ +/****** Allocate/Deallocate Radial Quadrature States ******/ +/**********************************************************/ #define INTXX_GD_RAD_IMPL(cname, qname, tname) \ int intxx_generate_##cname##_quad(intxx_quad_type* p) { \ @@ -457,14 +142,11 @@ int intxx_destroy_##cname##_quad(intxx_quad_type* p) { \ return intxx_destroy_impl(p); \ } -using becke_quad_type = IntegratorXX::Becke; -using mhl_quad_type = IntegratorXX::MurrayHandyLaming; -using ta_quad_type = IntegratorXX::TreutlerAhlrichs; -using mk_quad_type = IntegratorXX::MuraKnowles; - INTXX_GD_RAD_IMPL(becke, becke_quad_type, becke_traits_type); INTXX_GD_RAD_IMPL(mhl, mhl_quad_type, mhl_traits_type); INTXX_GD_RAD_IMPL(ta, ta_quad_type, ta_traits_type); INTXX_GD_RAD_IMPL(mk, mk_quad_type, mk_traits_type); +#undef INTXX_GD_RAD_IMPL + } From f0c21cbfd895cee4109ba931855d68d0145c3eb3 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 25 Sep 2023 10:18:07 -0700 Subject: [PATCH 18/22] Added dim info to C API, added additional meta data checks to C API UTs --- include/integratorxx/c_api.h | 1 + src/c_api/c_misc.cxx | 1 + src/c_api/c_primitive.cxx | 1 + src/c_api/c_radial.cxx | 1 + test/c_api.cxx | 38 +++++++++++++++++++++++++----------- 5 files changed, 31 insertions(+), 11 deletions(-) diff --git a/include/integratorxx/c_api.h b/include/integratorxx/c_api.h index df1cab9..d818b11 100644 --- a/include/integratorxx/c_api.h +++ b/include/integratorxx/c_api.h @@ -63,6 +63,7 @@ typedef struct { 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 diff --git a/src/c_api/c_misc.cxx b/src/c_api/c_misc.cxx index 639b797..de718de 100644 --- a/src/c_api/c_misc.cxx +++ b/src/c_api/c_misc.cxx @@ -6,6 +6,7 @@ extern "C" { void intxx_default_quad_info(intxx_quad_info_type* p) { p->number = 0; // Invalid p->kind = 0; // Invalid + p->dim = 0; // Invalid p->name = "Default"; // Default state p->generate = NULL; // Invalid p->destroy = NULL; // Invalid diff --git a/src/c_api/c_primitive.cxx b/src/c_api/c_primitive.cxx index 21fc568..6015348 100644 --- a/src/c_api/c_primitive.cxx +++ b/src/c_api/c_primitive.cxx @@ -43,6 +43,7 @@ int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { quad = quad & INTXX_PRMQ_MASK; p->number = quad; p->kind = INTXX_PRM_QUAD; + p->dim = 1; switch(quad) { case INTXX_PRMQ_UNIFORM: p->name = "UNIFORM"; diff --git a/src/c_api/c_radial.cxx b/src/c_api/c_radial.cxx index 076ff28..f736755 100644 --- a/src/c_api/c_radial.cxx +++ b/src/c_api/c_radial.cxx @@ -48,6 +48,7 @@ int intxx_get_radq_info(intxx_quad_info_type* p, int quad) { quad = quad & INTXX_RADQ_MASK; p->number = quad; p->kind = INTXX_RAD_QUAD; + p->dim = 1; switch(quad) { case INTXX_RADQ_BECKE: p->name = "BECKE"; diff --git a/test/c_api.cxx b/test/c_api.cxx index deea65d..7a5010f 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -26,61 +26,69 @@ TEST_CASE("C API") { using base_quad_type = QuadratureBase, std::vector>; std::unique_ptr base_quad = nullptr; + int quad_num; SECTION("Uniform") { using quad_type = UniformTrapezoid; - error = intxx_quad_init(&quad, INTXX_PRMQ_UNIFORM); + quad_num = INTXX_PRMQ_UNIFORM; name = "UNIFORM"; base_quad = std::make_unique(base_npts, 0.0, 1.0); } SECTION("Gauss-Legendre") { using quad_type = GaussLegendre; - error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSLEG); + quad_num = INTXX_PRMQ_GAUSSLEG; name = "GAUSS_LEGENDRE"; base_quad = std::make_unique(base_npts); } SECTION("Gauss-Lobatto") { using quad_type = GaussLobatto; - error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSLOB); + quad_num = INTXX_PRMQ_GAUSSLOB; name = "GAUSS_LOBATTO"; base_quad = std::make_unique(base_npts); } SECTION("Gauss-Chebyshev 1") { using quad_type = GaussChebyshev1; - error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_1); + quad_num = INTXX_PRMQ_GAUSSCHEB_1; name = "GAUSS_CHEBYSHEV_1"; base_quad = std::make_unique(base_npts); } SECTION("Gauss-Chebyshev 2") { using quad_type = GaussChebyshev2; - error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_2); + quad_num = INTXX_PRMQ_GAUSSCHEB_2; name = "GAUSS_CHEBYSHEV_2"; base_quad = std::make_unique(base_npts); } SECTION("Gauss-Chebyshev 2MOD") { using quad_type = GaussChebyshev2Modified; - error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_2MOD); + quad_num = INTXX_PRMQ_GAUSSCHEB_2MOD; name = "GAUSS_CHEBYSHEV_2MOD"; base_quad = std::make_unique(base_npts); } SECTION("Gauss-Chebyshev 3") { using quad_type = GaussChebyshev3; - error = intxx_quad_init(&quad, INTXX_PRMQ_GAUSSCHEB_3); + quad_num = INTXX_PRMQ_GAUSSCHEB_3; name = "GAUSS_CHEBYSHEV_3"; base_quad = std::make_unique(base_npts); } + // Initialize + error = intxx_quad_init(&quad, quad_num); REQUIRE(error == INTXX_SUCCESS); // Meta data REQUIRE(quad.info != NULL); REQUIRE(quad.npoints == -1); REQUIRE(quad._state_quad == NULL); + REQUIRE(quad.info->number == quad_num); + REQUIRE(quad.info->kind == INTXX_PRM_QUAD); + REQUIRE(quad.info->dim == 1); + REQUIRE(quad.info->generate != NULL); + REQUIRE(quad.info->destroy != NULL); REQUIRE(quad.info->ext_params.n == 0); REQUIRE(!strcmp(quad.info->name, name)); @@ -127,10 +135,11 @@ TEST_CASE("C API") { std::unique_ptr base_quad_default = nullptr; std::unique_ptr base_quad_scaled = nullptr; + int quad_num; SECTION("Becke") { using quad_type = Becke; using traits_type = typename quad_type::traits_type; - error = intxx_quad_init(&quad, INTXX_RADQ_BECKE); + quad_num = INTXX_RADQ_BECKE; name = "BECKE"; base_quad_default = std::make_unique(base_npts); base_quad_scaled = std::make_unique(base_npts,traits_type(RSCAL)); @@ -138,31 +147,38 @@ TEST_CASE("C API") { SECTION("MHL") { using quad_type = MurrayHandyLaming; - error = intxx_quad_init(&quad, INTXX_RADQ_MHL); + quad_num = INTXX_RADQ_MHL; name = "MURRAY_HANDY_LAMING"; base_quad_default = std::make_unique(base_npts); } SECTION("TA") { using quad_type = TreutlerAhlrichs; - error = intxx_quad_init(&quad, INTXX_RADQ_TA); + quad_num = INTXX_RADQ_TA; name = "TREUTLER_AHLRICHS"; base_quad_default = std::make_unique(base_npts); } SECTION("MK") { using quad_type = MuraKnowles; - error = intxx_quad_init(&quad, INTXX_RADQ_MK); + quad_num = INTXX_RADQ_MK; name = "MURA_KNOWLES"; base_quad_default = std::make_unique(base_npts); } + // Initialize + error = intxx_quad_init(&quad, quad_num); REQUIRE(error == INTXX_SUCCESS); // Meta data REQUIRE(quad.info != NULL); REQUIRE(quad.npoints == -1); REQUIRE(quad._state_quad == NULL); + REQUIRE(quad.info->number == quad_num); + REQUIRE(quad.info->kind == INTXX_RAD_QUAD); + REQUIRE(quad.info->dim == 1); + REQUIRE(quad.info->generate != NULL); + REQUIRE(quad.info->destroy != NULL); REQUIRE(!strcmp(quad.info->name, name)); From 8c02fc5ac9d6c5cb619ed71d016080dba4715c5b Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 25 Sep 2023 10:19:18 -0700 Subject: [PATCH 19/22] Added parameter UTs for remainer of radial quadratures --- test/c_api.cxx | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/c_api.cxx b/test/c_api.cxx index 7a5010f..5d3db00 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -147,23 +147,29 @@ TEST_CASE("C API") { SECTION("MHL") { using quad_type = MurrayHandyLaming; + using traits_type = typename quad_type::traits_type; quad_num = INTXX_RADQ_MHL; name = "MURRAY_HANDY_LAMING"; base_quad_default = std::make_unique(base_npts); + base_quad_scaled = std::make_unique(base_npts,traits_type(RSCAL)); } SECTION("TA") { using quad_type = TreutlerAhlrichs; + using traits_type = typename quad_type::traits_type; quad_num = INTXX_RADQ_TA; name = "TREUTLER_AHLRICHS"; base_quad_default = std::make_unique(base_npts); + base_quad_scaled = std::make_unique(base_npts,traits_type(RSCAL)); } SECTION("MK") { using quad_type = MuraKnowles; + using traits_type = typename quad_type::traits_type; quad_num = INTXX_RADQ_MK; name = "MURA_KNOWLES"; base_quad_default = std::make_unique(base_npts); + base_quad_scaled = std::make_unique(base_npts,traits_type(RSCAL)); } // Initialize From 775eaf8053ae3b5ea67290fd04041ba81687c91d Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 25 Sep 2023 11:12:25 -0700 Subject: [PATCH 20/22] Starting exposure of angular grids to C API --- src/c_api/CMakeLists.txt | 2 +- src/c_api/c_angular.cxx | 95 ++++++++++++++++++++++++++++++++++++++++ src/c_api/c_api.c | 2 +- src/c_api/c_misc.cxx | 1 + test/c_api.cxx | 87 ++++++++++++++++++++++++++++++++++++ 5 files changed, 185 insertions(+), 2 deletions(-) create mode 100644 src/c_api/c_angular.cxx diff --git a/src/c_api/CMakeLists.txt b/src/c_api/CMakeLists.txt index 16401c1..06c0cbb 100644 --- a/src/c_api/CMakeLists.txt +++ b/src/c_api/CMakeLists.txt @@ -1 +1 @@ -add_library( integratorxx c_api.c c_primitive.cxx c_radial.cxx c_misc.cxx ) +add_library( integratorxx c_api.c c_primitive.cxx c_radial.cxx c_misc.cxx c_angular.cxx ) diff --git a/src/c_api/c_angular.cxx b/src/c_api/c_angular.cxx new file mode 100644 index 0000000..0144272 --- /dev/null +++ b/src/c_api/c_angular.cxx @@ -0,0 +1,95 @@ +#include + +#include "c_internal.h" +#include "c_api_util.hpp" + + +extern "C" { + +using leb_quad_type = IntegratorXX::LebedevLaikov; +using delley_quad_type = IntegratorXX::Delley; +using wom_quad_type = IntegratorXX::Womersley; +using ab_quad_type = IntegratorXX::AhrensBeylkin; + +/*******************************************************/ +/****** 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; +} + +/*************************************************/ +/****** 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_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(p, ##__VA_ARGS__); \ +} \ +int intxx_destroy_##cname(intxx_quad_type* p) { \ + return intxx_destroy_impl(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 +} diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c index f6ea373..2032f4f 100644 --- a/src/c_api/c_api.c +++ b/src/c_api/c_api.c @@ -48,7 +48,7 @@ int intxx_quad_init(intxx_quad_type* p, int quad) { error = intxx_get_radq_info(finfo, quad); } else { // Angular by exclusion - // TODO: Get angular quadrature info + error = intxx_get_angq_info(finfo, quad); } if(error) return error; diff --git a/src/c_api/c_misc.cxx b/src/c_api/c_misc.cxx index de718de..f4d9e1e 100644 --- a/src/c_api/c_misc.cxx +++ b/src/c_api/c_misc.cxx @@ -8,6 +8,7 @@ void intxx_default_quad_info(intxx_quad_info_type* p) { p->kind = 0; // Invalid p->dim = 0; // Invalid p->name = "Default"; // Default state + p->ext_params.n = 0; // To initialize *everything* p->generate = NULL; // Invalid p->destroy = NULL; // Invalid } diff --git a/test/c_api.cxx b/test/c_api.cxx index 5d3db00..f7c51b5 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -262,6 +262,93 @@ TEST_CASE("C API") { } } + SECTION("Angular") { + const char* name; + using base_quad_type = QuadratureBase>, std::vector>; + std::unique_ptr base_quad_default = nullptr; + std::unique_ptr base_quad_other = nullptr; + + int default_npts, other_npts; + int quad_num; + SECTION("LebedevLaikov") { + using quad_type = LebedevLaikov; + quad_num = INTXX_ANGQ_LEB; + name = "LEBEDEV_LAIKOV"; + } + + SECTION("Delley") { + using quad_type = Delley; + quad_num = INTXX_ANGQ_DEL; + name = "DELLEY"; + } + + SECTION("AB") { + using quad_type = AhrensBeylkin; + quad_num = INTXX_ANGQ_AB; + name = "AHRENS_BEYLKIN"; + } + + SECTION("WOM") { + using quad_type = Womersley; + quad_num = INTXX_ANGQ_WOM; + name = "WOMERSLEY"; + } + + // Initialize + error = intxx_quad_init(&quad, quad_num); + REQUIRE(error == INTXX_SUCCESS); + + // Meta data + REQUIRE(quad.info != NULL); + REQUIRE(quad.npoints == -1); + REQUIRE(quad._state_quad == NULL); + REQUIRE(quad.info->number == quad_num); + REQUIRE(quad.info->kind == INTXX_ANG_QUAD); + REQUIRE(quad.info->dim == 3); + REQUIRE(quad.info->generate != NULL); + REQUIRE(quad.info->destroy != NULL); + REQUIRE(!strcmp(quad.info->name, name)); + + + // Check default parameters + REQUIRE(quad.info->ext_params.n == 0); + + // Get before set + int npts; + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_INVALID_OUT); + REQUIRE(npts == -1); + +#if 0 + // Set NPTS + error = intxx_quad_set_npts(&quad, base_npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(quad.npoints == base_npts); + + // Get NPTS + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_SUCCESS); + REQUIRE(npts == base_npts); + + // Check Quadrature Generation and Destruction + if(base_quad_default) { + intxx_generate_quad(&quad); + REQUIRE(quad._state_quad != NULL); + + // Check validity of the state + auto state_as_quad = reinterpret_cast(quad._state_quad); + REQUIRE(state_as_quad->npts() == npts); + for(auto i = 0; i < npts; ++ i) { + REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad_default->points()[i], 1e-15)); + REQUIRE_THAT(state_as_quad->weights()[i], Matchers::WithinAbs(name, base_quad_default->weights()[i], 1e-15)); + } + + intxx_destroy_quad(&quad); + REQUIRE(quad._state_quad == NULL); + } +#endif + } + intxx_quad_end(&quad); } From 16b382b102405cd220d0e3685177115b61b95c32 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 25 Sep 2023 11:34:36 -0700 Subject: [PATCH 21/22] Expose custom npts setters in C API for angular grids + UTs --- include/integratorxx/c_api.h | 1 + src/c_api/c_angular.cxx | 24 +++++++++++++++++++++++- src/c_api/c_api.c | 9 +++++++-- src/c_api/c_internal.h | 3 ++- src/c_api/c_misc.cxx | 5 ++++- src/c_api/c_primitive.cxx | 2 +- test/c_api.cxx | 26 +++++++++++++++++++++----- 7 files changed, 59 insertions(+), 11 deletions(-) diff --git a/include/integratorxx/c_api.h b/include/integratorxx/c_api.h index d818b11..a6e1609 100644 --- a/include/integratorxx/c_api.h +++ b/include/integratorxx/c_api.h @@ -72,6 +72,7 @@ typedef struct { 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 { diff --git a/src/c_api/c_angular.cxx b/src/c_api/c_angular.cxx index 0144272..76d572f 100644 --- a/src/c_api/c_angular.cxx +++ b/src/c_api/c_angular.cxx @@ -11,6 +11,7 @@ using delley_quad_type = IntegratorXX::Delley; using wom_quad_type = IntegratorXX::Womersley; using ab_quad_type = IntegratorXX::AhrensBeylkin; + /*******************************************************/ /****** Forward declare C-API for ANG Quadratures ******/ /*******************************************************/ @@ -57,6 +58,27 @@ int intxx_get_angq_info(intxx_quad_info_type* p, int 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 ******/ /*************************************************/ @@ -64,7 +86,7 @@ int intxx_get_angq_info(intxx_quad_info_type* p, int quad) { #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_destroy_##cname, intxx_##cname##_set_npts); \ } INTXX_NOPARAM_GET_INFO_IMPL(leb); diff --git a/src/c_api/c_api.c b/src/c_api/c_api.c index 2032f4f..abe3e34 100644 --- a/src/c_api/c_api.c +++ b/src/c_api/c_api.c @@ -87,8 +87,13 @@ int intxx_quad_set_npts(intxx_quad_type* p, int npts) { if(npts <= 0) return INTXX_INVALID_ARG; // TODO: Handle the case when NPTS is derived from params - p->npoints = npts; - return INTXX_SUCCESS; + if(p->info->set_npts) { + // Some quadratures have rules about how to set npts (e.g. angular grids) + return p->info->set_npts(p, npts); + } else { + p->npoints = npts; + return INTXX_SUCCESS; + } } int intxx_quad_get_npts(intxx_quad_type* p, int* npts) { diff --git a/src/c_api/c_internal.h b/src/c_api/c_internal.h index b7b62a8..748df83 100644 --- a/src/c_api/c_internal.h +++ b/src/c_api/c_internal.h @@ -13,7 +13,8 @@ int intxx_get_prmq_info(intxx_quad_info_type* p, int quad); int intxx_noparam_info(intxx_quad_info_type* p, int (*g)(intxx_quad_type*), - int (*d)(intxx_quad_type*)); + int (*d)(intxx_quad_type*), + int (*s)(intxx_quad_type*, int)); int intxx_radscal_info(intxx_quad_info_type* p, int (*sn)(intxx_quad_type*,const char*,double), diff --git a/src/c_api/c_misc.cxx b/src/c_api/c_misc.cxx index f4d9e1e..3279028 100644 --- a/src/c_api/c_misc.cxx +++ b/src/c_api/c_misc.cxx @@ -19,12 +19,14 @@ void intxx_default_quad_info(intxx_quad_info_type* p) { int intxx_noparam_info(intxx_quad_info_type* p, int (*g)(intxx_quad_type*), - int (*d)(intxx_quad_type*)) { + int (*d)(intxx_quad_type*), + int (*s)(intxx_quad_type*, int)) { if(p == NULL) return INTXX_NULL_INFOPTR; p->ext_params.n = 0; /// No External Parameters p->generate = g; p->destroy = d; + p->set_npts = s; return INTXX_SUCCESS; } @@ -54,6 +56,7 @@ int intxx_radscal_info(intxx_quad_info_type* p, p->generate = gquad; p->destroy = dquad; + p->set_npts = NULL; return INTXX_SUCCESS; } diff --git a/src/c_api/c_primitive.cxx b/src/c_api/c_primitive.cxx index 6015348..9cd866e 100644 --- a/src/c_api/c_primitive.cxx +++ b/src/c_api/c_primitive.cxx @@ -78,7 +78,7 @@ int intxx_get_prmq_info(intxx_quad_info_type* p, int quad) { #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_destroy_##cname, NULL); \ } INTXX_NOPARAM_GET_INFO_IMPL(uniform); diff --git a/test/c_api.cxx b/test/c_api.cxx index f7c51b5..cfc5b67 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -274,24 +274,32 @@ TEST_CASE("C API") { using quad_type = LebedevLaikov; quad_num = INTXX_ANGQ_LEB; name = "LEBEDEV_LAIKOV"; + default_npts = 302; + other_npts = 974; } SECTION("Delley") { using quad_type = Delley; quad_num = INTXX_ANGQ_DEL; name = "DELLEY"; + default_npts = 302; + other_npts = 974; } SECTION("AB") { using quad_type = AhrensBeylkin; quad_num = INTXX_ANGQ_AB; name = "AHRENS_BEYLKIN"; + default_npts = 372; + other_npts = 972; } SECTION("WOM") { using quad_type = Womersley; quad_num = INTXX_ANGQ_WOM; name = "WOMERSLEY"; + default_npts = 393; + other_npts = 969; } // Initialize @@ -319,17 +327,25 @@ TEST_CASE("C API") { REQUIRE(error == INTXX_INVALID_OUT); REQUIRE(npts == -1); -#if 0 - // Set NPTS - error = intxx_quad_set_npts(&quad, base_npts); + // Set NPTS (incorrect) + error = intxx_quad_set_npts(&quad, default_npts+1); + REQUIRE(error == INTXX_INVALID_ARG); + error = intxx_quad_get_npts(&quad, &npts); + REQUIRE(error == INTXX_INVALID_OUT); + REQUIRE(npts == -1); + + // Set NPTS (correct) + error = intxx_quad_set_npts(&quad, default_npts); REQUIRE(error == INTXX_SUCCESS); - REQUIRE(quad.npoints == base_npts); + REQUIRE(quad.npoints == default_npts); + // Get NPTS error = intxx_quad_get_npts(&quad, &npts); REQUIRE(error == INTXX_SUCCESS); - REQUIRE(npts == base_npts); + REQUIRE(npts == default_npts); +#if 0 // Check Quadrature Generation and Destruction if(base_quad_default) { intxx_generate_quad(&quad); From 7b7c3c633b8a679c0e06141f6b6ea372423f9444 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Mon, 25 Sep 2023 11:41:10 -0700 Subject: [PATCH 22/22] Check angular quadrature generation in C API --- test/c_api.cxx | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/test/c_api.cxx b/test/c_api.cxx index cfc5b67..2afe0a0 100644 --- a/test/c_api.cxx +++ b/test/c_api.cxx @@ -266,7 +266,7 @@ TEST_CASE("C API") { const char* name; using base_quad_type = QuadratureBase>, std::vector>; std::unique_ptr base_quad_default = nullptr; - std::unique_ptr base_quad_other = nullptr; + //std::unique_ptr base_quad_other = nullptr; int default_npts, other_npts; int quad_num; @@ -284,6 +284,8 @@ TEST_CASE("C API") { name = "DELLEY"; default_npts = 302; other_npts = 974; + base_quad_default = std::make_unique(default_npts); + //base_quad_other = std::make_unique(other_npts); } SECTION("AB") { @@ -292,6 +294,8 @@ TEST_CASE("C API") { name = "AHRENS_BEYLKIN"; default_npts = 372; other_npts = 972; + base_quad_default = std::make_unique(default_npts); + //base_quad_other = std::make_unique(other_npts); } SECTION("WOM") { @@ -300,6 +304,8 @@ TEST_CASE("C API") { name = "WOMERSLEY"; default_npts = 393; other_npts = 969; + base_quad_default = std::make_unique(default_npts); + //base_quad_other = std::make_unique(other_npts); } // Initialize @@ -339,13 +345,11 @@ TEST_CASE("C API") { REQUIRE(error == INTXX_SUCCESS); REQUIRE(quad.npoints == default_npts); - // Get NPTS error = intxx_quad_get_npts(&quad, &npts); REQUIRE(error == INTXX_SUCCESS); REQUIRE(npts == default_npts); -#if 0 // Check Quadrature Generation and Destruction if(base_quad_default) { intxx_generate_quad(&quad); @@ -355,14 +359,14 @@ TEST_CASE("C API") { auto state_as_quad = reinterpret_cast(quad._state_quad); REQUIRE(state_as_quad->npts() == npts); for(auto i = 0; i < npts; ++ i) { - REQUIRE_THAT(state_as_quad->points()[i], Matchers::WithinAbs(name, base_quad_default->points()[i], 1e-15)); - REQUIRE_THAT(state_as_quad->weights()[i], Matchers::WithinAbs(name, base_quad_default->weights()[i], 1e-15)); + REQUIRE(state_as_quad->points()[i] == base_quad_default->points()[i]); + REQUIRE(state_as_quad->weights()[i] == base_quad_default->weights()[i]); } intxx_destroy_quad(&quad); REQUIRE(quad._state_quad == NULL); } -#endif + } intxx_quad_end(&quad);