SymEngine.f90 repeatedly evaluate an expression

1

What’s the closest thing to lambdifying an expression using SymEngine.f90?

Say I have the following code:

use symengine
type(Basic) :: x, y, f
    x = Symbol("x") 
    f = parse("x**2+2*x")

I would I do the equivalent to

xx = [1.0,2.0,3.0]
do i = 1, size(xx)
   fnum = f(xx[i])
end do

So far, the best I could figure out is to parse the function, substitute the value of x and convert to numerical value inside the do loop, which can’t be efficient…

do i = 1, size(xx)
   f = parse("x**2+2*x")
   f = f%subs(x, RealDouble(xx[i]))
   fnum = f%dbl()
end do
2

I did pretty much the same thing here.
It’s probably not necessary to reparse the expression in the loop and you can call evalf to evaluate after the substitution.

real(r8) function compute(eq)
        character(*), intent(in) :: eq
        !private
        type(Basic) :: func
        integer :: j
        func = parse(eq)
        do j=1,11
            func = func%subs(var(j), val(j))
        end do
        func = func%evalf()
        compute = func%dbl()
    end function

I tested all equation parsers that I could find and Feq-parse might be more appropriate for what you want to do since it supports arrays directly.

1 Like
3

SymEngine supports fast lambdify, even using LLVM to compile, see the various supported API examples here: symengine/symengine/tests/eval/test_lambda_double.cpp at 17871adbd8366d18fc8f372f23f07e508571de46 · symengine/symengine · GitHub. So if this is what you want, we just need to expose this from Fortran.

4

@certik, that would be really great. In the mean time, instead of parsing before each evaluation, I will just copy the parsed expression.

        feval = f
        do j = 1, nvars
            feval = feval%subs(vars(j), vals(j))
        end do
        feval = feval%evalf()

This is not a very performance sensitive part of my code, so it will work fine for now.

1 Like
5

You can also try calmat (GitHub - hassaniriad/calmat: a Fortran Equation Parser Involving Matrix Operations). For your example, I can think of at least three ways to do this (depending on whether or not you are required to use a loop):

! Version 1: parse first the expresson "r = x^2 + 2*x", then evaluate in a do-loop
use calmat_m
implicit none

type   (ind_t)              :: indx, indr
integer(Ikind)              :: i, handleId
real   (Rkind), allocatable :: xx(:), res(:), r

xx = [ 1.0, 2.0, 3.0 ] ; allocate(res, mold = xx)

call calmat ( exprIn = "x = [] ; r = [] ;" ) !< a 1st call to calmat just to define in the
                                             !< calmat workspace the involved variables
                                             
call calmat_inquire ( "x", indx ) !< get the Id of the variable x 
call calmat_inquire ( "r", indr ) !< get the Id of the variable r

call calmat ( exprIn = "r = x^2 + 2*x ;", parseOnly = handleId ) !< parse the expression

do i = 1, size(xx)
   call calmat_copy ( from = xx(i), to = indx ) !< Set x to xx(i)
   call calmat ( evaluate = handleId )          !< Evaluate the exprssion
   call calmat_copy ( from = indr, to = r )     !< Get the result in the intrinsic r
   res(i) = r
end do

print*, "The result is (Version 1):"
print*,res
end

! Version 2: direct evaluation (using a for-loop in the expression)

use calmat_m
implicit none

real(Rkind), allocatable :: res(:)

call calmat ( exprIn = &
   "x = [1.,2.,3.] ; r = zeros(x) ; for i = 1:size(x) ; r(i) = x(i)^2 + 2*x(i) ; end for" ) 

call calmat_moveAlloc ( from = 'r', to = res ) !< Get the result in the intrinsic res

print*, "The result is (Version 2):"
print*,res
end


! Version 3: direct evaluation (using a vectorial form)
use calmat_m
implicit none

real(Rkind), allocatable :: res(:)

call calmat ( exprIn = "x = [1.0, 2.0, 3.0] ; r = x.^2 + 2*x ; " ) !< Note the element-wise operation

call calmat_moveAlloc ( from = 'r', to = res ) !< Get the result in the intrinsci res

print*, "The result is (Version 3):"
print*,res
end

(Note: I have tested calmat only on linux and macOS (intel) with gfortran, ifort and nagfor. If anyone wants to try it on other platforms or compilers, I’d be very interested!)

1 Like