6. C and Python bindings

quad gets a Gauss-Legendre module: the nodes and weights of the n-point rule, exact for polynomials of degree 2n - 1:

!< Gauss-Legendre quadrature: nodes and weights on [-1, 1].
module gauss_legendre
!< Gauss-Legendre quadrature: nodes and weights on [-1, 1].
use quad_rules, only : R8P
implicit none
private
public :: gauss_nodes

contains
   pure subroutine gauss_nodes(n, x, w)
   !< Nodes and weights of the n-point rule: the roots of the Legendre polynomial P_n, by Newton's method.
   integer,   intent(in)  :: n    !< Number of points.
   real(R8P), intent(out) :: x(n) !< Nodes, ascending.
   real(R8P), intent(out) :: w(n) !< Weights.
   real(R8P), parameter   :: pi = acos(-1._R8P) !< Pi.
   real(R8P)              :: z    !< Root estimate.
   real(R8P)              :: p    !< P_n(z).
   real(R8P)              :: dp   !< P_n'(z).
   integer                :: i, step !< Counters.

   do i = 1, (n + 1) / 2
      z = cos(pi * (i - 0.25_R8P) / (n + 0.5_R8P))
      do step = 1, 100
         call legendre(n, z, p, dp)
         z = z - p / dp
         if (abs(p / dp) < 1.e-15_R8P) exit
      enddo
      call legendre(n, z, p, dp) ! the derivative at the root itself, for the weight
      x(i) = -z
      x(n + 1 - i) = z
      w(i) = 2._R8P / ((1._R8P - z * z) * dp * dp)
      w(n + 1 - i) = w(i)
   enddo
   endsubroutine gauss_nodes

   pure subroutine legendre(n, z, p, dp)
   !< P_n(z) and its derivative, by the three-term recurrence.
   integer,   intent(in)  :: n  !< Degree.
   real(R8P), intent(in)  :: z  !< Abscissa.
   real(R8P), intent(out) :: p  !< P_n(z).
   real(R8P), intent(out) :: dp !< P_n'(z).
   real(R8P)              :: p1, p2 !< P_(j-1), P_(j-2).
   integer                :: j  !< Counter.

   p = 1._R8P
   p1 = 0._R8P
   do j = 1, n
      p2 = p1
      p1 = p
      p = ((2 * j - 1) * z * p1 - (j - 1) * p2) / j
   enddo
   dp = n * (z * p - p1) / (z * z - 1._R8P)
   endsubroutine legendre
endmodule gauss_legendre

To use it from C or Python, its procedures need a C interface: bind(c), C types, no assumed shapes. fobis bindings writes it, without touching the module:

$ fobis bindings quad_rules gauss_legendre
1 of 4 procedure(s) bound
  bound    gauss_legendre::gauss_nodes -> gauss_legendre_gauss_nodes
  skipped  quad_rules::midpoint: argument f: procedure arguments are not bound
  skipped  quad_rules::trapezoid: argument f: procedure arguments are not bound
  skipped  quad_rules::simpson: argument f: procedure arguments are not bound
written: bindings/quad_c.f90, bindings/quad.h, bindings/quad.py

The rules of quad_rules take a function as argument: that has no C counterpart here, and the report says why they are skipped. Three files are written: a Fortran shim module, a C header, a Python module.

quad_c.f90

$ cat bindings/quad_c.f90
module quad_c
!< bind(c) shims generated by fobis bindings: do not edit, regenerate.
use, intrinsic :: iso_c_binding
use gauss_legendre
use quad_rules, only : R8P
implicit none
private
public :: gauss_legendre_gauss_nodes
contains

subroutine gauss_legendre_gauss_nodes(n, x, w) bind(c, name="gauss_legendre_gauss_nodes")
  integer(c_int32_t), intent(in) :: n
  integer :: n_f
  real(c_double), intent(out) :: x(n)
  real(r8p) :: x_f(n)
  real(c_double), intent(out) :: w(n)
  real(r8p) :: w_f(n)
  n_f = int(n, kind(1))
  call gauss_nodes(n=n_f, x=x_f, w=w_f)
  x = real(x_f, c_double)
  w = real(w_f, c_double)
end subroutine gauss_legendre_gauss_nodes

end module quad_c

quad.h

$ cat bindings/quad.h
/* C interface of quad: generated by fobis bindings, do not edit. Link with the library built from the
   project's sources and quad_c.f90. Arrays are Fortran arrays: contiguous, column-major. */
#ifndef QUAD_H
#define QUAD_H
#include <stdbool.h>
#include <stdint.h>

#ifdef __cplusplus
extern "C" {
#endif

/* gauss_legendre::gauss_nodes(n intent(in), x(n) intent(out), w(n) intent(out)) */
void gauss_legendre_gauss_nodes(const int32_t *n, double *x, double *w);

#ifdef __cplusplus
}
#endif
#endif /* QUAD_H */

The shim converts each argument to and from its C type explicitly (R8P is a C double because it has 15 digits, not because of its name) and repeats the module's use statements, which give the kinds. A shared library of the shim and the modules it uses:

[python]
help      = The C and Python bindings, a shared library
template  = template-gnu
cflags    = -c -O2 -fPIC
src       = ./src/ ./bindings/
build_dir = ./build/python/
target    = bindings/quad_c.f90
mklib     = shared
output    = libquad_c.so

$ fobis build --mode python --quiet

Building src/romberg.f
Compiling src/romberg.f serially
Building bindings/quad_c.f90
Compiling src/quad_rules.F90 serially
Compiling src/gauss_legendre.F90 serially
Compiling bindings/quad_c.f90 serially
Linking build/python/libquad_c.so
Target bindings/quad_c.f90 has been successfully built

And from Python:

"""The 3-point Gauss-Legendre rule, from Fortran, through the generated ctypes wrapper."""

import sys

sys.path.insert(0, "bindings")
import quad  # bindings/quad.py, written by fobis bindings

quad.load("build/python/libquad_c.so")
x, w = quad.gauss_nodes(3, [0.0] * 3, [0.0] * 3)
print("nodes  ", " ".join(f"{v:+.15f}" for v in x))
print("weights", " ".join(f"{v:.15f}" for v in w))
# exact up to degree 2n - 1 = 5: the integral of x**4 on [-1, 1] is 2/5
print("integral of x**4:", sum(wi * xi**4 for xi, wi in zip(x, w, strict=True)))

$ python3 gauss.py
nodes   -0.774596669241483 +0.000000000000000 +0.774596669241483
weights 0.555555555555555 0.888888888888889 0.555555555555555
integral of x**4: 0.3999999999999998

The integral of x⁴ on [-1, 1] is 2/5: three points integrate it exactly, up to rounding.

What you learned

fobis bindings MODULE..., what can be bound and what cannot, building and calling the result. Reference: C and Python bindings.

Next: 7. Instructions for agents.