Symbolic matrices in sympy, print the trace and generate C code

user_8458838 0
Viewed 872

I am new to sympy an would like to print the trace of a symbolic matrix well as generate C code with the function ccode

Currently, I have the following:

import sympy as sp

# Creates a symbolic matrix (3x3) with symbol A
A = sp.MatrixSymbol('A', 3, 3)

# Trace of a matrix
traceM=sp.Trace(A)

# Generate C code
print(sp.ccode(traceM))

If I print (sp.pprint) matrix A, I will get:

In [49]:sp.pprint(sp.Matrix(A))

If I print the the trace of A, the following error appears

In [50]:sp.pprint(sp.Matrix(traceM))

TypeError: 
Data type not understood; expecting list of lists or lists of values.

I was hopping to get .

Additionally, If I try to generate C code from the trace I will get the following message

In [51]: print(sp.ccode(traceM))
// Not supported in C:
// Trace
Trace(A)

and I was hopping for:

A[0, 0]+A[1, 1]+A[2, 2]

Can anyone give me a hand with this?

Note: if I use a numpy function (traceM=numpy.trace(A)) it will give me the expected result... but I should be able to obtain the same with sympy...

Best Regards,

sympy matrix
1 Answers

So, I think the aim here would be to unroll that trace expression, and have it be replaced by the explicit sum. The only way I found to do that unrolling process was through the use of rewrite (I was hinted at this because the Trace class has a method called _eval_rewrite_as_Sum )

The module being used to generate the C source code is the codegen module (also see Aaron Meurer's nice codegen tutorial and its github repo).

This was tested on SymPy 1.7

import sympy as sp
from sympy.utilities.codegen import codegen
N=3
A = sp.MatrixSymbol('A', N, N)
traceM = sp.Trace(A).rewrite(sp.Sum)
[(c_name, c_code), (h_name, c_header)] = codegen(("f", traceM), "C89", "test", header=False, empty=False)
print(c_code)

The result was this:

#include "test.h"
#include <math.h>
double f(double *A) {
   double f_result;
   f_result = A[0] + A[4] + A[8];
   return f_result;
}

One thing to notice is that the 2D array A is accessed as a 1-dimensional array.

user_8275190
Related