Add Lambdify

CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 2.8)
+cmake_minimum_required(VERSION 3.10)
 enable_language(C Fortran)
 project(symengine.f90)
 set (CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/mod")

README.md

@@ -15,3 +15,53 @@ Install SymEngine.f90:
     cmake -DCMAKE_PREFIX_PATH=/some/path .
     make
     make install
+
+### Using fpm
+
+fpm does not know how SymEngine was built, so the libraries to link must be
+given with `--link-flag`: `-lsymengine`, followed, for a static SymEngine (the
+default), by the libraries it was built with (e.g. `-lgmp`, `-lLLVM`, `-lflint`,
+`-lmpfr`) and the C++ runtime. For example, for a SymEngine installed in
+`/some/path` and built with GMP and LLVM:
+
+On Linux (fpm passes `--link-flag` before the object files, so the GNU linker
+needs `--whole-archive` to keep SymEngine):
+
+    fpm test --flag "-DHAVE_SYMENGINE_LLVM" --link-flag "-L/some/path/lib \
+        -Wl,--whole-archive -lsymengine -Wl,--no-whole-archive \
+        -L$(llvm-config --libdir) -lgmp -lLLVM -lstdc++"
+
+On macOS, with GMP and LLVM from Homebrew:
+
+    fpm test --flag "-DHAVE_SYMENGINE_LLVM" --link-flag "-L/some/path/lib \
+        -L$(brew --prefix)/lib -L$(brew --prefix llvm)/lib -lsymengine -lgmp -lLLVM -lc++"
+
+`-DHAVE_SYMENGINE_LLVM` enables the LLVM backend of `Lambdify` and should only
+be passed if SymEngine was built with `WITH_LLVM=ON`.
+
+## Lambdify
+
+`Lambdify` turns expressions into numerical functions, like `Lambdify` in
+symengine.py:
+
+```fortran
+use symengine
+type(Basic) :: x, y
+type(Lambdify) :: f
+real(c_double) :: out(1), inps(2, 100), outs(1, 100)
+procedure(lambdify_function), pointer :: fn
+
+x = Symbol('x')
+y = Symbol('y')
+call f%init([ptr(x), ptr(y)], [ptr(x**2 + y**2)])   ! or arrays of type(Basic)
+call f%eval([1.0d0, 2.0d0], out)                   ! one point
+call f%eval(inps, outs)                            ! inps(2, n) -> outs(1, n)
+```
+
+Optional arguments of `init`: `cse` (common subexpression elimination),
+`backend` (`'lambda'`, the default, or `'llvm'`) and `opt_level` (LLVM only).
+With `backend='llvm'`, `call f%get_function(fn)` returns the compiled function
+as a procedure pointer, called as `call fn(inps, outs)`. The LLVM backend
+requires SymEngine built with `WITH_LLVM=ON` (with fpm, see above). A
+`Lambdify` object cannot be copied, and with the `'lambda'` backend each thread
+needs its own object.

src/CMakeLists.txt

@@ -9,9 +9,13 @@ set(SRC
     constants.f90
     dense_matrix.f90
     sets.f90
+    lambdify.F90
     symengine.f90
     )
 add_library(symengine_f90 ${SRC})
 target_link_libraries(symengine_f90 ${SYMENGINE_LIBRARIES})
+if (HAVE_SYMENGINE_LLVM)
+    target_compile_definitions(symengine_f90 PRIVATE HAVE_SYMENGINE_LLVM)
+endif()
 
 install(TARGETS symengine_f90 DESTINATION lib)

src/lambdify.F90

@@ -0,0 +1,276 @@
+module symengine_lambdify
+
+! Turns symbolic expressions into numerical functions, like Lambdify in
+! symengine.py:
+!
+!     type(Lambdify) :: f
+!     call f%init([ptr(x), ptr(y)], [ptr(x**2 + y**2)])
+!     call f%eval([1.0d0, 2.0d0], outs)          ! one point
+!     call f%eval(inps, outs)                    ! inps(2, n), outs(1, n)
+!
+! Two backends are available:
+!  - 'lambda' (default): a tree of closures, always available.
+!  - 'llvm': machine code generated by the LLVM JIT, only if SymEngine was
+!    built with LLVM. get_function returns it as a procedure pointer.
+!
+! A Lambdify object owns a C++ object and cannot be copied. With the 'lambda'
+! backend, use one object per thread.
+
+use iso_c_binding, only: c_int, c_long, c_double, c_size_t, c_ptr, c_null_ptr, c_funptr, &
+                         c_associated, c_f_procpointer
+use symengine_interface
+use symengine_basic
+use exceptions
+
+implicit none
+
+integer, parameter :: BACKEND_NONE = 0, BACKEND_LAMBDA = 1, BACKEND_LLVM = 2
+
+type Lambdify
+    type(c_ptr), private :: ptr = c_null_ptr
+    integer, private :: backend = BACKEND_NONE
+contains
+    procedure :: lambdify_init_ptr, lambdify_init_basic
+    generic :: init => lambdify_init_ptr, lambdify_init_basic
+    procedure :: num_inputs => lambdify_num_inputs
+    procedure :: num_outputs => lambdify_num_outputs
+    procedure :: lambdify_eval_point, lambdify_eval_batch
+    generic :: eval => lambdify_eval_point, lambdify_eval_batch
+    procedure :: get_function => lambdify_get_function
+    procedure :: free => lambdify_free
+    procedure, private :: lambdify_assign
+    generic :: assignment(=) => lambdify_assign
+    final :: lambdify_final
+end type
+
+! Interface of the function returned by get_function
+abstract interface
+    subroutine lambdify_function(inps, outs) bind(c)
+        import :: c_double
+        real(c_double), dimension(*), intent(in) :: inps
+        real(c_double), dimension(*), intent(out) :: outs
+    end subroutine
+end interface
+
+private
+public :: Lambdify, lambdify_function
+
+contains
+
+subroutine lambdify_init_ptr(this, args, exprs, cse, backend, opt_level)
+    class(Lambdify), intent(inout) :: this
+    type(c_ptr), dimension(:), intent(in) :: args, exprs
+    logical, intent(in), optional :: cse
+    character(len=*), intent(in), optional :: backend
+    integer, intent(in), optional :: opt_level
+
+    type(c_ptr) :: vargs, vexprs
+    integer(c_int) :: perform_cse, level, exception
+    integer :: i
+
+    call this%free()
+
+    exception = 0
+    perform_cse = 0
+    if (present(cse)) then
+        if (cse) perform_cse = 1
+    end if
+    level = 3
+    if (present(opt_level)) level = int(opt_level, c_int)
+
+    this%backend = BACKEND_LAMBDA
+    if (present(backend)) then
+        select case (backend)
+        case ('lambda')
+            this%backend = BACKEND_LAMBDA
+        case ('llvm')
+#ifdef HAVE_SYMENGINE_LLVM
+            this%backend = BACKEND_LLVM
+#else
+            error stop "Lambdify: SymEngine was built without LLVM support"
+#endif
+        case default
+            error stop "Lambdify: backend must be 'lambda' or 'llvm'"
+        end select
+    end if
+
+    vargs = c_vecbasic_new()
+    vexprs = c_vecbasic_new()
+    do i = 1, size(args)
+        call handle_exception(c_vecbasic_push_back(vargs, args(i)))
+    end do
+    do i = 1, size(exprs)
+        call handle_exception(c_vecbasic_push_back(vexprs, exprs(i)))
+    end do
+
+    select case (this%backend)
+    case (BACKEND_LAMBDA)
+        this%ptr = c_lambda_real_double_visitor_new()
+        exception = c_lambda_real_double_visitor_init(this%ptr, vargs, vexprs, perform_cse)
+#ifdef HAVE_SYMENGINE_LLVM
+    case (BACKEND_LLVM)
+        this%ptr = c_llvm_double_visitor_new()
+        exception = c_llvm_double_visitor_init(this%ptr, vargs, vexprs, perform_cse, level)
+#endif
+    end select
+
+    call c_vecbasic_free(vargs)
+    call c_vecbasic_free(vexprs)
+    if (exception /= 0) then
+        call this%free()
+        call handle_exception(int(exception, c_long))
+    end if
+end subroutine
+
+subroutine lambdify_init_basic(this, args, exprs, cse, backend, opt_level)
+    class(Lambdify), intent(inout) :: this
+    class(Basic), dimension(:), intent(in) :: args, exprs
+    logical, intent(in), optional :: cse
+    character(len=*), intent(in), optional :: backend
+    integer, intent(in), optional :: opt_level
+
+    type(c_ptr), dimension(size(args)) :: pargs
+    type(c_ptr), dimension(size(exprs)) :: pexprs
+    integer :: i
+
+    do i = 1, size(args)
+        pargs(i) = args(i)%ptr
+    end do
+    do i = 1, size(exprs)
+        pexprs(i) = exprs(i)%ptr
+    end do
+    call this%lambdify_init_ptr(pargs, pexprs, cse, backend, opt_level)
+end subroutine
+
+function lambdify_num_inputs(this) result(res)
+    class(Lambdify), intent(in) :: this
+    integer :: res
+
+    res = 0
+    select case (this%backend)
+    case (BACKEND_LAMBDA)
+        res = int(c_lambda_real_double_visitor_get_num_inputs(this%ptr))
+#ifdef HAVE_SYMENGINE_LLVM
+    case (BACKEND_LLVM)
+        res = int(c_llvm_double_visitor_get_num_inputs(this%ptr))
+#endif
+    end select
+end function
+
+function lambdify_num_outputs(this) result(res)
+    class(Lambdify), intent(in) :: this
+    integer :: res
+
+    res = 0
+    select case (this%backend)
+    case (BACKEND_LAMBDA)
+        res = int(c_lambda_real_double_visitor_get_num_outputs(this%ptr))
+#ifdef HAVE_SYMENGINE_LLVM
+    case (BACKEND_LLVM)
+        res = int(c_llvm_double_visitor_get_num_outputs(this%ptr))
+#endif
+    end select
+end function
+
+! Evaluates at npts points. inps and outs are stored point after point.
+subroutine lambdify_call(this, inps, outs, npts)
+    class(Lambdify), intent(inout) :: this
+    real(c_double), dimension(*), intent(in) :: inps
+    real(c_double), dimension(*), intent(out) :: outs
+    integer, intent(in) :: npts
+
+    select case (this%backend)
+    case (BACKEND_LAMBDA)
+        call c_lambda_real_double_visitor_call_batch(this%ptr, outs, inps, int(npts, c_size_t))
+#ifdef HAVE_SYMENGINE_LLVM
+    case (BACKEND_LLVM)
+        call c_llvm_double_visitor_call_batch(this%ptr, outs, inps, int(npts, c_size_t))
+#endif
+    case default
+        error stop "Lambdify: eval called before init"
+    end select
+end subroutine
+
+! Evaluates at one point: inps(num_inputs), outs(num_outputs)
+subroutine lambdify_eval_point(this, inps, outs)
+    class(Lambdify), intent(inout) :: this
+    real(c_double), dimension(:), intent(in) :: inps
+    real(c_double), dimension(:), intent(out) :: outs
+
+    integer :: nin, nout
+
+    if (this%backend == BACKEND_NONE) error stop "Lambdify: eval called before init"
+    nin = this%num_inputs()
+    nout = this%num_outputs()
+    if (size(inps) /= nin .or. size(outs) /= nout) then
+        error stop "Lambdify: eval called with arrays of the wrong size"
+    end if
+    call lambdify_call(this, inps, outs, 1)
+end subroutine
+
+! Evaluates at several points: inps(num_inputs, npts), outs(num_outputs, npts)
+subroutine lambdify_eval_batch(this, inps, outs)
+    class(Lambdify), intent(inout) :: this
+    real(c_double), dimension(:, :), intent(in) :: inps
+    real(c_double), dimension(:, :), intent(out) :: outs
+
+    integer :: nin, nout
+
+    if (this%backend == BACKEND_NONE) error stop "Lambdify: eval called before init"
+    nin = this%num_inputs()
+    nout = this%num_outputs()
+    if (size(inps, 1) /= nin .or. size(outs, 1) /= nout .or. size(inps, 2) /= size(outs, 2)) then
+        error stop "Lambdify: eval called with arrays of the wrong shape"
+    end if
+    call lambdify_call(this, inps, outs, size(inps, 2))
+end subroutine
+
+! Returns the JIT compiled function (backend='llvm' only). It can be called
+! concurrently from several threads, and is valid as long as this object is
+! alive and not re-initialized.
+subroutine lambdify_get_function(this, fn)
+    class(Lambdify), intent(in) :: this
+    procedure(lambdify_function), pointer, intent(out) :: fn
+#ifdef HAVE_SYMENGINE_LLVM
+    type(c_funptr) :: cfn
+#endif
+
+    nullify(fn)
+    if (this%backend /= BACKEND_LLVM) then
+        error stop "Lambdify: get_function requires backend='llvm'"
+    end if
+#ifdef HAVE_SYMENGINE_LLVM
+    cfn = c_llvm_double_visitor_get_function(this%ptr)
+    call c_f_procpointer(cfn, fn)
+#endif
+end subroutine
+
+subroutine lambdify_free(this)
+    class(Lambdify), intent(inout) :: this
+
+    if (c_associated(this%ptr)) then
+        select case (this%backend)
+        case (BACKEND_LAMBDA)
+            call c_lambda_real_double_visitor_free(this%ptr)
+#ifdef HAVE_SYMENGINE_LLVM
+        case (BACKEND_LLVM)
+            call c_llvm_double_visitor_free(this%ptr)
+#endif
+        end select
+    end if
+    this%ptr = c_null_ptr
+    this%backend = BACKEND_NONE
+end subroutine
+
+subroutine lambdify_final(this)
+    type(Lambdify), intent(inout) :: this
+    call this%free()
+end subroutine
+
+subroutine lambdify_assign(a, b)
+    class(Lambdify), intent(inout) :: a
+    class(Lambdify), intent(in) :: b
+    error stop "Lambdify: objects cannot be copied, call init instead"
+end subroutine
+
+end module

src/symengine.f90

@@ -11,6 +11,7 @@ module symengine
 use constants
 use dense_matrix
 use sets
+use symengine_lambdify
 
 implicit none
 
@@ -27,6 +28,7 @@ module symengine
 public :: DenseMatrix, transpose, ones, zeros, eye
 public :: emptyset, universalset, complexes, reals, rationals, integers, set_union, set_intersection
 public :: interval, finiteset
+public :: Lambdify, lambdify_function
 
 contains