Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
2  
I don't have time to write a full answer, but a couple of things that might help are HornerForm for polynomials and the OptimizeExpression code used by Compile. For the latter, see 1 and 2 and linked questions. –  Simon Jan 14 '13 at 7:19
    
Hi Simon, thank you very much for your comment. Could you please give me more detailed information about OptimizeExpression? I checked the help and only find the C code example. Thanks, Tang Laoya –  Tang Laoya Jan 14 '13 at 7:32
1  
@user58211 I've written up a quick example. Also, can you edit your profile so that it displays your name/alias? People are often more willing to help you if you're more than just a number! –  Simon Jan 14 '13 at 9:44
    
Did you really write FormatType -> FotranForm in your code? That's not gonna work. –  librik Jan 14 '13 at 12:17
add comment

2 Answers 2

up vote 22 down vote accepted

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.

share|improve this answer
    
Dear all, Since x have 6 roots, how to let it output all these roots to Fortran code, such as x(1)=... x(2)=... ... Thanks –  Tang Laoya Feb 8 at 7:16
    
Furthermore, I'd like to output the Fortran code to a file, by the For loop, how to do? –  Tang Laoya Feb 8 at 10:54
add comment

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.

share|improve this answer
1  
……Because it's not free? –  xzczd Feb 9 at 5:12
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.