# Generate C code from symbolic computation in Mathematica

I am doing some tensor multiplication operations and would like to have Mathematica generate C code so that I can use it later.

I want to supply two 2nd order tensors(A and B) and get the entries of a 4th order tensor as a function of the inputs A and B.

However, I am having problems getting the output.

Currently, I have:

ClearAll[a,b]

(* define tensor a and b as symmetric*)
SetAttributes[a,Orderless];
SetAttributes[b,Orderless];

tensorA = Array[a,{3,3},0];

tensorB = Array[b,{3,3},0];

g= 1-4*Det[tensorA];

(* define my operation *)
result = TensorProduct[tensorA,tensorB]+g*TensorTranspose[TensorProduct[tensorA,tensorB],{1,3,2,4}]

(* Generate C code *)

(*How can I generate C code?*)

CAssign[result]


Additionally, to generate better C code, should I define:

g= 1-4*Det[tensorA];


or

g= 1.0-4.0*Det[tensorA];


With the first I get the result of CAssign as fractions, with the second as floating point. Does this make a difference?

The expected output should be something like:

// entries of the input tensors
double var1 = a(0,0)
double var5 = a(1,1)
double var9 = a(2,2)
double var10 = b(0,0)
double var14 = b(1,1)
double var18 = b(2,2)

result(0,0,0,0)=var1 * var10
result(0,0,0,1)=...
result(0,0,0,2)=...


and so on

Best Regards!

• Also, you might want to learn about CForm and the undocumented ExperimentOptimizeExpression. (Search this site for it.) Commented Dec 22, 2020 at 11:55
• Please, be a bit more specific about what sort of C code you expect. You want to generate just an expression, a function or a callable library? It is impossible to answer your question without knowing what exactly you want to do with the generated piece of code. Commented Dec 22, 2020 at 11:57
• Yes, use floating point numbers wherever possible; type casts do not come for free. One notable exception is the second argument of Power: When possible, use an integer there. For example, Compile is then able to unroll this to a couple of multiplications which are by an order of magnitude faster than pow in C. Commented Dec 22, 2020 at 11:58
• What exactly do you want to use as input and return types? Or do you want to return the result as reference? Commented Dec 22, 2020 at 12:09
• I just want the entries of tensor (to be written as a function of the two entry tensors) The return type will be a 4th order tensor structure whos entries will be given by the operations with the two input tensors
– user75941
Commented Dec 22, 2020 at 12:12

Maybe this does approximately what you want. It is not a general solution but specifically written for this question. No quarantees, neither for correctness nor compilability.

Needs["SymbolicC"];

Block[{a, b, var, r},

MyN[code_] :=
N[code] /. {Times[-1., s__] :> -s, Times[1., s__] :> s};
ClearAll[part];
SetAttributes[part, NHoldRest];
With[{format = Format},
part /: format[part[a_, {i_, j_, k_, l_}], CForm] :=
With[{ii = (i - 1), jj = (j - 1), kk = (k - 1), ll = (l - 1),
aa = a},
HoldForm[aa[[ii, jj, kk, ll]]]
];
part /: format[part[a_, {i_, j_}], CForm] :=
With[{ii = (i - 1), jj = (j - 1), aa = a},
HoldForm[aa[[ii, jj]]]
];
];
tensorA = Table[part[a, {i, j}], {i, 1, 3}, {j, 1, 3}];
tensorA =
UpperTriangularize[tensorA] +
Transpose[UpperTriangularize[tensorA, 1]];
tensorB = Table[part[b, {i, j}], {i, 1, 3}, {j, 1, 3}];
tensorB =
UpperTriangularize[tensorB] +
Transpose[UpperTriangularize[tensorB, 1]];
g = 1 - 4*Det[tensorA];

(*define my operation*)

result = TensorProduct[tensorA, tensorB] +
g*TensorTranspose[TensorProduct[tensorA, tensorB], {1, 3, 2, 4}];

(*optimize expression*)

expr0 = ExperimentalOptimizeExpression[MyN@result,
"OptimizationLevel" -> 2, "OptimizationSymbol" -> var] /.
CompoundExpression -> List /. Power[a_, 2] :> HoldForm[a a];

(*declare local variables, convert to CForm expressions, and turn \
Set into CAssign*)
expr1 = expr0 /. HoldPattern[Set[a_, b_]] :>
With[{
lhs = ToString[a, CForm],
rhs = ToString[b, CForm]
},
CAssign["double " <> lhs, rhs]
];

(*convert return value into a set of CAssign operations*)

expr1[[1, 2, -1]] = MapIndexed[
CAssign[ToString[part[r, #2], CForm], ToString[#1, CForm]] &,
expr1[[1, 2, -1]],
{4}
];

(*generate the C code*)
GenerateCode[
CBlock[expr1[[1, 2]]],
Indent -> 1
]
]
`
• I am not sure what you mean and why you want to do that. Commented Dec 22, 2020 at 21:35