In my previous blog post I wrote about the new layer of testing that became available when the tensor module got functionality to analyze and report the structure of indexed expressions. Since then, I have improved the C printer and the C code generator by removing hacks that used to provide such functionality and replaced it with calls to the more robust methods defined in the tensor module.

I put up these improvements for review, and got valuable feedback that I spent the rest of the time implementing. In addition, I have addressed Ronan’s comments to my last blog post, and this has lead to stylistic improvements:

- Now an indexed expression is written with [] instead of (). This is more pythonic, as [] is used with all the builtin sequences. This is also consistent with numpy.
- Some classes have been renamed: IndexedElement is now Indexed, and the old Indexed class is now called IndexedBase. The new terminology should be less confusing.

The C code generator has become quite versatile and capable of handling complicated array expressions. Here is a little demonstration of my current codegen_C3 branch:

In [1]: A, x, y = map(IndexedBase, ['A', 'x', 'y'])
In [2]: i = Idx('i', m)
In [3]: j = Idx('j', n)
In [4]: print ccode(A[i, i], assign_to=z) #trace of a matrix
z = 0;
for (int i=0; i<m; i++){
z = z + A[i + i*m];
}
In [5]: print ccode(A[i, j]*x[j] + y[i], assign_to=y[i]) #matrix-vector 1
for (int i=0; i<m; i++){
for (int j=0; j<n; j++){
y[i] = x[j]*A[j + i*n] + y[i];
}
}
In [6]: print ccode(A[i, j]*x[j], assign_to=y[i]) #matrix-vector 2
for (int i=0; i<m; i++){
y[i] = 0;
}
for (int i=0; i<m; i++){
for (int j=0; j<n; j++){
y[i] = x[j]*A[j + i*n] + y[i];
}
}

The code printer takes care of introducing and initialization of accumulator variables if needed, not the codegen module as previaously. The initialization is skipped if the left hand side variable is also present on the right hand side (as in statement In[5]).

Finally, you may have noticed that this post was a little bit delayed. This is because I was traveling in the weekend and didn’t get as much time for blogging as I had planned. I’ll make up for it by blogging twice this week.

### Like this:

Like Loading...

*Related*

I’m curious about how your Fortran and C code generators work and what applications they would be useful in. I’ve written a small DSL in scala for index notation (for solid mechanics) and I was going to write a simple code generator for that, but it seems like you have a pretty good setup in python already. Is your code merged into the official numpy branch, or maintained only in your branch?

It has been included in SymPy. If you want to generate loops over arrays with SymPy, you must use the Tensor Module, which is currently only available in the git repository. You may also be interested in the codegen and autowrap utilities (also included with SymPy).

Thanks for quick reply. I noticed in this particular post you are using ccode. I have been playing around with it and I have noticed some function seem to get printed even though they are not in the “known_functions” list.

i.e.

print ccode(A[i,j]*x[j]*exp(3), assign_to=y[i])

I was wondering why that works.

Also I noticed Matrices cannot be printed. I am somewhat interested in supporting this so it would be nice if you could give me some hints for adding that in. For example I would like to do something like the following.

double xi[] = {1,2,3,1,1,2};

double xj[] = {1,2,3,2,3,3};

for( l= 0;l < 6; l++)

for( j = 0; j < 3 ; j++)

C[l] += F[j][xi[l]]*F[j][xj[l]];

Where F is 3×3 matrix and C is symmetric and can be reduced to 6 entries using voigt notation. Likewise, calculating CC = C[i,j]*C[i,j] in your python notation can be done on voigt represented matrices as so:

double vg[3][3] = {{1,4,5},{4,2,6},{5,6,3}};

for ( i = 0;i<3; i++)

for( j=0; j<3 ; j++)

CC += C[vg[i][j]] * C[vg[i][j]];

As you are the author of ccode(?), what are your thoughts on modifying or adding to ccode such that ccode could generate what is shown above should you happen to know that C is a symmetric matrix apriori.

Perhaps the known_functions dict got an unfortunate name. It is actually a mapping of common functions that are known to have different names in C and Sympy, such as ceiling() and ceil(). Functions that are not in the dict are printed without mapping.

The functionality you describe would be extremely cool! Perhaps it could be triggered by a keyword option, a la

`symmetric=True`

, to the`IndexedBase`

object? I think that would be the appropriate place to store such an assumption. Then hopefully`get_contraction_structure()`

in the sympy.tensor.index_methods module can be improved so that it detects and reports contractions where such optimizations are possible. If this can be done, I believe it will be easy to implement the actual code printing. If you decide to go for it, it is a good idea to open a feature request on the issue tracker, to serve as a central place to discuss design ideas.PS! I am not the author of the entire ccode, I just worked on it this summer to implement array printing, and some other stuff.

Oh, and about

`double vg[3][3] = {{1,4,5},{4,2,6},{5,6,3}};`

I think this can be done by just implementing ccode._print_Matrix(). It will be called by ccode() whenever a sympy.matrices.matrices.Matrix object is encountered.