Generate C code from symbolic computation in Mathematica

Question

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!

Answer 1

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 = Experimental`OptimizeExpression[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
  ]
 ]