[cig-commits] [commit] master: py: Simplify Lagrange derivatives a bit. (d8f09c2)

[cig_noreply at geodynamics.org]

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

On branch  : master
Link       : https://github.com/geodynamics/specfem1d/compare/2ccfba13481e26296eebf034518cd31d2c727906...33b512de01491f2c70de206939855d95b22febf4

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

commit d8f09c22a6835b505027a48073525a430791a506
Author: Elliott Sales de Andrade <esalesde at physics.utoronto.ca>
Date:   Mon Jan 12 16:54:25 2015 -0500

    py: Simplify Lagrange derivatives a bit.
    
    Extract some loop-invariants and use some product functions.


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

d8f09c22a6835b505027a48073525a430791a506
 Python_version/gll.py | 48 +++++++++++++++++++++++++++++++-----------------
 1 file changed, 31 insertions(+), 17 deletions(-)

diff --git a/Python_version/gll.py b/Python_version/gll.py
index 06d141e..7f755df 100644
--- a/Python_version/gll.py
+++ b/Python_version/gll.py
@@ -56,25 +56,39 @@ def lagrange_derivative(ksiGLL):
     """Calculates the values of the derivative of the Lagrange polynomials
     at the GLL points"""
     N = len(ksiGLL) - 1
-    deriv=np.zeros((N+1,N+1),dtype='d')
-    for i in range(N+1):
-        for j in range(N+1):
+    deriv = np.zeros((N + 1, N + 1))
+
+    # AKA: diffs[i,j] = ksiGLL[i] - ksiGLL[j]
+    diffs = ksiGLL[:, np.newaxis] - ksiGLL[np.newaxis, :]
+
+    for i in range(N + 1):
+        # Exclude i since ksiGLL[i] - ksiGLL[i] is obviously 0.
+        prod1 = np.product(diffs[i, :i]) * np.product(diffs[i, i + 1:])
+        prod1 = 1 / prod1
+
+        for j in range(N + 1):
             if i == 0 and j == 0:
-                deriv[i,j] = -N*(N+1.)/4.
+                deriv[i, j] = -N * (N + 1) / 4.0
+
             elif i == N and j == N:
-                deriv[i,j] = N*(N+1.)/4.
-            elif i == j:
-                deriv[i,j] = 0.
-            else:
-                prod1=1.
-                for m in range(N+1):
-                    if m!=i:
-                        prod1 *= 1/(ksiGLL[i]-ksiGLL[m])
-                prod2=1.
-                for k in range(N+1):
-                    if k!=j and k!=i:
-                        prod2 *= (ksiGLL[j]-ksiGLL[k])
-                deriv[i,j]=prod1*prod2
+                deriv[i, j] = N * (N + 1) / 4.0
+
+            elif i < j:
+                # Just like prod1, but excluding i and j.
+                prod2 = (np.product(diffs[j, :i]) *
+                         np.product(diffs[j, i + 1:j]) *
+                         np.product(diffs[j, j + 1:]))
+
+                deriv[i, j] = prod1 * prod2
+
+            elif i > j:
+                # Just like prod1, but excluding j and i.
+                prod2 = (np.product(diffs[j, :j]) *
+                         np.product(diffs[j, j + 1:i]) *
+                         np.product(diffs[j, i + 1:]))
+
+                deriv[i, j] = prod1 * prod2
+
     return deriv