[cig-commits] r4466 - cs/benchmark/trunk/benchmark

[baagaard at geodynamics.org]

Author: baagaard
Date: 2006-09-04 13:14:46 -0700 (Mon, 04 Sep 2006)
New Revision: 4466

Modified:
   cs/benchmark/trunk/benchmark/ProjectionQuadPts.py
Log:
Vectorized loops in ProjectionQuadPts.py.

Modified: cs/benchmark/trunk/benchmark/ProjectionQuadPts.py
===================================================================
--- cs/benchmark/trunk/benchmark/ProjectionQuadPts.py	2006-09-02 10:11:53 UTC (rev 4465)
+++ cs/benchmark/trunk/benchmark/ProjectionQuadPts.py	2006-09-04 20:14:46 UTC (rev 4466)
@@ -157,13 +157,10 @@
 
         fieldProj = numpy.zeros((nelems, nquadpts, ndims),
                                 dtype=numpy.Float32)
+        fieldLocal = solnfield[simplices,:]
+        tmp = dot(basisfns_quad, fieldLocal)
+        fieldProj[:] = tmp.astype(numpy.Float32)
 
-        for ielem in xrange(nelems):
-            nodeIndices = simplices[ielem]
-            fieldLocal = solnfield[nodeIndices,:]
-            tmp = dot(basisfns_quad, fieldLocal)
-            fieldProj[ielem] = tmp.astype(numpy.Float32)
-
         group = tablewaiter.createGroup(self.workspace,
                                         self.workspace.root,
                                         "projection/data/%s/snapshot%d" % \
@@ -213,8 +210,6 @@
         'error'.
         """
 
-        from numpy import dot, sqrt, reshape
-
         dataA = tablewaiter.getNode(self.workspace.root,
                                     "projection/data/%s/snapshot%d/%s" % \
                                     (labelA, timestep, field))[:]
@@ -224,7 +219,7 @@
         (quadpts, quadwts) = \
                   self.discretization.getQuadPts(self.quadratureOrder)
         nquadpts = quadwts.shape[0]
-        quadwts = reshape(quadwts, (1,nquadpts))
+        quadwts = numpy.reshape(quadwts, (1,nquadpts))
         jacobian = self.workspace.root.integration.jacobian
         ndims = self.workspace.root._v_attrs.spatial_dim
 
@@ -232,13 +227,18 @@
         localerr = numpy.zeros( (nelems,), dtype=numpy.Float32)
         globalerr = 0.0
         globalvol = 0.0
-        for ielem in xrange(nelems):
-            diff = reshape(dataA[ielem] - dataB[ielem], (nquadpts, ndims))
-            localvol = jacobian[ielem] * sum(quadwts)
-            err = jacobian[ielem] * sum(sum(dot(quadwts, diff**2)))
-            localerr[ielem] = sqrt(err / localvol)
-            globalerr += err
-            globalvol += localvol
+
+        dataA = numpy.reshape(dataA, (nelems, nquadpts, ndims))
+        dataB = numpy.reshape(dataB, (nelems, nquadpts, ndims))
+        diff = dataA - dataB
+        localvol = jacobian * sum(quadwts)
+        err = jacobian * numpy.sum(numpy.sum(numpy.dot(quadwts, diff**2),
+                                             axis=2),
+                                   axis=0)
+        localerr[:] = numpy.sqrt(err / localvol).astype(numpy.Float32)
+        globalerr = sum(err)
+        globalvol = sum(localvol)
+        
         group = tablewaiter.createGroup(self.workspace,
                                         self.workspace.root,
                                         "difference/%s_%s/snapshot%d" % \
@@ -248,6 +248,7 @@
             dataset = getattr(group, field)
             dataset._f_remove()
         self.workspace.createArray(group, field, localerr)
+        from math import sqrt
         return sqrt(globalerr/globalvol)
 
 
@@ -290,7 +291,7 @@
         """
         Setup quadrature information.
         """
-        self._info.log("Seeting up integration information...")
+        self._info.log("Setting up integration information...")
 
         workspace = self.workspace
         vertices = workspace.root.geometry.vertices[:]
@@ -326,12 +327,9 @@
         ndims = len(quadpts[0])
         
         points = numpy.zeros((nelems, nquadpts, ndims), dtype=numpy.Float32)
-        for ielem in xrange(nelems):
-            j = mapFromRef[ielem]
-            origin = mapOrigin[ielem]
-            for iquad in xrange(nquadpts):   # TODO: vectorize this loop
-                pts = origin + dot(j, quadpts[iquad])
-                points[ielem, iquad] = pts.astype(numpy.Float32)
+        for iquad in xrange(nquadpts):
+            pts = mapOrigin + dot(mapFromRef, quadpts[iquad])
+            points[:,iquad] = pts.astype(numpy.Float32)
 
         group = tablewaiter.createGroup(workspace, workspace.root,
                                         "projection")