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.