New tests and summary

.gitignore

@@ -18,3 +18,4 @@ tests/oxlog/bosou
 tests/oxlog/bosonf90
 tests/oxlog/eelsou
 tests/phintlog/eelsf90
+tests/phint optimization/phintdriver

tests/phint optimization/README & RESULTS

+These programs compare different versions of phint.
+
+Phint is calculated for 
+  0 <= phi <= π
+  0.2 <= u <= 2.4
+  a = 1.73205
+
+Phint: original from Lambin
+Phint1: new formula 
+Phint2: new formula, use conjugate complex
+Phint3: new formula, partial formula ?
+Phint4: new formula, partial formula ?
+Phint5: new formula, 
+Phint6: new formula, consider conjugate complex property
+Phint7: new formula, consider conjugate complex property, extract common sub-expressions 
+Phint8: new formula, consider conjugate complex property, extract common sub-expressions, treat phi = π
+Phint9: new formula, consider conjugate complex property, extract common sub-expressions, treat phi = π, use rational approximation of original Lambin code
+Phint10: new formula, consider conjugate complex property, no common sub-expression, treat phi = π, leave everything to the compiler.
+
+Log:
+
+>% gfortran -O0 phintdriver.f90 phint.f90 phint1.f90 phint2.f90 phint3.f90 phint4.f90 phint5.f90 phint6.f90 phint7.f90 phint8.f90 phint9.f90 phint10.f90 -o phintdriver
+>% phintdriver  
+Time for phint:  9.277 s
+Time for phint1: 17.785 s
+Time for phint2: 14.707 s
+Time for phint3:  8.899 s
+Time for phint4: 17.914 s
+Time for phint5: 25.697 s
+Time for phint6: 12.598 s
+Time for phint7:  6.052 s
+Time for phint8:  6.455 s
+Time for phint9:  6.594 s
+Time for phint10:  7.036 s
+
+>% gfortran -O7 phintdriver.f90 phint.f90 phint1.f90 phint2.f90 phint3.f90 phint4.f90 phint5.f90 phint6.f90 phint7.f90 phint8.f90 phint9.f90 phint10.f90 -o phintdriver
+>% phintdriver
+Time for phint:  8.648 s
+Time for phint1: 13.632 s
+Time for phint2: 13.302 s
+Time for phint3:  8.099 s
+Time for phint4: 11.682 s
+Time for phint5: 14.071 s
+Time for phint6:  8.823 s
+Time for phint7:  5.859 s
+Time for phint8:  5.008 s
+Time for phint9:  5.023 s
+Time for phint10:  4.963 s
+
+Summary:
+Considering the conjugate complex property is a major plus. As expected.
+Manually extracting common sub-expression is a plus with -O0, but no advantage, maybe even a small disadvantage with -O7. Not really expected.
+The rational approximation of the original Lambin code does not give any advantage.

tests/phint optimization/dooptimization.sh

-#!/bin/sh -v
-for number in 001 002 003 004 005 006 007 008 009 010 011 012 013 014 015 016 017 018 019 020 021 022 023 024 025 ; do
-    echo '*** File number' $number '***'
-
-    cp ../inputFiles/'eelsin'$number eelsin
-    ./eelsf90 >/dev/null
-    cp phint.log logFiles/'phint'$number'.log'
-
-    rm -f eelsin eelsou
-done

tests/phint optimization/phint10.f90

+double precision function phint10(phi, a, u)
+
+! ******************************************************************
+! *                                                                *
+! * evaluate the integral from zero to phi of                      *
+! *                                                                *
+! *                 u                 2                            *
+! *  ( ----------------------------- )  dphi                       *
+! *                          2     2                               *
+! *    (1 - a * u * cos(phi))  +  u                                *
+! *                                                                *
+! * for 0 <= phi <= pi , u >= 0 and a >= 0                         *
+! *                                                                *
+! ******************************************************************
+
+  implicit none
+
+  double precision, intent(in) :: phi
+  double precision, intent(in) :: a
+  double precision, intent(in) :: u
+
+  double precision :: e, pi, tanPhi2, sinPhi, cosPhi, Term1
+
+  double complex :: preFactor
+
+  if (phi == 0) then
+  	phint10 = 0
+  else
+    e = a*u
+    pi = 4 * atan (1.0d0)
+	preFactor = cmplx((e**2 - 1)/u + 2*u, 3) / sqrt(cmplx(e**2 - 1 + u**2, 2*u))**3
+	if (pi - phi < 1e-6) then
+	  phint10 = pi / 2 * real(preFactor)
+	else
+	  tanPhi2 = tan(phi/2)
+	  sinPhi = sin(phi)
+	  cosPhi = cos(phi)
+      Term1 = real(preFactor * atan(cmplx(u, 1 + e) / sqrt(cmplx(e**2 - 1 + u**2, 2*u)) * tanPhi2))
+      
+      phint10 = 																					&
+        Term1 +																					&
+          0.5 * e * (e**2 - 1 + 3*u**2 - e*(e**2 - 1 + u**2)*cosPhi) * sinPhi /					&
+            (((e**2 - 1)**2 + 2*(1 + e**2)*u**2 + u**4) * (1 + u**2 + e*cosPhi*(-2 + e*cosPhi)))
+     endif
+
+  endif
+
+  return
+end function phint10

