# Replacement Patterns in FortranForm

I have recently asked a question concerning string replacement rules for FortranForm. This question has been answered by Carl Woll - thanks for that!

Now as I was implementing this I was thinking of overcoming number definitions and simply replace real numbers with something like 1. --> 1._rkind. This replacement should work for all fractions, powers, etc.

I have found a post where this task has been tackled but unfortunately my Mathematica skills are very basic and I do not really understand what is going on there. Anyway, I have tried and this is what I was coming up so far:

G = -4/3 pi (-RAi + RAip1) (RAi^2 + RAi RAip1 + RAip1^2) rhoAi;
testReplace[someEq_] := Module[{trep},
Clear[fort, fortranreal];
fort[x_Real] := fortranreal[x];
fort[x_] := x;
trep = StringReplace[ToString@MapAll[fort, FortranForm[#]], Shortest["fortranreal(" ~~ s : __ ~~ ")"] :> s <> "_rkind"] & /@ someEq
];


This leaves me with the (unfortunately not quite useful) output:

(-4*pi*(-RAi + RAip1)*(RAi**2 + RAi*RAip1 + RAip1**2)*rhoAi)/3.


because obviously 4 is not an integer and would cause problems with the Fortran compiler. Next thing I've tried was replacing someEq with N[someEq] in the Module output which then leaves me with:

-1.3333333333333333_rkind pi RAi**2 + RAi*RAip1 + RAip1**2 RAip1 + RAi*-1._rkind rhoAi


where all the brackets are stripped away and e.g. -RAi gets replaced with RAi*-1._rkind.

What I am trying to achieve is evaluating all fractions and then attach _rkind to all real values. Can anyone tell me what am I doing wrong?

• FWIW the first expression is fine. Your Fortran compiler will automatically upcast the 4 to the type of pi. You do have to watch if you end up with actual integer division (integer type both numerator and denominator) in your expressions. – george2079 Mar 8 '18 at 12:55
• Hey! Thanks, but unfortunately this doesn't apply to every Fortran compiler. Moreover the 3. is missing something like 3._rkind as well. – Florian Ragossnig Mar 8 '18 at 16:19
• Yes, fortran standard requires automatic upcasting. specifying a kind for 3. is also unnecessary because its an integer value. Only if you have actual decimals ( say 0.3 ) , then you need to specify the type if you want it to be other than single precision. – george2079 Mar 8 '18 at 17:07
• Thanks for that! I obviously was lacking knowledge here! I will run a few more tests and post an answer once that is finished. – Florian Ragossnig Mar 9 '18 at 9:18

After running a few tests and digging deeper into Fortran standards (thanks @george2079 for pushing this into the right direction), I can finally post an aswer.

The method above is right but will cause problems when adding negative reals. As soon as we plug in e.g. 4./3. it will get converted to a real (1.33) and because the constant literal of 1.33 != 1.33_rkind the method attaches the _rkind to the number. Unfortunately Mathematica's Times[] will shuffle the number to the very end of the output. For positive values this does not matter but in case of negative values this will cause problems as we end up with e.g.:

G = -4./3 x;
testReplace[G]


Out:

x*-1.3333333333333_rkind >>> read update!


In order to fix this I've simply stopped Times[] by doing this with ClearAttributes[Times, Orderless]. So for completion the whole thing now looks like this:

Unprotect[Times];
ClearAttributes[Times, Orderless];
Protect[Times];

testReplace[someEq_] := Module[{trep}, Clear[fort, fortranreal];
fort[x_Real] := fortranreal[x];
fort[x_] := x;
trep = StringReplace[ToString@MapAll[fort, FortranForm[#]],
Shortest["fortranreal(" ~~ s : __ ~~ ")"] :> s <> "_rkind"] & /@ {someEq}[[1]]
];

G = -4./3 pi (-RAi + RAip1) (RAi^2 + RAi RAip1 + RAip1^2) rhoAi;
testReplace[G]


Out:

-1.333333333333333_rkind*pi*(-RAi + RAip1)*(RAi**2 + RAi*RAip1 + RAip1**2)*rhoAi


## UPDATE

Actually you there is no need of fixing the reordering here as a*-b, a+-b and all other variations are allowed (thanks @george for the comment!) under Fortran (ifort and gfortran), it's just a bit better for readability.

• the x*-1.3.. is actually legal. Looks weird though, I cant blame you for wanting to fix it. Another approach here might be to wrap parenthesis around negative numbers. – george2079 Mar 9 '18 at 16:48
• I actually did that. By the way I've tried the x*-b and x+-b thing and it actually works. The compiler gives a warning though but that doesn't matter. I'm doing all this because I have some complicated equations and need their derivatives and the only way to get them (in multiple different geometries) is by using Mathematica. So readability is not really an issue here. – Florian Ragossnig Mar 9 '18 at 18:32