diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 9d960a6f4eb920f920fb7398c1816e447e97e05e..6ee07ae45d90f7fa66cf6d8ca7e219d5cde2017b 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -53,6 +53,19 @@ Python 3.5 POCL:
   except:
   - tags
 
+Python 3.6 POCL:
+  script:
+  - export PY_EXE=python3.6
+  - export PYOPENCL_TEST=portable
+  - export EXTRA_INSTALL="numpy mako"
+  - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh
+  - ". ./build-and-test-py-project.sh"
+  tags:
+  - python3.6
+  - pocl
+  except:
+  - tags
+
 Documentation:
   script:
   - EXTRA_INSTALL="numpy mako"
diff --git a/examples/curve-pot.py b/examples/curve-pot.py
index 9d8678f5e6d962c5656fa572c06d48b3a167a932..f42122088dbb21c12b589804566618fd35e1ad7c 100644
--- a/examples/curve-pot.py
+++ b/examples/curve-pot.py
@@ -229,6 +229,7 @@ def draw_pot_figure(aspect_ratio,
     else:
         # {{{ 3D plots
 
+
         plotval_vol = vol_pot.real
         plotval_c = curve_pot.real
 
@@ -264,7 +265,8 @@ def draw_pot_figure(aspect_ratio,
 if __name__ == "__main__":
     draw_pot_figure(aspect_ratio=1, nsrc=100, novsmp=100, helmholtz_k=(15+4j)*0.3,
             what_operator="D", what_operator_lpot="D", force_center_side=1)
-    pt.savefig("eigvals-ext-nsrc100-novsmp100.pdf")
+
+#    pt.savefig("eigvals-ext-nsrc100-novsmp100.pdf")
     #pt.clf()
     #draw_pot_figure(aspect_ratio=1, nsrc=100, novsmp=100, helmholtz_k=0,
     #        what_operator="D", what_operator_lpot="D", force_center_side=-1)
diff --git a/sumpy/cse.py b/sumpy/cse.py
index 3318fba023eccbc7673fab31db5d31b0cfeb4035..573c78c93961c8fe0e69c9650c4df7d10133e98b 100644
--- a/sumpy/cse.py
+++ b/sumpy/cse.py
@@ -187,11 +187,32 @@ class FuncArgTracker(object):
         from collections import defaultdict
         count_map = defaultdict(lambda: 0)
 
-        for arg in argset:
-            for func_i in self.arg_to_funcset[arg]:
+        # Sorted by size to make best use of the performance hack below.
+        funcsets = sorted((self.arg_to_funcset[arg] for arg in argset), key=len)
+
+        for funcset in funcsets[:-threshold+1]:
+            for func_i in funcset:
                 if func_i >= min_func_i:
                     count_map[func_i] += 1
 
+        for i, funcset in enumerate(funcsets[-threshold+1:]):
+            # When looking at the tail end of the funcsets list, items below
+            # this threshold in the count_map don't have to be considered
+            # because they can't possibly be in the output.
+            count_map_threshold = i + 1
+
+            # We pick the smaller of the two containers to iterate over to
+            # reduce the number of items we have to look at.
+            (smaller_funcs_container,
+             larger_funcs_container) = sorted([funcset, count_map], key=len)
+
+            for func_i in smaller_funcs_container:
+                if count_map[func_i] < count_map_threshold:
+                    continue
+
+                if func_i in larger_funcs_container:
+                    count_map[func_i] += 1
+
         return dict(
             (k, v) for k, v in count_map.items()
             if v >= threshold)
@@ -275,14 +296,14 @@ def match_common_args(func_class, funcs, opt_subs):
     from sumpy.tools import OrderedSet
 
     for i in range(len(funcs)):
-        common_arg_candidates = arg_tracker.get_common_arg_candidates(
+        common_arg_candidates_counts = arg_tracker.get_common_arg_candidates(
                 arg_tracker.func_to_argset[i], i + 1, threshold=2)
 
         # Sort the candidates in order of match size.
         # This makes us try combining smaller matches first.
         common_arg_candidates = OrderedSet(sorted(
-                common_arg_candidates.keys(),
-                key=lambda k: (common_arg_candidates[k], k)))
+                common_arg_candidates_counts.keys(),
+                key=lambda k: (common_arg_candidates_counts[k], k)))
 
         while common_arg_candidates:
             j = common_arg_candidates.pop(last=False)
diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index 0121cf6bcd198e0c04821d06b7110020e8428cc7..e29ac5167051e5d781a9fdb331d97fb3fef84df2 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -27,7 +27,6 @@ from six.moves import range, zip
 import loopy as lp
 import numpy as np
 from pymbolic.mapper import IdentityMapper, CSECachingMapperMixin
-import sumpy.symbolic as sym
 from sumpy.symbolic import pymbolic_real_norm_2
 from pymbolic.primitives import make_sym_vector
 from pymbolic import var
@@ -51,6 +50,7 @@ PDE kernels
 .. autoclass:: LaplaceKernel
 .. autoclass:: HelmholtzKernel
 .. autoclass:: StokesletKernel
+.. autoclass:: StressletKernel
 
 Derivatives
 -----------
@@ -264,11 +264,10 @@ class ExpressionKernel(Kernel):
         if self.dim != len(dist_vec):
             raise ValueError("dist_vec length does not match expected dimension")
 
-        import sumpy.symbolic as sym
-        expr = expr.subs(dict(
-            (sym.Symbol("d%d" % i), dist_vec_i)
+        expr = expr.subs([
+            ("d%d" % i, dist_vec_i)
             for i, dist_vec_i in enumerate(dist_vec)
-            ))
+            ])
 
         return expr
 
@@ -429,7 +428,7 @@ class StokesletKernel(ExpressionKernel):
     init_arg_names = ("dim", "icomp", "jcomp", "viscosity_mu_name")
 
     def __init__(self, dim, icomp, jcomp, viscosity_mu_name="mu"):
-        """
+        r"""
         :arg viscosity_mu_name: The argument name to use for
                 dynamic viscosity :math:`\mu` the then generating functions to
                 evaluate this kernel.
@@ -494,12 +493,11 @@ class StokesletKernel(ExpressionKernel):
 
 
 class StressletKernel(ExpressionKernel):
-    init_arg_names = ("dim", "icomp", "jcomp", "viscosity_mu_name",
-        "stresslet_vector_name")
+    init_arg_names = ("dim", "icomp", "jcomp", "kcomp", "viscosity_mu_name")
 
-    def __init__(self, dim=None, icomp=None, jcomp=None, viscosity_mu_name="mu",
-                        stresslet_vector_name="stresslet_vec"):
-        """
+    def __init__(self, dim=None, icomp=None, jcomp=None, kcomp=None,
+                        viscosity_mu_name="mu"):
+        r"""
         :arg viscosity_mu_name: The argument name to use for
                 dynamic viscosity :math:`\mu` the then generating functions to
                 evaluate this kernel.
@@ -508,23 +506,17 @@ class StressletKernel(ExpressionKernel):
 
         if dim == 2:
             d = make_sym_vector("d", dim)
-            n = make_sym_vector(stresslet_vector_name, dim)
             r = pymbolic_real_norm_2(d)
             expr = (
-                sum(n[axis]*d[axis] for axis in range(dim))
-                *
-                d[icomp]*d[jcomp]/r**4
+                d[icomp]*d[jcomp]*d[kcomp]/r**4
                 )
             scaling = 1/(var("pi"))
 
         elif dim == 3:
             d = make_sym_vector("d", dim)
-            n = make_sym_vector(stresslet_vector_name, dim)
             r = pymbolic_real_norm_2(d)
             expr = (
-                sum(n[axis]*d[axis] for axis in range(dim))
-                *
-                d[icomp]*d[jcomp]/r**5
+                d[icomp]*d[jcomp]*d[kcomp]/r**5
                 )
             scaling = -3/(4*var("pi"))
 
@@ -535,9 +527,9 @@ class StressletKernel(ExpressionKernel):
             raise RuntimeError("unsupported dimensionality")
 
         self.viscosity_mu_name = viscosity_mu_name
-        self.stresslet_vector_name = stresslet_vector_name
         self.icomp = icomp
         self.jcomp = jcomp
+        self.kcomp = kcomp
 
         ExpressionKernel.__init__(
                 self,
@@ -547,42 +539,26 @@ class StressletKernel(ExpressionKernel):
                 is_complex_valued=False)
 
     def __getinitargs__(self):
-        return (self._dim, self.icomp, self.jcomp, self.viscosity_mu_name,
-                      self.stresslet_vector_name)
+        return (self._dim, self.icomp, self.jcomp, self.kcomp,
+                      self.viscosity_mu_name)
 
     def update_persistent_hash(self, key_hash, key_builder):
         key_hash.update(type(self).__name__.encode())
         key_builder.rec(key_hash, (
-            self.dim, self.icomp, self.jcomp, self.viscosity_mu_name,
-            self.stresslet_vector_name))
+            self.dim, self.icomp, self.jcomp, self.kcomp,
+            self.viscosity_mu_name))
 
     def __repr__(self):
