[VB.NET 4.0] Can you speed up this function? :)

[DragonQ]

I have some code for calculating a "polytrope" (a good model for some types of star). I was reading some optimisation threads on this forum today and it made me curious enough to try to make it faster. It was already probably the most efficient program I'd made since it was my first truly multithreaded application (that actually benefited hugely from it).

Moving "Math.PI" into a constant so it's only called once didn't make any measurable difference in the time taken. However, replacing two occurrences of "r ^ 2" with a new variable initialised to "r ^ 2", which was then used twice, improved the speed of the code by 35%!

I don't actually need this code as it's from an old program but I'm just curious to see if there are any other efficiency improvements one could make?

P.S. This is running in 64-bit mode so I don't think changing any of the Doubles into Singles will improve speed.


Public Const G As Double = 0.00000000006674 'gravitational constant
Public Const Pi As Double = Math.PI 'Pi

Public Function CalculatePolytrope(ByVal n As Double, ByVal StepSize As Double, ByVal Rho As Double, ByVal P As Double) As Double()

Dim K As Double = P / (Rho ^ ((n + 1) / n)) 'constant of polytrope model
Dim m As Double = 0 'mass
Dim r As Double = 0 - StepSize 'radius
Dim dP_dr As Double 'dP/dr
Dim dm_dr As Double 'dm/dr
Dim rSquared As Double 'r^2

'Process loop until Rho or P are equal to (or less than) zero:
Do
'Increment r
r += StepSize

'This prevents dP/dr and dm/dr becoming NaN ("not a number") when r = 0
If r = 0 Then
dP_dr = 0
dm_dr = 0
Else
rSquared = r ^ 2
dP_dr = -(G * m * Rho) / rSquared
dm_dr = 4 * Pi * Rho * rSquared
End If

'Calculate new mass and pressure
m += StepSize * dm_dr
P += StepSize * dP_dr

'This prevents square rooting a minus number, which would make Rho NaN ("not a number")
If P > 0 Then
Rho = (P / K) ^ (n / (n + 1))
Else
Rho = 0
End If
Loop Until Rho <= 0 OrElse P <= 0

'Create an array of doubles to hold our final M and R then return them to the calling subroutine.
Dim Results() As Double = {m, r}
Return Results

End Function

[si_the_geek]

A minor improvement would be to alter lines 19 to 30, because m and p will only be altered if r is not 0. As such, you could use this to get the same behaviour:
If r <> 0 Then
rSquared = r ^ 2
dP_dr = -(G * m * Rho) / rSquared
dm_dr = 4 * Pi * Rho * rSquared
'Calculate new mass and pressure
m += StepSize * dm_dr
P += StepSize * dP_dr
End If...and maybe even eliminate dP_dr and dm_dr.


While storing Pi didn't do you much good, storing 4*Pi (which is what you actually use) is likely to be better.


As n doesn't change inside the loop, it would be a good idea to store (n / (n + 1)) to a variable too.


Due to the conditions/values used in lines 32 to 38, you could simplify that a bit too, eg:
'This prevents square rooting a minus number, which would make Rho NaN ("not a number")
If P <= 0 Then Exit Do
Rho = (P / K) ^ (n / (n + 1))
Loop Until Rho <= 0

[DragonQ]

Thanks for the suggestions! I've read that using >= or <= is slower than > or < so I've tried to avoid using them in this new version. By setting Rho to -1 instead of 0 if "P > 0" is false, I can change the loop condition to "Rho < 0", which should be a bit quicker.

Below is what I have now. This version is 4% faster (averaged over 5 tests) than the version I posted above (which was in turn 35% faster than the code I started with):

'
Public Const G As Double = 0.00000000006674 'gravitational constant
Public Const Pi_x4 As Double = Math.PI * 4 '4 Pi

Public Function CalculatePolytrope(ByVal n As Double, ByVal StepSize As Double, ByVal Rho As Double, ByVal P As Double) As Double()

Dim PolytropeIndex As Double = n / (n + 1) 'n / (n + 1)
Dim K As Double = P / (Rho ^ (1 / PolytropeIndex)) 'constant of polytrope model
Dim m As Double = 0 'mass
Dim r As Double = 0 - StepSize 'radius
Dim dP_dr As Double 'dP/dr
Dim dm_dr As Double 'dm/dr
Dim rSquared As Double 'r^2

'Process loop until Rho or P are equal to (or less than) zero:
Do
'Increment r
r += StepSize

'This prevents dP/dr and dm/dr becoming NaN ("not a number") when r = 0
If r = 0 Then
dP_dr = 0
dm_dr = 0
Else
rSquared = r ^ 2
dP_dr = -(G * m * Rho) / rSquared
dm_dr = Pi_x4 * Rho * rSquared
End If