tests/phint optimization/phint8.f90

@@ -40,11 +40,11 @@ double precision function phint8(phi, a, u)
 	  tanPhi2 = tan(phi/2)
 	  sinPhi = sin(phi)
 	  cosPhi = cos(phi)
-	  Term1 = real(preFactor	* atan(cmplx(u,1+e) / denominator * tanPhi2))
+	  Term1 = real(preFactor * atan(cmplx(u,1+e) / denominator * tanPhi2))
 	
-	  phint8 = 																	&
-	    Term1 +																	&
-		  0.5 * e * (e21 + 3*u2 - e*(e21 + u2)*cosPhi) * sinPhi /				&
+	  phint8 = 																		&
+	    Term1 +																		&
+		  0.5 * e * (e21 + 3*u2 - e*(e21 + u2)*cosPhi) * sinPhi /					&
 		    ((e21**2 + 2*(1 + e2)*u2 + u2**2) * (1 + u2 + e*cosPhi*(-2 + e*cosPhi)))
     endif
 

tests/phint optimization/phint9.f90

+double precision function phint9(phi, a, u)
+
+! ******************************************************************
+! *                                                                *
+! * evaluate the integral from zero to phi of                      *
+! *                                                                *
+! *                 u                 2                            *
+! *  ( ----------------------------- )  dphi                       *
+! *                          2     2                               *
+! *    (1 - a * u * cos(phi))  +  u                                *
+! *                                                                *
+! * for 0 <= phi <= pi , u >= 0 and a >= 0                         *
+! *                                                                *
+! ******************************************************************
+
+  implicit none
+
+  double precision, intent(in) :: phi
+  double precision, intent(in) :: a
+  double precision, intent(in) :: u
+
+  double precision :: e, e2, e21, pi, u2, zr, esr, tanPhi2, sinPhi, cosPhi, Term1
+
+  double complex :: nominator, denominator, preFactor
+
+  if (phi == 0) then
+  	phint9 = 0
+  else
+    e = a*u
+	e2 = e*e
+    u2 = u*u
+	e21 = e2 - 1
+    pi = 4 * atan (1.0d0)
+	nominator   =      cmplx(e21/u + 2*u ,3)
+	denominator = sqrt(cmplx(e21   +   u2,2*u))
+	preFactor   = nominator / denominator**3
+	if (pi - phi < 1e-6) then
+	  phint9 = pi / 2 * real(preFactor)
+	else
+	  tanPhi2 = tan(phi/2)
+	  sinPhi = sin(phi)
+	  cosPhi = cos(phi)
+      if (u < 1.0d0 .and. e < 1.0d-02 * (1.0d0 + u2)) then
+        zr = 1.0d0 + u2
+        esr = e / zr
+        phint9 = u2 / zr**2 * ((( (4.0d0 / 3.0d0) * (2.0d0 + cosPhi**2) * sinPhi * (5.0d0 - 3 * u2) *  &
+                esr + (phi + cosPhi * sinPhi) * (5.0d0 - u2)) * esr + 4 * sinPhi) * esr + phi)
+      else
+        Term1 = real(preFactor * atan(cmplx(u,1+e) / denominator * tanPhi2))
+        
+        phint9 = 																		&
+          Term1 +																		&
+            0.5 * e * (e21 + 3*u2 - e*(e21 + u2)*cosPhi) * sinPhi /						&
+              ((e21**2 + 2*(1 + e2)*u2 + u2**2) * (1 + u2 + e*cosPhi*(-2 + e*cosPhi)))
+       endif
+     endif
+
+  endif
+
+  return
+end function phint9

tests/phint optimization/phintdriver.f90

@@ -3,7 +3,8 @@ program phintdriver
   implicit none
 
   double precision :: phi, a, u, result
-  double precision :: phint, phint1, phint2, phint3, phint4, phint5, phint6, phint7, phint8
+  double precision :: phint, phint1, phint2, phint3, phint4, phint5
+  double precision :: phint6, phint7, phint8, phint9, phint10
   double precision :: step_phi, step_u
   integer :: i, j, number
   
@@ -114,4 +115,26 @@ program phintdriver
   call cpu_time(finish)
   print '("Time for phint8: ", f6.3, " s")',finish-start
   
+  call cpu_time(start)
+  do i = 0, number
+    phi = i * step_phi
+    do j = 0, number
+      u = 0.2 + j * step_u
+      result = phint9(phi, a, u)
+    enddo
+  enddo
+  call cpu_time(finish)
+  print '("Time for phint9: ", f6.3, " s")',finish-start
+  
+  call cpu_time(start)
+  do i = 0, number
+    phi = i * step_phi
+    do j = 0, number
+      u = 0.2 + j * step_u
+      result = phint10(phi, a, u)
+    enddo
+  enddo
+  call cpu_time(finish)
+  print '("Time for phint10: ", f6.3, " s")',finish-start
+  
   end program phintdriver