Hi,
I finally needed to generate Fortran code from a SymPy expression. I
couldn't find any documentation for it. Is there any? I found the
FCodeGen object, but
it has no useful docstring.... So I looked at tests and those are
good. Here is my real life code that I needed today:
from StringIO import StringIO
from sympy.utilities.codegen import FCodeGen, Routine
from sympy import var, legendre, integrate, exp
var("l r12 r1 r2 D")
f = (2*l+1) / (2*r1*r2) * integrate(exp(-r12/D)*legendre(l,
(r1**2-r12**2+r2**2) / (2*r1*r2)), (r12, r1-r2, r1+r2))
routines = []
for _l in range(3):
expr = f.subs(l, _l).doit().simplify()
routines.append(Routine("V%d" % _l, expr))
code_gen = FCodeGen()
output = StringIO()
code_gen.dump_f95(routines, output, "file", False, False)
source = output.getvalue()
print source
This prints:
REAL*8 function V0(D, r1, r2)
implicit none
REAL*8, intent(in) :: D
REAL*8, intent(in) :: r1
REAL*8, intent(in) :: r2
V0 = (-D*exp(-(r1 + r2)/D) + D*exp(-(r1 - r2)/D))/(2*r1*r2)
end function
REAL*8 function V1(D, r1, r2)
implicit none
REAL*8, intent(in) :: D
REAL*8, intent(in) :: r1
REAL*8, intent(in) :: r2
V1 = 3*D*(-D**2*exp(2*r2/D) + D**2 - D*r1*exp(2*r2/D) + D*r1 + D*r2*exp( &
2*r2/D) + D*r2 + r1*r2*exp(2*r2/D) + r1*r2)*exp(-r1/D - r2/D)/(2* &
r1**2*r2**2)
end function
REAL*8 function V2(D, r1, r2)
implicit none
REAL*8, intent(in) :: D
REAL*8, intent(in) :: r1
REAL*8, intent(in) :: r2
V2 = 5*D*(9*D**4*exp(2*r2/D) - 9*D**4 + 9*D**3*r1*exp(2*r2/D) - 9*D**3* &
r1 - 9*D**3*r2*exp(2*r2/D) - 9*D**3*r2 + 3*D**2*r1**2*exp(2*r2/D &
) - 3*D**2*r1**2 - 9*D**2*r1*r2*exp(2*r2/D) - 9*D**2*r1*r2 + 3*D &
**2*r2**2*exp(2*r2/D) - 3*D**2*r2**2 - 3*D*r1**2*r2*exp(2*r2/D) - &
3*D*r1**2*r2 + 3*D*r1*r2**2*exp(2*r2/D) - 3*D*r1*r2**2 + r1**2*r2 &
**2*exp(2*r2/D) - r1**2*r2**2)*exp(-r1/D - r2/D)/(2*r1**3*r2**3)
end function
It looks pretty good. I would like to use "real(dp)" instead of
"REAL*8", as well as probably create the whole module and just put
"implicit none" there just one time. But otherwise it looks good, it
compiles and does the job. Nice!
I just found the documentation here:
http://docs.sympy.org/dev/modules/utilities/codegen.html
So I think I should use the codegen() function. Here is how to use it
to do the same as above:
from sympy.utilities.codegen import codegen
from sympy import var, legendre, integrate, exp
var("l r12 r1 r2 D")
f = (2*l+1) / (2*r1*r2) * integrate(exp(-r12/D)*legendre(l,
(r1**2-r12**2+r2**2) / (2*r1*r2)), (r12, r1-r2, r1+r2))
routines = []
for _l in range(3):
expr = f.subs(l, _l).doit().simplify()
routines.append(("V%d" % _l, expr))
[(f_name, f_code), (interface_name, f_interface)] = codegen(routines,
"F95", "potential", header=False)
print f_code
The f_interface is the f95 interface to the functions above. I don't
use such things in Fortran, I just create a module and then import it
from the module using "use moduly, only: V0". See here for more
details:
https://github.com/qsnake/qsnake/wiki/Fortran-best-practices
But the f_code does what it should. So I think I'll add the
functionality in there, as an optional parameter, to create a regular
Fortran module.
Ondrej
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