'Calculate new mass and pressure
m += StepSize * dm_dr
P += StepSize * dP_dr

'This prevents square rooting a minus number, which would make Rho NaN ("not a number")
If P > 0 Then
Rho = (P / K) ^ PolytropeIndex
Else
Rho = -1
End If
Loop Until Rho < 0

'Create an array of doubles to hold our final M and R then return them to the calling subroutine.
Dim Results() As Double = {m, r}
Return Results

End Function

Here's a couple of other things I've tried without measurable effect:

Making the "G" variable negative so it doesn't need to be made negative each time dP_dr is calculated.
Moving the negative sign from the dP_dr calculation to the "P += ..." line (change to "P -= ...")

And here's some more things I've tried that all make the code slower (bizzarly):

Changing the "If r = 0" condition to "If r <> 0" or "If r > 0" or "If Not r = 0" and swapping the If/Else sections around.
Changing "Rho = -1" to "Exit Do".
Putting the "m += ..." and "P += ..." lines inside the Else section of the "If r = 0" condition.
Getting rid of dP_dr and dm_dr in conjunction with the previous change.
Commenting out the unneeded "dP_dr = 0" and "dm_dr = 0" lines.

Strange...

[Zach_VB6]

Instead of initializing a variable to return, you can return directly: Return New Double() {m, r}

Instead of rSquared = r ^ 2 tryr=r*r

(ya I know it's minor)

[DragonQ]

Woah...replacing "r ^ 2" with "r * r" improved speed by a further 69%! Ridiculous. That means compared to what I started with, the code is now ~140% faster (~2.4x). :D

EDIT: After some further trial-and-error, the fastest I can get so far is a ~165% improvement (2.65x). Here's the code for this:

'
Public Const G As Double = 0.00000000006674 'gravitational constant
Public Const Pi_x4 As Double = Math.PI * 4 '4π

Public Function CalculatePolytrope(ByVal n As Double, ByVal StepSize As Double, ByVal Rho As Double, ByVal P As Double) As Double()

Dim PolytropeIndex As Double = n / (n + 1) 'n / (n + 1)
Dim K As Double = P / (Rho ^ (1 / PolytropeIndex)) 'constant of polytrope model
Dim m As Double = 0 'mass
Dim r As Double = 0 - StepSize 'radius
Dim dP_dr As Double 'dP/dr (excluding Rho term)
Dim rSquared As Double 'r^2
Dim StepSize_x_Rho As Double 'StepSize x Rho

'Process loop until Rho or P are equal to (or less than) zero:
Do
'Increment r
r += StepSize

'This prevents dP/dr and dm/dr becoming NaN ("not a number") when r = 0
If r > 0 Then
rSquared = r * r
dP_dr = (G * m) / rSquared
StepSize_x_Rho = StepSize * Rho

'Calculate new mass and pressure
m += (StepSize_x_Rho * Pi_x4 * rSquared)
P -= (StepSize_x_Rho * dP_dr)
End If

'This prevents square rooting a minus number, which would make Rho NaN ("not a number")
If P > 0 Then
Rho = (P / K) ^ PolytropeIndex
Else
Rho = -1
End If
Loop Until Rho < 0

'Create an array of doubles to hold our final M and R then return them to the calling subroutine.
Return {m, r}

End Function

[si_the_geek]

Is it possible for the the Rho value to ever be 0 or less from the calculation? If so are you sure that your loop exit condition is apt?

In terms of speed I suspect the quickest code for it would be Loop Until Rho = 0 , but that could lead to a change in behaviour.

[DragonQ]

That is true, the condition could be changed to "Rho = 0" without a change in behaviour. It doesn't improve speed though, unfortunately.

The original program managed to perform the calculation with a range of input values that I originally made the code for in 414 s (5000 runs, 8 threads at a time, i7 920 @ 3.8 GHz). The current version (post 5) does it in 154 s. :)

[minitech]

I don't know if this does anything for speed, but stop initializing variables to 0, and replace 0 - stepSize with just -stepSize, it looks better :)

At the end instead of setting Rho to -1 you could use GoTo and maybe take out the loop, I wonder how that would work?

For conciseness, take out all those unnecessary parentheses (like around Stepsize_x_Rho * Pi_x4 * rSquared and G * m and Stepsize_x_Rho * dp_dr).

[Nightwalker83]

Actually, r ^ 2 = r * r. Sorry.

Yeah, I tried explaining it but I wrote the equation incorrectly thinking r*r*r = r ^2 when it was meant to be r *r = r ^2. I deleted my original post because it was too long-winded and confusing. Although, the above still doesn't seem to explain why r*r would be faster than r^2?

Edit:

Maybe the fact that DragonQ had rSquared dimmed as double had something to do with it.

[minitech]

^ is an operator that implicitly uses Doubles on both sides and performs a multitude of completely unnecessary operations to accomodate decimal and negative powers. Multiplication, being a fundamental binary operation, is much faster.

[dbasnett]

@Dragon - post sample data in the following format please

Dim smplData() As Double = New Double() {n, StepSize, Rho, P, M_answer, R_answer, _
n, StepSize, Rho, P, M_answer, R_answer, _
n, StepSize, Rho, P, M_answer, R_answer}


replacing the variables with the actual values. Two or three sets should be good, more if you want.

[DragonQ]

Dim smplData() As Double = New Double() {1, 1000, 4.6E+1, 4.4E+10, 2.005282E+028, 6.995690E+008, _
1, 1000, 4.6E+3, 4.4E+12, 2.005340E+027, 6.995200E+007, _
1, 1000, 4.6E+5, 4.4E+14, 2.005925E+026, 6.991000E+006}

[dbasnett]

Maybe a little better, but it is hard to tell...


Private Sub Button2_Click(sender As System.Object, _
e As System.EventArgs) Handles Button2.Click

Dim sd() As Double = New Double() {1, 1000, 46.0, 44000000000.0, 2.0052816957994727E+28, 699569000.0, _
1, 1000, 4600.0, 4400000000000.0, 2.0053401856028196E+27, 69952000.0, _
1, 1000, 460000.0, 440000000000000.0, 2.0059251372079435E+26, 6991000.0}
Dim stpw As New Stopwatch
stpw.Start()
For y As Integer = 1 To 20
For x As Integer = 0 To sd.Length - 1 Step 6
Dim ans() As Double = CalculatePolytrope(sd(x), sd(x + 1), sd(x + 2), sd(x + 3))
If ans(0) = sd(x + 4) AndAlso ans(1) = sd(x + 5) Then
'Debug.WriteLine("m={0} r={1}", ans(0), ans(1))
Else
Stop
End If
Next
Next y
stpw.Stop()
'Debug.WriteLine(stpw.ElapsedMilliseconds.ToString("n0"))
TextBox1.Text = stpw.ElapsedMilliseconds.ToString("n0")
End Sub

Public Const G As Double = 0.00000000006674# 'gravitational constant
Public Const Pi_x4 As Double = Math.PI * 4 '4π

Public Function CalculatePolytrope(ByVal n As Double, _
ByVal StepSize As Double, _
ByVal Rho As Double, _
ByVal P As Double) As Double()

Dim PolytropeIndex As Double = n / (n + 1) 'n / (n + 1)
Dim K As Double = P / (Rho ^ (1 / PolytropeIndex)) 'constant of polytrope model
Dim m As Double 'mass
Dim r As Double = -StepSize 'radius
Dim dP_dr As Double 'dP/dr (excluding Rho term)
Dim rSquared As Double 'r^2
Dim StepSize_x_Rho As Double 'StepSize x Rho

'Process loop until Rho or P are equal to (or less than) zero:
Do
'Increment r
r += StepSize

'This prevents dP/dr and dm/dr becoming NaN ("not a number") when r = 0
If r > 0 Then
rSquared = r * r
dP_dr = (G * m) / rSquared
StepSize_x_Rho = StepSize * Rho

'Calculate new mass and pressure
m += (StepSize_x_Rho * Pi_x4 * rSquared)
P -= (StepSize_x_Rho * dP_dr)
If P <= 0 Then Exit Do
End If

Rho = (P / K) ^ PolytropeIndex
Loop Until Rho < 0

'Create an array of doubles to hold our final M and R then return them to the calling subroutine.
Return {m, r}

End Function

[DragonQ]

Nice try but it's no quicker. :(

It's probably a good idea to give you the settings used for my benchmark which gives me 4.10-4.15 seconds with my current code (the slowest of the 9 original calculations the program was intended for):

'
Dim StartTime as DateTime = DateTime.Now

Dim Results As Double() = CalculatePolytrope(1, 1000, 46.0, 440000000000000.0)

Dim TimeElapsed As TimeSpan = DateTime.Now.Subtract(StartTime)
txtTimeElapsed.Text = Format((TimeElapsed.TotalMilliseconds / 1000), "##0.00").ToString & " s"

[leinad31]

(a/b)^x same as a^x / b^x

K ^ PolytropeIndex is a constant

[DragonQ]

Good suggestion but again that idea ends up slowing the program down (4.25-4.30 s vs 4.10-4.15 s).

[leinad31]

I see.

I'm curious though why you didn't use a formula based on dM/dr since you are after mass anyway. isn't it possible? Maybe the algorithm and computation would be more efficient with that approach.