gaussjacobiquad

A permissively licensed set of Fortran routines for the Gauss Jacobi Quadrature

https://github.com/haozeke/gaussjacobiquad

Keywords

Repository

A permissively licensed set of Fortran routines for the Gauss Jacobi Quadrature

readme.org

* GaussJacobiQuad [[https://zenodo.org/badge/latestdoi/667604312][file:https://zenodo.org/badge/667604312.svg]] [[https://github.com/HaoZeke/GaussJacobiQuad/actions/workflows/build_test.yml][file:https://github.com/HaoZeke/GaussJacobiQuad/actions/workflows/build_test.yml/badge.svg]] [[https://github.com/HaoZeke/GaussJacobiQuad/actions/workflows/build_docs.yml][file:https://github.com/HaoZeke/GaussJacobiQuad/actions/workflows/build_docs.yml/badge.svg]]


** About
A permissively licensed modern implementation of Gauss-Jacobi quadrature which returns the weights and nodes over the standard interval ~[-1, 1]~.
** Usage
The most automated approach is to use the ~conda~ environment and ~fpm~ build:
#+begin_src bash
micromamba create -f environment.yml
micromamba activate gaussjacquad
fpm build
#+end_src

An analytic result can be obtained from the ~scripts~ folder.
#+begin_src bash
python scripts/sympy_gauss_jac.py --npts  --alpha  --beta  --n_dig 
fpm test
#+end_src

Running the implemented recursion based Gauss-Jacobi can be done via:
#+begin_src bash
fpm run gjp_quad_rec --   
fpm run gjp_quad --    
#+end_src

Currently the only supported methods are "rec" and "gw" (Golub Welsch).
*** Meson Support
A ~meson~ build backend is also present, which makes it easier to incorporate as subprojects of projects other than those supported by ~fpm~.

#+begin_src bash
meson setup bbdir
meson compile -C bbdir
./bbdir/gjp_quad    
#+end_src

*** External Fortran libaries
We support and encourage users to generate single file versions of the algorithms here to include in their code bases. This can be done with the ~scripts/add_headers.py~ script:
#+begin_src bash
# Strips comments by default
python scripts/export_single.py --modulename "gjp_algo665"
python scripts/export_single.py --modulename "gjp_gw" --keep-comments
#+end_src
These can be dropped into any code base or compiled as is into a shared library.
#+begin_src bash
gfortran dist/gjp_algo665_single.f90 --shared -fPIC -o libgjp_algo665.so
#+end_src
*** Interfaces
*** C/C++ header interface
We provide a header only interface, which bypasses passing strings between
Fortran and C. In order to do this, there is some duplication logic in
~GaussJacobiQuad~ and ~GaussJacobiQuadCInterp~.

There is also a CLI interface to the ~C~ bound interface, ~c_cli_gjpq~. However,
this will not be compiled by ~fpm~.

#+begin_src bash
meson setup bbdir
meson compile -C bbdir
./bbdir/c_cli_gjpq    
#+end_src

The ~C~ CLI might be more pleasant in that decimals do not need to be provided
explicitly for ~alpha~ and ~beta~.
*** f2py generated interface
The ~ISO_C_BINDING~ variant of the code is also used for a ~python~
interface generated with ~f2py~. It is easiest to use with the new ~meson~ back-end:
#+begin_src bash
cd interfaces/PyInterface
f2py -c --backend meson gjquadpy.pyf \
    --dep lapack \
    ../../src/GaussJacobiQuadCCompat.f90 \
    ../../src/GaussJacobiQuad.f90 \
    ../../src/gjp_constants.f90 \
    ../../src/gjp_types.f90 \
    ../../src/gjp_rec.f90 \
    ../../src/gjp_common.f90 \
    ../../src/gjp_lapack.f90 \
    ../../src/gjp_gw.f90 \
    ../../src/gjp_algo665.f90
#+end_src
Once compiled, the ~gjpquad_cli.py~ script can be used to run the code:
#+begin_src bash
python gjquad_cli.py --npts 5 --alpha 2 --beta 3
Root: -6.90457750126761027E-01 Weight: 2.74101780663370022E-02
Root: -3.26519931349000647E-01 Weight: 2.12917860603648035E-01
Root:  8.23378495520349085E-02 Weight: 4.39084379443950568E-01
Root:  4.75178870612831761E-01 Weight: 3.22206565472217876E-01
Root:  7.92794294644228348E-01 Weight: 6.50476830805121059E-02
#+end_src

** Notes
- The ~rec~ method fails for high values of ~beta~ so the ~gw~ method
should be used in such situations.
- ~algo665~ is an in-place variant of ~gw~ and is much faster when many points are needed
** Benchmarks
A very preliminary set can be run once all the interfaces have been compiled:
#+begin_src bash
# From $GITROOT
mv interfaces/PyInterface/gjquadpy*.so .
hyperfine --warmup 10 --export-markdown gjp_benchmarks.md \
    'bbdir/gjp_quad 5 2. 3. "gw"' \
    'bbdir/c_cli_gjpq 5 2 3 gw' \
    'PYTHONPATH=$(pwd) python interfaces/PyInterface/gjquad_cli.py --npts 5 --alpha 2 --beta 3' \
    'python scripts/scipy_gauss_jac.py --npts 5 --alpha 2 --beta 3' \
    'python scripts/sympy_gauss_jac.py --npts 5 --alpha 2 --beta 3 --n_dig 15'
#+end_src

Which gives:
#+begin_src bash
ummary
  bbdir/c_cli_gjpq 5 2 3 gw ran
    1.03 ± 0.45 times faster than bbdir/gjp_quad 5 2. 3. "gw"
   39.08 ± 11.48 times faster than PYTHONPATH=$(pwd) python interfaces/PyInterface/gjquad_cli.py --npts 5 --alpha 2 --beta 3
   55.71 ± 16.51 times faster than python scripts/scipy_gauss_jac.py --npts 5 --alpha 2 --beta 3
   86.73 ± 24.86 times faster than python scripts/sympy_gauss_jac.py --npts 5 --alpha 2 --beta 3 --n_dig 15
hyperfine --warmup 10 --export-markdown gjp_benchmarks.md       32.65s user 61.31s system 469% cpu 20.004 total
#+end_src

Or in other words:

| Command                                                     | Mean [ms]    | Min [ms] | Max [ms] | Relative      |
|-------------------------------------------------------------+--------------+----------+----------+---------------|
| ~gjp_quad 5 2. 3. "gw"~                                     | 2.6 ± 0.8    |      1.8 |      7.8 | 1.03 ± 0.45   |
| ~c_cli_gjpq 5 2 3 gw~                                       | 2.5 ± 0.7    |      1.8 |      8.4 | 1.00          |
| ~gjquad_cli.py --npts 5 --alpha 2 --beta 3~                 | 97.5 ± 6.3   |     95.4 |    130.8 | 39.08 ± 11.48 |
| ~scipy_gauss_jac.py --npts 5 --alpha 2 --beta 3~            | 139.0 ± 10.6 |    132.4 |    173.4 | 55.71 ± 16.51 |
| ~sympy_gauss_jac.py --npts 5 --alpha 2 --beta 3 --n_dig 15~ | 216.4 ± 1.7  |    214.7 |    219.7 | 86.73 ± 24.86 |

Which suggests what one might suspect, that there is a large overhead in calling
~python~ , and that the ~C~ and ~Fortran~ variants are almost exactly as fast as
each other. However, the ~f2py~ variant is still way faster than the existing
~python~ implementations.

#+begin_src bash
hyperfine --warmup 10 --export-markdown gjp_benchmarks.md \
    'PYTHONPATH=$(pwd) python interfaces/PyInterface/gjquad_cli.py --npts 5 --alpha 2 --beta 3' \
    'python scripts/scipy_gauss_jac.py --npts 5 --alpha 2 --beta 3'
Summary
  PYTHONPATH=$(pwd) python interfaces/PyInterface/gjquad_cli.py --npts 5 --alpha 2 --beta 3 ran
    1.38 ± 0.11 times faster than python scripts/scipy_gauss_jac.py --npts 5 --alpha 2 --beta 3
#+end_src
** Development
*** Developing locally
A ~pre-commit~ job is setup on CI to enforce consistent styles, so it is best to
set it up locally as well (using [[https://pypa.github.io/pipx][pipx]] for isolation):

#+begin_src sh
# Run before commiting
pipx run pre-commit run --all-files
# Or install the git hook to enforce this
pipx run pre-commit install
#+end_src
*** Updating licenses
When the headers in the sources need to be updated modify ~add_headers.py~ and run:
#+begin_src sh
python scripts/add_headers.py --dirs src/ interfaces/ --ftypes "f90,f77" --cchar '!'
python scripts/add_headers.py --dirs interfaces --ftypes "c,h" --cchar '//'
python scripts/add_headers.py --dirs interfaces scripts --ftypes "py" --cchar '#'
#+end_src
Remember to do this before exporting the code into other projects (e.g. [[https://github.com/atomic-solvers/featom][featom]]).
* License
MIT.
*** Citation
If you use this library (including the interfaces) please remember to cite it as:

#+begin_quote
Rohit Goswami. (2023). HaoZeke/GaussJacobiQuad: GaussJacobiQuad I (v0.1.0). Zenodo. https://doi.org/10.5281/ZENODO.8285112
#+end_quote

Or use the ~bibtex~ entry:
#+begin_src bibtex
@software{rgGaussQuad23,
  author       = {Rohit Goswami},
  title        = {HaoZeke/GaussJacobiQuad: GaussJacobiQuad I},
  month        = aug,
  year         = 2023,
  publisher    = {Zenodo},
  version      = {v0.1.0},
  doi          = {10.5281/zenodo.8285112},
  url          = {https://doi.org/10.5281/zenodo.8285112}
}
#+end_src
An ArXiv --> JOSS paper is in the works.

Owner

Citation (CITATION.cff)

cff-version: 1.1.0
message: "If you use this software, please cite it as below."
authors:
  - family-names: Goswami
    given-names: Rohit
    orcid: https://orcid.org/0000-0002-2393-8056
  - family-names: Čertík
    given-names: Ondřej
    orcid: https://orcid.org/0000-0003-3968-3614
title: HaoZeke/GaussJacobiQuad: GaussJacobiQuad I
version: v0.1.0
date-released: 2023-08-26

GitHub Events

Issues and Pull Requests

Last synced: 11 months ago