Feedback wanted on library design

1

I’ve written a Python C API to the molecular geometry optimization package DL-FIND. The design of DL-FIND requires the user of the package to supply a number of functions at compile time that are later called by the library during runs: dlf_error, dlf_get_gradient, dlf_get_hessian and so on. This design is obviously not compatible with a Python library where the user could dynamically provide different versions of these functions at runtime. My solution to this was to construct dummy functions on the Fortran side that just pass on their arguments to a similarly named procedure imported from the module mod_globals.

subroutine dlf_get_gradient(nvar, coords, energy, gradient, iimage, kiter, status)
  use mod_globals, only: dlf_get_gradient_
  use dlf_parameter_module, only: rk

  implicit none
  integer, intent(in) :: nvar ! number of xyz variables (3*nat)
  real(rk), intent(in) :: coords(nvar) ! coordinates
  real(rk), intent(out) :: energy ! energy
  real(rk), intent(out) :: gradient(nvar) ! gradient
  integer, intent(in) :: iimage ! current image (for NEB)
  integer, intent(in) :: kiter ! flag related to microiterations
  integer, intent(out) :: status ! return code

  call dlf_get_gradient_(nvar, coords, energy, gradient, iimage, kiter, status)
end subroutine

mod_globals stores all the procedures as procedure pointers to the functions obtained from the C (Python) side.

...
procedure(c_dlf_get_gradient), pointer :: dlf_get_gradient_ => null()
...

The abstract interfaces for these functions are defined in the mod_api module

  abstract interface
    subroutine c_dlf_get_gradient(nvar, coords, energy, gradient, iimage, kiter, status)
      import c_double, c_int

      implicit none
      integer(c_int), intent(in), value :: nvar
      real(c_double), intent(in) :: coords(nvar)
      real(c_double), intent(out) :: energy
      real(c_double), intent(out) :: gradient(nvar)
      integer(c_int), intent(in), value :: iimage
      integer(c_int), intent(in), value :: kiter
      integer(c_int), intent(out) :: status
    end subroutine
  end interface

Finally, the way the procedure pointers in mod_globals get associated to the functions from the C side is by arguments as c_funptr to the function api_dl_find that is called from the C/Python side. c_f_procpointer is used to convert from C pointers to Fortran pointers.

subroutine api_dl_find(... , c_dlf_get_gradient_, ...) bind(c)
...
  type(c_funptr), intent(in), value :: ..., c_dlf_get_gradient_, .... ! Functions received from C side
...
 call c_f_procpointer(c_dlf_get_gradient_, dlf_get_gradient_)
...

This design for the API seems to work well in my tests with different energy and gradient calculators, but also feels rather hacky. If anyone has any ideas for how to improve or simplify it would be greatly appreciated.

I’m building Python wheels on GitHub Actions, linking statically on Windows and “repairing” the wheels on Linux and MacOS to include the required libraries. Once I make an initial release I will continue with conda-forge packages as well.

Some references:

DL-FIND is also available via Py-ChemShell, but I wanted something that is:

  • Leaner and more modular
  • Compatible with arbitrary energy and gradient codes
  • Available to install as a dependency via package managers
1 Like
2

There are certainly some discussions around about designing interfaces to optimization libraries, best checkout the thread at

A central topic is the design and separation of the objective function callback and the solver object. I collected a couple of links at Callback function design · Issue #16 · grimme-lab/nlopt-f · GitHub.

2 Likes
3

A quick question: does your actual code include bind(C) clause for the interfaces? e.g., for c_dlf_get_gradient?

1 Like
4

I’d be very cautious placing the pointers in a global scope. This means your Python library can’t be used in a parallel context (assuming that DL-FIND is thread-safe to begin with…). You should document the “parallel status” explicitly in your README.

1 Like
5

As opposed to data pointers, how do procedure pointers affect thread safety and/or using the library in a parallel context?

1 Like
6

It appears my understanding of how procedure pointers work is incomplete. Wouldn’t two or more threads trying to change the same procedure pointer potentially lead to a race?

The example I built to test this, shows no wrong behavior using a global pointer. (I only inspected small n sizes, and didn’t attempt to verify that multiple threads are actually being used.)

(No) Race

Compile with: $ gfortran -Wall -fopenmp -O2 polyrace.f90

module poly

implicit none
private

public :: the_poly, poly1, poly2, poly3, poly4, wp

integer, parameter :: wp = kind(1.0d0)

abstract interface
  function poly_interface(x) result(y)
    import wp
    real(wp), intent(in) :: x
    real(wp) :: y
  end function
end interface

procedure(poly_interface), pointer :: the_poly => null()

contains

  function poly1(x) result(y)
    real(wp), intent(in) :: x
    real(wp) :: y
    print *, "calling poly1 with x =", x
    y = 1.0_wp + x
  end function

  function poly2(x) result(y)
    real(wp), intent(in) :: x
    real(wp) :: y
    print *, "calling poly2 with x =", x
    y = 1.0_wp + x*(1.0_wp + x)
  end function

  function poly3(x) result(y)
    real(wp), intent(in) :: x
    real(wp) :: y
    print *, "calling poly3 with x =", x
    y = 1.0_wp + x*(1.0_wp + x*(1.0_wp + x))
  end function

  function poly4(x) result(y)
    real(wp), intent(in) :: x
    real(wp) :: y
    print *, "calling poly4 with x =", x
    y = 1.0_wp + x*(1.0_wp + x*(1.0_wp + x*(1.0_wp + x)))
  end function

end module

program polyrace

  use poly 
  implicit none

  integer, parameter :: n = 1000

  real(wp) :: y(n), r
  integer :: i 
  
  !$omp parallel do schedule(static) shared(the_poly,y)
  do i = 1, size(y)

    call random_number(r)
    if (r < 0.25_wp) then
      the_poly => poly1
    else if (r < 0.5_wp) then
      the_poly => poly2
    else if (r < 0.75_wp) then
      the_poly => poly3
    else
      the_poly => poly4
    end if 

    y(i) = the_poly(r)

  end do
  !$omp end parallel do

end program
1 Like
7

It does not. After some experimentation, compilation seemed to crash when I included it, if I recall correctly. I’m not sure if it’s because I am using the c_f_procpointer method. The GCC documentation that I followed omits bind(c).

8

The problem is that DL-FIND expects to find the callback functions in the global scope. I don’t really want to change the design of DL-FIND itself (to be able to upgrade to the latest version with minimal changes) and storing the pointers in a module scope was the only solution that I could come up with where they would be accessible from both (a) the function called from the C side and (b) the dummy compiled callback. I will investigate the thread safety, but it could potentially work well with multiprocessing where each Python process would get its own copy of the library. DL-FIND has internal parallelization via MPI which I have not investigated.

9

Thanks everyone who gave me feedback on this project. As a result I made a number of changes that I think improved the library a lot:

  • Added bind(c) to the abstract interfaces on the suggestion of @FortranFan. This seems to be an omission in the example in the GCC documentation.
  • Reworked most of the factory functions into decorators on the suggestion of @ivanpribec
  • Clarified the thread safety in the README.md on the suggestion of @ivanpribec. Seems to crash with multiple threads but runs fine with multiple processes.

I opted not to make any greater changes to the design of the API to be as close as possible to DL-FIND itself. This saves me from having to write a detailed documentation as I can just refer to the documentation of DL-FIND.