As I mentioned in my previous post, I did not want to implement array arguments for the C code generator by simply taking the same route as I did with the Fortran generator. In hindsight it was clear that I could do a much better implementation, and I’ve spent the week figuring out exactly how it should be done, and putting together a working implementation. By the way it’s been beneficial to do the C and Fortran implementations in parallel, as it provides a good opportunity to test out alternative implementations.

The principal weakness of the Fortran code generator is that the support of array arguments has only limited testability. It is currently testable only on two levels:

1) The properties and behavior of IndexedElement instances,

2) The output from the Fortran code generator.

The missing piece in there, is the testability of complex expressions composed of IndexedElement instances. It is important to be able to test the consistency of indices in an expression, as it is very easy to create ambigous or even meaningless expessions of IndexedElement objects. For instance, the innocent looking expression: (which is considered undefined in my development branch)

x(i) + x(j)

would be ambigous at best. It could be (1) a scalar being the sum of components i and j of array x or (2) a symmetric rank 2 object a(i, j) where every component is made by adding components i and j of x. To avoid this kind of ambiguities, it is necessary to implement checks of index consistency. And it is preferable to implement it in the indexed module so that these checks are available whenever the module is in use. This gives the bonus of another level of testability:

3) The contraction structure of indexed expressions and the algorithms used to analyze the expressions.

With this in place, the reliability of array based code from the C code generator depends directly on how robust the analysis of indexed expressions is. The third level of testability is extremely important to make Sympy generated code a viable option for scientific applications.

In my branch codegen_C on github, these thoughts have materialized in a new file sympy/tensor/indexed_methods.py, as well as in modifications of the C code printer. I have implemented testable high-level interfaces that can analyze expressions and provide information about the index structure. The code printers should now be able to generate loops over arrays in a much more reliable fashion.

### Like this:

Like Loading...

*Related*

I’m rather confused by your terminology. The fact that you’re apparently using parentheses instead of brackets to denote indexation doesn’t help.

In particular, I don’t understand what IndexedElement really means. Obviously, it’s not an element of a tensor, otherwise an expression like a[i] + a[j] wouldn’t cause any trouble. I think it behaves like a tensor with named axes but, if that’s the case, I wonder how to access the actual elements of the tensor.

Thanks for your comment, you raise two very good points. Brackets would be more pythonic, and it would be a good idea to change that.

I also see your point about IndexedElement not really being an element, and the notion of “named axes” is accurate. Perhaps it would be better to just call it Indexed, and denote the stem as IndexedStem.

About the actual elements, the idea is that if all indices are “fixed”, it represents a scalar element. But I am not sure yet how to represent a fixed index. One way is to have a separate class FixedIndex, and another is to use an Idx with range (None, None). The latter is tempting, but I also think it could be useful to interpret that as just an unspecified range.

edit: A Symbol or Integer is also a good choice to represent a fixed index.