on., 10.11.2010 kl. 10.56 -0800, skrev Omar Awile:
> Hi all,
Hi
>
> I've recently started playing around with sympy. I am mainly
> interested in the fcode printer module as I would like to be able to
> write symbolic expressions and generate on the fly fortran code that
> is inserted into a fortran subroutine that is part of a larger
> parallel code.
>
> After trying out a couple simple examples that compiled nicely into
> fortran code I tried an example that is slightly more complicated, but
> I can't seem to figure out how to write this properly in sympy. So
> here is what I would like to express
>
> $dw_p = \frac{V_p}{e^2} \sum_{q=1}^{N}{(w_p - w_q) \mu (x_p - x_q)}
> \forall p=1,...,N$
>
> In fortran this would be something like
>
> do p=1,n
> do q=1,n
> dw(p) = dw(p) + (w(p) - w(q))*mu*(x(p)-x(q))
> enddo
> dw(p) = dw(p)*(V(p)/(e**2))
> enddo
This loop is currently too complex for the fcode printer, at least
without some hacking. Maybe we can figure something out, and this may
give us an idea on how to make the printers more sophisticated.
The exception is raised in a routine that tries to figure out which
indices (that is Idx-instances) imply summations, and which are fixed.
Repeated indices imply a summation, and this is all figured out by a
recursion into the subexpressions. A decision had to be made about how
to treat an expression like
(y[i] + x[j])*z[i]
I see two alternatives:
1) (sum_i y[i]*z[i]) + x[j]*z[i]
2) Undefined
The reason I did not choose option 1) in the current implementation, is
because this leads to the summation of objects with 0 indices and 2
indices. Since this may not be expected by the user, i preferred to err
on the cautious side, raising the exception whenever a sum contains
terms with different sets of indices. However, a long term goal is to
implement all broadcasting rules of numpy so that we can directly
produce and verify valid numpy expressions. This will give us a well
defined behaviour for option 1, so at some point we can probably disable
the exception safely.
This is the reason you see the error. Here is a relly ugly workaround
that works only because functions and arrays in Fortran look the same:
In order to avoid the index conformance requirements, you can use
functions and symbols instead of IndexedBase and Idx. Accordingly, I
modified your script to:
from sympy import IndexedBase, Idx, Symbol, Function
import sympy
w = IndexedBase('w')
dw = IndexedBase('dw')
x = IndexedBase('x')
N = Symbol('N',integer=True)
q = Idx('q',N)
p = Idx('p',N)
V = IndexedBase('V')
mu = Symbol('mu')
nu = Symbol('nu')
eps = Symbol('eps')
# expr = (w[q]-w[p])*mu*(x[q]-x[p])
wfunc = Function('w')
xfunc = Function('x')
psymb = Symbol('p')
expr = (w[q]-wfunc(psymb))*mu*(x[q]-xfunc(psymb))
sympy.printing.print_fcode(expr, assign_to=dw[p])
This gives me the loop:
do p = 1, N
do q = 1, N
dw(p) = mu*(-w(p) + w(q))*(-x(p) + x(q)) + dw(p)
end do
end do
In order to put in the V(p)/e**2, you can use the same trick, but for
optimization, you should probably do it in a separate call to fcode.
Now, the workaround is ugly, and it may be difficult to automate. So
maybe what we really need in order to deal with complex indexed
expressions, is a better way to represent summation.
Øyvind
>
> I tried to accomplish this with following python/sympy code:
> http://codepad.org/O7L70QJV
>
> But when I try to run this I get an exception:
> IndexConformanceException: Indices are not consistent: -w[p] + w[q]
>
> I must be doing something wrong, but I don't quite get what..
> (actually the code I tried doesn't even contain the $\frac{V_p}{e^2}$,
> but I don't know how to express this properly, apart from using
> sympy.sum, but the fcode printer doesn't like sum.)
>
> I was looking a bit through the documentation but wasn't sure were to
> look actually and I even bugged the devs in the IRC channel who
> adviced me to post this to the mailinglist. So I hope I'm not asking a
> too obvious question...
>
> Thanks in advance!
>
> omar
>
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to
sympy+unsubscr...@googlegroups.com.
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
