How can I get Mathematica to produce better Fortran code?

Question

I am new to Mathematica and I am trying to get to it to produce some expressions in Fortran code. However, it seems that Mathematica will output duplicated expressions, that is, expressions the need not be calculated more than once. Is there any good method to define a rule so that Mathematica will recognize duplicated expressions, assign them to a new variable, and then optimize the Fortran code?

P.S:

The following is some code I wrote to test Mathematica's Fortran capabilities. The solution of x, y and z contains many duplicated expressions, for example, (-36*a + 20*a**3 - 216*b)**2.

sol = Solve[{x^2 + y^2 + z^2 == 1, x + y + z == a , x*y*z == b}, {x, y, z}];
xx = x /. sol
yy = y /. sol
zz = z /. sol
Print["Writing Fortran Code . . . : / "];
SetDirectory["F:\\tang\\mathtest"];
strm = OpenWrite["test.f90", FormatType -> FotranForm,  PageWidth -> 70];
(* write subroutine of invisopar*)
WriteString[strm, "subroutine test(x,y,z,a,b)\n"];
WriteString[strm, "implicit none\n"];
WriteString[strm, "real*8::x,y,z,a,b\n"];
nroot = Length[xx];
For[ii = 1, ii <= nroot,
  WriteString[strm,
   "x = " <> ToString[FortranForm[xx[[ii]]]] <> "\n"]; ii++];

For[ii = 1, ii <= nroot,
  WriteString[strm,
   "y = " <> ToString[FortranForm[yy[[ii]]]] <> "\n"]; ii++];

For[ii = 1, ii <= nroot,
  WriteString[strm,
   "z = " <> ToString[FortranForm[zz[[ii]]]] <> "\n"]; ii++];

WriteString[strm, "end subroutine\n"];
Close[strm];
Print["Finished Writing Fortran Code . . . : / "];

When the test.f90 is outputted, the first solution of x is:

x = a/3. + (-6 + 2*a**2)/(3.2*0.6666666666666666*(-36*a + 20*a*3 + Sqrt(4(-6 + 2*a*2)*3 + (-36*a + 20*a**3 - 216*b)**2) - 216*b)**0.3333333333333333) - (-36*a + 20*a*3 + Sqrt(4(-6 + 2*a*2)*3 + (-36*a + 20*a**3 - 216*b)**2) - 216*b)**0.3333333333333333/(6.2*0.3333333333333333)

I would like that the code is as follows:

tmp0 = 216*b
tmp1 = -36*a + 20*a**3
tmp2 = (tmp1 - tmp0 )**2
tmp3 = -6 + 2*a**2
tmp4 = 4*tmp3**3
tmp5 = Sqrt(tmp4 + tmp2 )
tmp6 = tmp1 + tmp5 - tmp0
tmp7 = tmp6**0.3333333333333333
tmp8 = 2**0.6666666666666666
x = a/3. + tmp3 /(3.*tmp8 *tmp7 ) - tmp7/(6.*tmp8)

I noticed that maple has the capability for optimizing Fortran; for example, the following maple code can output optimized Fortran code:

with(codegen, fortran)
A := array(1 .. 2, 1 .. 2, symmetric); 
A[1, 1] := log(x); A[1, 2] := 1-log(x); A[2, 2] := 2-log(x);
print(A);
fortran(A, optimized, mode = double)

So I want to know whether Mathematica also can do similar optimization.

Answer 1

OK, here's a quick go using Experimental`OptimizeExpression, which is used internally by Compile. I'm not claiming that this is a polished solution to your problem, but hopefully it can be used as a place to start. (Also, forgive any Fortran mistakes, I haven't looked at Fortran for a long time...)

xx = x /. Solve[{x^2 + y^2 + z^2 == 1, x + y + z == a, x*y*z == b}, {x, y, z}] // First;

Quiet[Remove["tmp*"]]  (* otherwise the numbers after 'tmp' keep going up *)
Module[{out}, 
  out = Experimental`OptimizeExpression[xx, "OptimizationSymbol" -> tmp] /. 
    HoldPattern[Experimental`OptimizedExpression[Block[{tempVars__}, 
      CompoundExpression[defs : Set[_, _] ..., expr_]]]] :> Block[{
        Set = ToString[#1] <> " = " <> ToString[FortranForm[#2]] &}, List[defs, expr]];
{"subroutine test(x,a,b)", "implicit none", "real*8::x,a,b",
 Most@out, "x = " <> ToString@FortranForm@N@Last@out, 
 "end subroutine"}
] // Flatten // ExportString[#, "Table"] &

subroutine test(x,a,b)
implicit none
real*8::x,a,b
tmp4 = a**2
tmp5 = 2*tmp4
tmp6 = -6 + tmp5
tmp7 = -36*a
tmp8 = a**3
tmp9 = 20*tmp8
tmp12 = -216*b
tmp10 = tmp6**3
tmp11 = 4*tmp10
tmp13 = tmp12 + tmp7 + tmp9
tmp14 = tmp13**2
tmp15 = tmp11 + tmp14
tmp16 = Sqrt(tmp15)
tmp17 = tmp12 + tmp16 + tmp7 + tmp9
x = 0.3333333333333333*a - 0.13228342099734997*tmp17**0.3333333333333333 + \
    (0.20998684164914552*tmp6)/tmp17**0.3333333333333333
end subroutine

This can easily be modified to process and output more than just the first solution for x and can also be made to export to a file instead of a string. It's not as compact as hand code would make it, some of the tmp definitions, I would join together, but it's not bad. Also, note that the other solutions for x contain complex numbers and might have to be treated a little more carefully.

Answer 2

Why not use MathCode?

MathCode makes it possible to develop prototypes in the interactive Mathematica environment which can be automatically translated to fast production code in C++ or Fortran90 and, if necessary, linked to external applications.