From 26d01a06882274fad6584c827caa160fc9620526 Mon Sep 17 00:00:00 2001 From: hmb Date: Sat, 28 Sep 2019 20:36:41 -0700 Subject: [PATCH 1/2] Implemented qawc wrapper --- src/gsl/mod.rs | 3 ++ src/gsl/qawc.rs | 84 +++++++++++++++++++++++++++++++++++++++++++++++++ src/gsl/test.rs | 21 ++++++++++++- 3 files changed, 107 insertions(+), 1 deletion(-) create mode 100644 src/gsl/qawc.rs diff --git a/src/gsl/mod.rs b/src/gsl/mod.rs index 4611bd4..58ba07a 100644 --- a/src/gsl/mod.rs +++ b/src/gsl/mod.rs @@ -52,6 +52,9 @@ pub use self::qng::QNG; mod qag; pub use self::qag::{QAG, QAGRule}; +mod qawc; +pub use self::qawc::QAWC; + mod qags; pub use self::qags::QAGS; diff --git a/src/gsl/qawc.rs b/src/gsl/qawc.rs new file mode 100644 index 0000000..b3c4823 --- /dev/null +++ b/src/gsl/qawc.rs @@ -0,0 +1,84 @@ +use ::bindings; +use ::{IntegrationResult, Integrator, Real}; +use ::ffi::LandingPad; +use ::traits::{IntegrandInput, IntegrandOutput}; + +use super::{make_gsl_function, GSLIntegrationError, GSLIntegrationWorkspace}; + +/// Quadrature Adaptive Weighed integration for Cauchy principal values +/// Calculates the principal value of a function f(x) with a singular weight +/// 1 / (x - c) centered at a point c +#[derive(Debug, Clone)] +pub struct QAWC { + range_low: Real, + range_high: Real, + c: Real, + wkspc: GSLIntegrationWorkspace +} + +impl QAWC { + /// Creates a new QAGS with enough memory for `nintervals` subintervals. + /// This will create a QAG to integrate the range [0, 1]. + /// To change the integration bounds, see `with_range`. + /// The singularity is set to 0.5 by default and can be changed with `with_singularity` + pub fn new(nintervals: usize) -> Self{ + QAWC { + range_low: 0.0, + range_high: 1.0, + c: 0.5, + wkspc: GSLIntegrationWorkspace::new(nintervals) + } + } + + /// Discards the old workspace and allocates a new one with enough memory + /// for `nintervals` subintervals. + pub fn with_nintervals(self, nintervals: usize) -> Self{ + QAWC { + wkspc: GSLIntegrationWorkspace::new(nintervals), + ..self + } + } + + pub fn with_range(self, range_low: Real, range_high: Real) -> Self { + QAWC { range_low, range_high, ..self } + } + + pub fn with_singularity(self, c: Real) -> Self{ + QAWC {c, ..self} + } +} + +impl Integrator for QAWC { + type Success = IntegrationResult; + type Failure = GSLIntegrationError; + + fn integrate B>(&mut self, fun: F, epsrel: Real, epsabs: Real) -> Result + where A: IntegrandInput, + B: IntegrandOutput + { + let mut value: Real = 0.0; + let mut error: Real = 0.0; + + let mut lp = LandingPad::new(fun); + let retcode = unsafe { + let mut gslfn = make_gsl_function(&mut lp, self.range_low, self.range_high)?; + bindings::gsl_integration_qawc(&mut gslfn.function, + self.range_low, self.range_high, + self.c, + epsabs, epsrel, + self.wkspc.nintervals, + self.wkspc.wkspc, + &mut value, + &mut error) + }; + lp.maybe_resume_unwind(); + + if retcode != bindings::GSL_SUCCESS { + Err(GSLIntegrationError::GSLError(retcode.into())) + } else { + Ok(IntegrationResult { + value, error + }) + } + } +} \ No newline at end of file diff --git a/src/gsl/test.rs b/src/gsl/test.rs index 950f8ae..ee0703e 100644 --- a/src/gsl/test.rs +++ b/src/gsl/test.rs @@ -1,7 +1,7 @@ //use std::intrinsics::unchecked_div; use ::Real; use ::Integrator; -use super::{GSLIntegrationError, QNG, QAG, QAGS, QAGP}; +use super::{GSLIntegrationError, QNG, QAG, QAGS, QAGP, QAWC}; fn nan(_: Real) -> Real { ::std::f64::NAN @@ -124,3 +124,22 @@ fn test_qagp_singularity() { .expect("integration should succeed"); assert!((res2.value - 2f64).abs() <= res2.error); } + +#[test] +fn test_qawc(){ + let mut qawc = QAWC::new(1000) + .with_range(-1., 1.) + .with_singularity(0.0); + // Trivial example + let res1 = qawc.integrate(|x: Real| x, 1e-6, 1e-10) + .expect("integration should succeed"); + assert!((res1.value - 2.0f64).abs() <= res1.error); + + let mut qawc = qawc.with_range(0.0, 2.0) + .with_singularity(1.0); + // Analytical value: 2 E SinhIntegral[1] ~ 5.74781168 + let res2 = qawc.integrate(|x: Real| x.exp(), 1e-6, 1e-10) + .expect("integration should succeed"); + println!("Result: {}, err:{}", res2.value, res2.error); + assert!((res2.value - 5.747811685312522940687587).abs() <= res2.error); +} From 987f81a21011fd029d826b4a8974e2f50e04b268 Mon Sep 17 00:00:00 2001 From: hmb Date: Sat, 23 Nov 2019 01:28:38 -0800 Subject: [PATCH 2/2] Added ability to disable linking to gslcblas --- Cargo.toml | 1 + build.rs | 4 +++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index 9c83ec7..e6127fc 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -19,6 +19,7 @@ travis-ci = { repository = "jhod0/integrators", branch = "master" } default = ["cuba", "gsl"] cuba = [] gsl = [] +no_gslcblas=[] [build-dependencies] bindgen = "0.43.*" diff --git a/build.rs b/build.rs index d71d08f..dcae299 100644 --- a/build.rs +++ b/build.rs @@ -15,7 +15,9 @@ fn main() { if env::var("CARGO_FEATURE_GSL").is_ok() { println!("cargo:rustc-link-lib=gsl"); - println!("cargo:rustc-link-lib=gslcblas"); + if env::var("CARGO_FEATURE_NO_GSLCBLAS").is_err(){ + println!("cargo:rustc-link-lib=gslcblas"); + } bindings = bindings.header("gsl_wrapper.h"); any_bindings = true; }