-        return "StressletKnl%dD_%d%d[%s]" % (self.dim, self.icomp, self.jcomp,
-                self.stresslet_vector_name)
+        return "StressletKnl%dD_%d%d%d" % (self.dim, self.icomp, self.jcomp,
+                self.kcomp)
 
     def get_args(self):
         return [
                 KernelArgument(
                     loopy_arg=lp.ValueArg(self.viscosity_mu_name, np.float64),
-                    ),
-                KernelArgument(
-                        loopy_arg=lp.GlobalArg(self.stresslet_vector_name,
-                            None,
-                            shape=(self.dim, "nsources"),
-                            dim_tags="sep,C"))
+                    )
                 ]
 
-    def get_code_tranformer(self):
-        from sumpy.codegen import VectorComponentRewriter
-        vcr = VectorComponentRewriter([self.stresslet_vector_name])
-        from pymbolic.primitives import Variable
-        via = _VectorIndexAdder(self.stresslet_vector_name, (Variable("isrc"),))
-
-        def transform(expr):
-            return via(vcr(expr))
-
-        return transform
-
     mapper_method = "map_stresslet_kernel"
 
 # }}}
@@ -739,12 +715,12 @@ class DirectionalTargetDerivative(DirectionalDerivative):
         dim = len(bvec)
         assert dim == self.dim
 
