code generation for array

[TIANJIAO SUN]

Dear all,

I'm playing around with generating some C code for numerical analysis and simulation. The library is good and easy to work with in general, only thing I'm not sure so far is how sympy is treating multi-dimensional arrays.

Simple example of what I'm looking for:

sympy code: U[i, j] = U[i, j-1] + U[i, j-2]

C code: U[i][j] = U[i][j-1] + U[i][j-2]

As you see it's pretty much direct substitution of symbols. However in order to handle >2D array, I need to replace printing.ccode._print_Indexed() with something like below:

def _print_Indexed(self, expr):
    # Array base and append indices
    output = self._print(expr.base.label) + ''.join([ '[' + self._print(x) + ']' for x in expr.indices ])
    return output

The current implementation seems to try to flatten the multi-dimensional array into 1D vector, so I will get a 1D array back and I need to supply the IndexedBase with its shape:

for i in reversed(range(expr.rank)):
    elem += expr.indices[i]*offset
    offset *= dims[i]
return "%s[%s]" % (self._print(expr.base.label), self._print(elem))

Can someone shed some light on why it is done this way? Thanks vm.

Regards,

-TS

[Björn Dahlgren]

On Monday, 1 June 2015 16:09:13 UTC+2, TIANJIAO SUN wrote:

/.../

Can someone shed some light on why it is done this way? Thanks vm.

I think it's mostly because people have been working with contiguous arrays (e.g. NumPy arrays).
Even though multidimensional arrays in C offer great flexibility, the associated pointer
indirection often makes it hard/impossible for compilers to do all their optimizations
(e.g. autovectorization etc.).

If you have a use case for pointer-to-potiner-to... type of arrays in C you may want to introduce a
keyword argument to ccode enabling this behaviour.

Note that you may want to base your work on this PR: https://github.com/sympy/sympy/pull/9314

Best,
Björn

[TIANJIAO SUN]

Many thanks.
It definitely make sense to vectorize for performance. In this case the reason I want to keep the dimensions is that I only want to generate code for the body of the loops, and create the looping separately. I feel that gives more flexibility in handling boundary conditions, and I want to retain some of the structure of the problem in the code as well.
Thanks for the pointer to the PR, tensor stuff is always exciting!
Cheers,
-TS