From a15653799fd99eff2aa35737b2b183b9f64bd827 Mon Sep 17 00:00:00 2001 From: 9il Date: Fri, 7 Feb 2020 10:27:30 +0700 Subject: [PATCH] add mir.quadrature --- .gitignore | 6 + .travis.yml | 16 ++ README.md | 12 +- dub.sdl | 13 ++ dub.selections.json | 11 ++ meson.build | 60 ++++++ meson_options.txt | 1 + source/mir/quadrature.d | 339 +++++++++++++++++++++++++++++++++ subprojects/cblas.wrap | 4 + subprojects/lapack.wrap | 4 + subprojects/mir-algorithm.wrap | 4 + subprojects/mir-blas.wrap | 4 + subprojects/mir-core.wrap | 4 + subprojects/mir-lapack.wrap | 4 + 14 files changed, 480 insertions(+), 2 deletions(-) create mode 100644 .travis.yml create mode 100644 dub.sdl create mode 100644 dub.selections.json create mode 100644 meson.build create mode 100644 meson_options.txt create mode 100644 source/mir/quadrature.d create mode 100644 subprojects/cblas.wrap create mode 100644 subprojects/lapack.wrap create mode 100644 subprojects/mir-algorithm.wrap create mode 100644 subprojects/mir-blas.wrap create mode 100644 subprojects/mir-core.wrap create mode 100644 subprojects/mir-lapack.wrap diff --git a/.gitignore b/.gitignore index 74b926f..9789c6c 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,9 @@ docs/ # Code coverage *.lst +mir-integral-test-library +build +doc +docs +subprojects/*/ +mir-integral-test-test_travis diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..8c89789 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,16 @@ +sudo: false + +language: d + +d: + - dmd + - ldc + +addons: + apt: + packages: + - liblapack-dev + - libblas-dev + +script: + - dub test -c test_travis diff --git a/README.md b/README.md index cd1023d..9c282f6 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,10 @@ -# mir-integral -Integration routines +# Numerical Integration Routines + +## Contents + +### `mir.quadrature` + + - Gauss-Hermite Quadrature + - Gauss-Jacobi Quadrature + - Gauss-Laguerre Quadrature + - Gauss-Legendre Quadrature diff --git a/dub.sdl b/dub.sdl new file mode 100644 index 0000000..00e53b9 --- /dev/null +++ b/dub.sdl @@ -0,0 +1,13 @@ +name "mir-integral" +description "Integration routines" +authors "Ilya Yaroshenko" +copyright "Copyright © 2020 Symmetry Investments and Kaleidic Associates" +license "BSL-1.0" +dependency "mir-lapack" version=">=1.2.2" + +configuration "library" { +} + +configuration "test_travis" { + subConfiguration "mir-lapack" "blas" +} diff --git a/dub.selections.json b/dub.selections.json new file mode 100644 index 0000000..c79d754 --- /dev/null +++ b/dub.selections.json @@ -0,0 +1,11 @@ +{ + "fileVersion": 1, + "versions": { + "cblas": "2.0.4", + "lapack": "0.3.0", + "mir-algorithm": "3.7.17", + "mir-blas": "1.1.9", + "mir-core": "1.0.2", + "mir-lapack": "1.2.2" + } +} diff --git a/meson.build b/meson.build new file mode 100644 index 0000000..f2f30da --- /dev/null +++ b/meson.build @@ -0,0 +1,60 @@ +project('mir-integral', 'd', 'c', version : '1.0.0', license: 'BSL-1.0') + +cblas_d_dep = dependency('cblas-d', fallback : ['cblas', 'cblas_dep']) +lapack_d_dep = dependency('lapack-d', fallback : ['lapack', 'lapack_dep']) +mir_lapack_dep = dependency('mir-lapack', fallback : ['mir-lapack', 'mir_lapack_dep']) +mir_algorithm_dep = dependency('mir-algorithm', fallback : ['mir-algorithm', 'mir_algorithm_dep']) +mir_blas_dep = dependency('mir-blas', fallback : ['mir-blas', 'mir_blas_dep']) + +required_deps = [ + cblas_d_dep, + lapack_d_dep, + mir_algorithm_dep, + mir_blas_dep, + mir_lapack_dep, +] + +mir_integral_dir = include_directories('source/') + +mir_integral_src = [ + 'source/mir/quadrature.d', +] + +mir_integral_lib = library(meson.project_name(), + mir_integral_src, + include_directories: mir_integral_dir, + install: true, + version: meson.project_version(), + dependencies: required_deps, +) + +mir_integral_dep = declare_dependency( + link_with: [mir_integral_lib], + include_directories: mir_integral_dir, + dependencies: required_deps, +) + +install_subdir('source/', + strip_directory : true, + install_dir: 'include/d/' + meson.project_name(), +) + +import('pkgconfig').generate(mir_integral_lib, + description: 'Integration Routines', + subdirs: 'd/' + meson.project_name(), +) + +if get_option('with_test') + + mir_lapack_test_exe = executable(meson.project_name() + '-test', + mir_integral_src, + include_directories: mir_integral_dir, + d_unittest: true, + d_module_versions: ['mir_test'], + link_args: '-main', + dependencies: required_deps, + ) + + test(meson.project_name() + '-test', mir_lapack_test_exe) + +endif diff --git a/meson_options.txt b/meson_options.txt new file mode 100644 index 0000000..27602cf --- /dev/null +++ b/meson_options.txt @@ -0,0 +1 @@ +option('with_test', type : 'boolean', value : false) diff --git a/source/mir/quadrature.d b/source/mir/quadrature.d new file mode 100644 index 0000000..f1e5d77 --- /dev/null +++ b/source/mir/quadrature.d @@ -0,0 +1,339 @@ +/++ +This module contains betterC compatible quadrature computation routines. ++/ +module mir.quadrature; + +import mir.math.common: sqrt, exp; +import mir.math.constant: PI, LN2; + +@safe pure nothrow: + +/++ +Gauss-Hermite Quadrature + +Params: + x = (out) user-allocated quadrature nodes in ascending order length of `N` + w = (out) user-allocated corresponding quadrature weights length of `N` + work = (temp) user-allocated workspace length of greate or equal to `(N + 1) ^^ 2` + +Returns: 0 on success, `xSTEQR` LAPACK code on numerical error. + +See_also: $(LREF gaussQuadratureWork); ++/ +size_t gaussHermiteQuadrature(T)( + scope T[] x, + scope T[] w, + scope T[] work) @nogc +in { + assert(x.length == w.length); + if (x.length) + assert(work.length >= x.length ^^ 2); +} +do { + enum T mu0 = sqrt(PI); + foreach (i; 0 .. x.length) + { + x[i] = 0; + w[i] = T(0.5) * i; + } + return gaussQuadratureImpl!T(x, w, work, mu0, true); +} + +/// +unittest +{ + import mir.math.common; + import mir.ndslice.allocation; + + auto n = 5; + auto x = new double[n]; + auto w = new double[n]; + auto work = new double[(n + 1) ^^ 2]; + + gaussHermiteQuadrature(x, w, work); + + static immutable xc = + [-2.02018287, + -0.95857246, + 0. , + 0.95857246, + 2.02018287]; + + static immutable wc = + [0.01995324, + 0.39361932, + 0.94530872, + 0.39361932, + 0.01995324]; + + foreach (i; 0 .. n) + { + assert(x[i].approxEqual(xc[i])); + assert(w[i].approxEqual(wc[i])); + } +} + +/++ +Gauss-Jacobi Quadrature + +Params: + x = (out) user-allocated quadrature nodes in ascending order length of `N` + w = (out) user-allocated corresponding quadrature weights length of `N` + work = (temp) user-allocated workspace length of greate or equal to `(N + 1) ^^ 2` + alpha = parameter '> -1' + beta = parameter '> -1' + +Returns: 0 on success, `xSTEQR` LAPACK code on numerical error. + +See_also: $(LREF gaussQuadratureWork); ++/ +size_t gaussJacobiQuadrature(T)( + scope T[] x, + scope T[] w, + scope T[] work, + T alpha, + T beta) @nogc +in { + assert(T.infinity > alpha && alpha > -1); + assert(T.infinity > beta && beta > -1); + assert(x.length == w.length); + if (x.length) + assert(work.length >= x.length ^^ 2); +} +do { + if (x.length == 0) + return 0; + auto s = alpha + beta; + auto d = beta - alpha; + import core.stdc.tgmath: lgamma; + auto mu0 = exp(T(LN2) * (s + 1) + (lgamma(alpha + 1) + lgamma(beta + 1) - lgamma(s + 2))); + x[0] = d / (s + 2); + const sd = s * d; + foreach (i; 1 .. x.length) + { + const m_i = T(1) / i; + const q = (2 + s * m_i); + x[i] = sd * m_i ^^ 2 / (q * (2 + (s + 2) * m_i)); + w[i] = 4 * (1 + alpha * m_i) * (1 + beta * m_i) * (1 + s * m_i) + / ((2 + (s + 1) * m_i) * (2 + (s - 1) * m_i) * q ^^ 2); + } + return gaussQuadratureImpl!T(x, w, work, mu0); +} + +/// +unittest +{ + import mir.math.common; + import mir.ndslice.allocation; + + auto n = 5; + auto x = new double[n]; + auto w = new double[n]; + auto work = new double[(n + 1) ^^ 2]; + + gaussJacobiQuadrature(x, w, work, 2.3, 3.6); + + static immutable xc = + [-0.6553677 , + -0.29480426, + 0.09956621, + 0.47584565, + 0.78356514]; + + static immutable wc = + [0.02262392, + 0.19871672, + 0.43585107, + 0.32146619, + 0.0615342 ]; + + foreach (i; 0 .. n) + { + assert(x[i].approxEqual(xc[i])); + assert(w[i].approxEqual(wc[i])); + } +} + +/++ +Gauss-Laguerre Quadrature + +Params: + x = (out) user-allocated quadrature nodes in ascending order length of `N` + w = (out) user-allocated corresponding quadrature weights length of `N` + work = (temp) user-allocated workspace length of greate or equal to `(N + 1) ^^ 2` + alpha = (optional) parameter '> -1' + +Returns: 0 on success, `xSTEQR` LAPACK code on numerical error. + +See_also: $(LREF gaussQuadratureWork); ++/ +size_t gaussLaguerreQuadrature(T)( + scope T[] x, + scope T[] w, + scope T[] work, + T alpha = 0) @nogc +in { + assert(T.infinity > alpha && alpha > -1); + assert(x.length == w.length); + if (x.length) + assert(work.length >= x.length ^^ 2); +} +do { + import core.stdc.tgmath: tgamma; + auto mu0 = tgamma(alpha + 1); + foreach (i; 0 .. x.length) + { + x[i] = 2 * i + (1 + alpha); + w[i] = i * (i + alpha); + } + return gaussQuadratureImpl!T(x, w, work, mu0); +} + +/// +unittest +{ + import mir.math.common; + import mir.ndslice.allocation; + + auto n = 5; + auto x = new double[n]; + auto w = new double[n]; + auto work = new double[(n + 1) ^^ 2]; + + gaussLaguerreQuadrature(x, w, work); + + static immutable xc = + [ 0.26356032, + 1.41340306, + 3.59642577, + 7.08581001, + 12.64080084]; + + static immutable wc = + [5.21755611e-01, + 3.98666811e-01, + 7.59424497e-02, + 3.61175868e-03, + 2.33699724e-05]; + + foreach (i; 0 .. n) + { + assert(x[i].approxEqual(xc[i])); + assert(w[i].approxEqual(wc[i])); + } +} + +/++ +Gauss-Legendre Quadrature + +Params: + x = (out) user-allocated quadrature nodes in ascending order length of `N` + w = (out) user-allocated corresponding quadrature weights length of `N` + work = (temp) user-allocated workspace length of greate or equal to `(N + 1) ^^ 2` + +Returns: 0 on success, `xSTEQR` LAPACK code on numerical error. + +See_also: $(LREF gaussQuadratureWork); ++/ +size_t gaussLegendreQuadrature(T)( + scope T[] x, + scope T[] w, + scope T[] work) @nogc +in { + assert(x.length == w.length); + if (x.length) + assert(work.length >= x.length ^^ 2); +} +do { + if (x.length == 0) + return 0; + enum mu0 = 2; + x[0] = 0; + foreach (i; 1 .. x.length) + { + const m_i = T(1) / i; + x[i] = 0; + w[i] = 1 / (4 - m_i ^^ 2); + } + return gaussQuadratureImpl!T(x, w, work, mu0, true); +} + +/// +unittest +{ + import mir.math.common; + import mir.ndslice.allocation; + + auto n = 5; + auto x = new double[n]; + auto w = new double[n]; + auto work = new double[(n + 1) ^^ 2]; + + gaussLegendreQuadrature(x, w, work); + + static immutable xc = + [-0.90617985, + -0.53846931, + 0. , + 0.53846931, + 0.90617985]; + + static immutable wc = + [0.23692689, + 0.47862867, + 0.56888889, + 0.47862867, + 0.23692689]; + + foreach (i; 0 .. n) + { + assert(x[i].approxEqual(xc[i])); + assert(w[i].approxEqual(wc[i])); + } +} + +private size_t gaussQuadratureImpl(T)( + scope T[] alpha_x, + scope T[] beta_w, + scope T[] work, + double mu0, + bool symmetrize = false) @nogc +in { + assert(alpha_x.length == beta_w.length); + if (alpha_x.length) + assert(work.length >= (alpha_x.length + 1) ^^ 2); + foreach (ref b; beta_w[1 .. $]) + assert (T.infinity > b && b > 0); +} +do { + pragma(inline, false); + auto n = alpha_x.length; + if (n == 0) + return n; + foreach (ref b; beta_w[1 .. n]) + b = b.sqrt; + auto nq = n ^^ 2; + import mir.ndslice.slice: sliced; + import mir.ndslice.topology: canonical; + import mir.lapack: steqr; + auto z = work[0 .. nq].sliced(n, n); + auto info = steqr('I', alpha_x.sliced, beta_w[1 .. $].sliced, z.canonical, work[nq .. $].sliced); + foreach (i; 0 .. n) + beta_w[i] = z[i, 0] ^^ 2 * mu0; + if (symmetrize) + { + auto h = n / 2; + alias x = alpha_x; + alias w = beta_w; + foreach (i; 0 .. h) + { + x[i] = -(x[n - (i + 1)] = T(0.5) * (x[n - (i + 1)] - x[i])); + w[i] = +(w[n - (i + 1)] = T(0.5) * (w[n - (i + 1)] + w[i])); + } + if (n % 2) + { + x[h] = 0; + } + } + return info; +} diff --git a/subprojects/cblas.wrap b/subprojects/cblas.wrap new file mode 100644 index 0000000..509178c --- /dev/null +++ b/subprojects/cblas.wrap @@ -0,0 +1,4 @@ +[wrap-git] +directory=cblas +url=https://github.com/DlangScience/cblas.git +revision=head \ No newline at end of file diff --git a/subprojects/lapack.wrap b/subprojects/lapack.wrap new file mode 100644 index 0000000..b73e2be --- /dev/null +++ b/subprojects/lapack.wrap @@ -0,0 +1,4 @@ +[wrap-git] +directory=lapack +url=https://github.com/libmir/lapack.git +revision=head \ No newline at end of file diff --git a/subprojects/mir-algorithm.wrap b/subprojects/mir-algorithm.wrap new file mode 100644 index 0000000..94e57aa --- /dev/null +++ b/subprojects/mir-algorithm.wrap @@ -0,0 +1,4 @@ +[wrap-git] +directory=mir-algorithm +url=https://github.com/libmir/mir-algorithm.git +revision=head \ No newline at end of file diff --git a/subprojects/mir-blas.wrap b/subprojects/mir-blas.wrap new file mode 100644 index 0000000..9b7501b --- /dev/null +++ b/subprojects/mir-blas.wrap @@ -0,0 +1,4 @@ +[wrap-git] +directory=mir-blas +url=https://github.com/libmir/mir-blas.git +revision=head \ No newline at end of file diff --git a/subprojects/mir-core.wrap b/subprojects/mir-core.wrap new file mode 100644 index 0000000..8f2070b --- /dev/null +++ b/subprojects/mir-core.wrap @@ -0,0 +1,4 @@ +[wrap-git] +directory=mir-core +url=https://github.com/libmir/mir-core.git +revision=head \ No newline at end of file diff --git a/subprojects/mir-lapack.wrap b/subprojects/mir-lapack.wrap new file mode 100644 index 0000000..b919c04 --- /dev/null +++ b/subprojects/mir-lapack.wrap @@ -0,0 +1,4 @@ +[wrap-git] +directory=mir-lapack +url=https://github.com/libmir/mir-lapack.git +revision=head \ No newline at end of file