3
$\begingroup$

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!

$\endgroup$
5
  • $\begingroup$ Also, you might want to learn about CForm and the undocumented Experiment`OptimizeExpression. (Search this site for it.) $\endgroup$ Commented Dec 22, 2020 at 11:55
  • $\begingroup$ 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. $\endgroup$ Commented Dec 22, 2020 at 11:57
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Dec 22, 2020 at 11:58
  • $\begingroup$ What exactly do you want to use as input and return types? Or do you want to return the result as reference? $\endgroup$ Commented Dec 22, 2020 at 12:09
  • $\begingroup$ 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 $\endgroup$
    – user75941
    Commented Dec 22, 2020 at 12:12

1 Answer 1

4
$\begingroup$

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
  ]
 ]
$\endgroup$
1
  • $\begingroup$ I am not sure what you mean and why you want to do that. $\endgroup$ Commented Dec 22, 2020 at 21:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.