-        from sumpy.symbolic import make_sym_vector
-        dir_vec = make_sym_vector(self.dir_vec_name, dim)
+        from sumpy.symbolic import make_sympy_vector
+        dir_vec = make_sympy_vector(self.dir_vec_name, dim)
 
         # bvec = tgt-center
-        return sym.Add(*tuple(dir_vec[axis]*expr.diff(bvec[axis])
-                for axis in range(dim)))
+        return sum(dir_vec[axis]*expr.diff(bvec[axis])
+                for axis in range(dim))
 
     mapper_method = "map_directional_target_derivative"
 
@@ -770,12 +746,12 @@ class DirectionalSourceDerivative(DirectionalDerivative):
         dimensions = len(avec)
         assert dimensions == self.dim
 
-        from sumpy.symbolic import make_sym_vector
-        dir_vec = make_sym_vector(self.dir_vec_name, dimensions)
+        from sumpy.symbolic import make_sympy_vector
+        dir_vec = make_sympy_vector(self.dir_vec_name, dimensions)
 
         # avec = center-src -> minus sign from chain rule
-        return sym.Add(*tuple(-dir_vec[axis]*expr.diff(avec[axis])
-                for axis in range(dimensions)))
+        return sum(-dir_vec[axis]*expr.diff(avec[axis])
+                for axis in range(dimensions))
 
     mapper_method = "map_directional_source_derivative"
 
@@ -909,8 +885,8 @@ class KernelDimensionSetter(KernelIdentityMapper):
         return StressletKernel(self.dim,
                 kernel.icomp,
                 kernel.jcomp,
-                viscosity_mu_name=kernel.viscosity_mu_name,
-                stresslet_vector_name=kernel.stresslet_vector_name)
+                kernel.kcomp,
+                viscosity_mu_name=kernel.viscosity_mu_name)
 
 # }}}