[cig-commits] [commit] devel: updates routine lagrange_any() for slightly better performance (b47b76b)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/specfem3d_globe

On branch  : devel
Link       : https://github.com/geodynamics/specfem3d_globe/compare/0e5b55c6f30be94583639fd325373eecd6facc6d...8be3e0b0267c8d4cf5af3bc26e8903da17bc4fd1

>---------------------------------------------------------------

commit b47b76b5905b2229f050eae12978f6ccc05015a0
Author: daniel peter <peterda at ethz.ch>
Date:   Mon Nov 24 15:38:33 2014 +0100

    updates routine lagrange_any() for slightly better performance


>---------------------------------------------------------------

b47b76b5905b2229f050eae12978f6ccc05015a0
 src/shared/lagrange_poly.f90 | 35 ++++++++++++++++++++++++++---------
 1 file changed, 26 insertions(+), 9 deletions(-)

diff --git a/src/shared/lagrange_poly.f90 b/src/shared/lagrange_poly.f90
index b762f23..49a684b 100644
--- a/src/shared/lagrange_poly.f90
+++ b/src/shared/lagrange_poly.f90
@@ -38,32 +38,49 @@
   double precision,dimension(NGLL),intent(in) :: xigll
   double precision,dimension(NGLL),intent(out) :: h,hprime
 
+  ! local parameters
   integer :: dgr,i,j
-  double precision :: prod1,prod2
+  double precision :: prod1,prod2,prod3
+  double precision :: prod2_inv
+  double precision :: sum
+  double precision :: x0,x
+
+! note: this routine is hit pretty hard by the mesher, optimizing the loops here will be beneficial
 
   do dgr = 1,NGLL
 
     prod1 = 1.0d0
     prod2 = 1.0d0
+
+    ! lagrangian interpolants
+    x0 = xigll(dgr)
     do i = 1,NGLL
       if (i /= dgr) then
-        prod1 = prod1*(xi-xigll(i))
-        prod2 = prod2*(xigll(dgr)-xigll(i))
+        x = xigll(i)
+        prod1 = prod1*(xi-x)
+        prod2 = prod2*(x0-x)
       endif
     enddo
-    h(dgr)=prod1/prod2
 
-    hprime(dgr) = 0.0d0
+    ! takes inverse to avoid additional divisions 
+    ! (multiplications are cheaper than divisions)
+    prod2_inv = 1.d0/prod2
+
+    h(dgr) = prod1 * prod2_inv
+
+    ! first derivatives
+    sum = 0.0d0
     do i = 1,NGLL
       if (i /= dgr) then
-        prod1 = 1.0d0
+        prod3 = 1.0d0
         do j = 1,NGLL
-          if (j /= dgr .and. j /= i) prod1 = prod1*(xi-xigll(j))
+          if (j /= dgr .and. j /= i) prod3 = prod3*(xi-xigll(j))
         enddo
-        hprime(dgr) = hprime(dgr)+prod1
+        sum = sum + prod3
       endif
     enddo
-    hprime(dgr) = hprime(dgr)/prod2
+
+    hprime(dgr) = sum * prod2_inv
 
   enddo