diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 3be8173c10a9b7de5e03f6444b3159b16b108407..70c3c1a9498feb57dee299d5b5f9815838df3e35 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -8,9 +8,10 @@ Python 2.7:
   - python2.7
   except:
   - tags
-Python 3.5:
+
+Python 3.6:
   script:
-  - py_version=3.5
+  - py_version=3.6
   - EXTRA_INSTALL=numpy
   - 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"
@@ -18,6 +19,7 @@ Python 3.5:
   - python3.5
   except:
   - tags
+
 Python 2.6:
   script:
   - py_version=2.6
@@ -28,6 +30,7 @@ Python 2.6:
   - python2.6
   except:
   - tags
+
 Documentation:
   script:
   - EXTRA_INSTALL="numpy"
@@ -37,3 +40,12 @@ Documentation:
   - python3.5
   only:
   - master
+
+Flake8:
+  script:
+  - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/prepare-and-run-flake8.sh
+  - ". ./prepare-and-run-flake8.sh pyopencl test"
+  tags:
+  - python3.5
+  except:
+  - tags
diff --git a/modepy/__init__.py b/modepy/__init__.py
index 0c180587b526029bddadf9f1d2d872325421e85b..72ee4090515f5d39d4e1242b7782b8e84c55c250 100644
--- a/modepy/__init__.py
+++ b/modepy/__init__.py
@@ -40,7 +40,7 @@ from modepy.quadrature.xiao_gimbutas import XiaoGimbutasSimplexQuadrature
 from modepy.quadrature.vioreanu_rokhlin import VioreanuRokhlinSimplexQuadrature
 from modepy.quadrature.grundmann_moeller import GrundmannMoellerSimplexQuadrature
 
-from modepy.version import VERSION_TEXT as __version__
+from modepy.version import VERSION_TEXT as __version__  # noqa: N811
 
 __all__ = [
         "__version__",
diff --git a/modepy/modal_decay.py b/modepy/modal_decay.py
index d56c30adaf85781e9cb719d19681b4f0fe87dba1..9004123baf80ce42e25ce079826193c43e1cae5e 100644
--- a/modepy/modal_decay.py
+++ b/modepy/modal_decay.py
@@ -42,8 +42,6 @@ The method implemented in this module follows this article:
 """
 
 
-
-
 def make_mode_number_vector(mode_order_tuples, ignored_modes):
     node_cnt = len(mode_order_tuples)
 
@@ -56,7 +54,6 @@ def make_mode_number_vector(mode_order_tuples, ignored_modes):
     return mode_number_vector
 
 
-
 def create_decay_baseline(mode_number_vector, n):
     """Create a vector of modal coefficients that exhibit 'optimal'
     (:math:`k^{-N}`) decay.
@@ -65,29 +62,25 @@ def create_decay_baseline(mode_number_vector, n):
     zeros = mode_number_vector == 0
 
     modal_coefficients = mode_number_vector**(-n)
-    modal_coefficients[zeros] = 1 # irrelevant, just keeps log from NaNing
+    modal_coefficients[zeros] = 1  # irrelevant, just keeps log from NaNing
 
     modal_coefficients /= la.norm(modal_coefficients)
 
     return modal_coefficients
 
 
-
-
 def get_decay_fit_matrix(mode_number_vector, ignored_modes, weight_vector):
     a = np.zeros((len(mode_number_vector), 2), dtype=np.float64)
-    a[:,0] = weight_vector
-    a[:,1] = weight_vector * np.log(mode_number_vector)
+    a[:, 0] = weight_vector
+    a[:, 1] = weight_vector * np.log(mode_number_vector)
 
     if ignored_modes == 0:
-        assert not np.isfinite(a[0,1])
-        a[0,1] = 0
+        assert not np.isfinite(a[0, 1])
+        a[0, 1] = 0
 
     return la.pinv(a)
 
 
-
-
 def skyline_pessimize(modal_values):
     nelements, nmodes = modal_values.shape
 
@@ -106,8 +99,6 @@ def skyline_pessimize(modal_values):
     return result
 
 
-
-
 def fit_modal_decay(coeffs, dims, n, ignored_modes=1):
     """Fit a curve to the modal decay on each element.
 
@@ -117,7 +108,7 @@ def fit_modal_decay(coeffs, dims, n, ignored_modes=1):
     :arg ignored_modes: the number of modal coefficients to ignore at the
         beginning. The default value of '1' ignores the constant.
     :returns: a tuple *(expn, constant)* of arrays of length *num_elements*,
-        where the modal decay is fit to the curve 
+        where the modal decay is fit to the curve
         ``constant * total_degree**exponent``.
 
     ``-exponent-1`` can be used as a rough indicator of how many continuous
@@ -131,7 +122,8 @@ def fit_modal_decay(coeffs, dims, n, ignored_modes=1):
 
     coeffs_squared = skyline_pessimize(coeffs**2)
 
-    mode_number_vector_int = make_mode_number_vector(mode_order_tuples, ignored_modes)
+    mode_number_vector_int = make_mode_number_vector(
+            mode_order_tuples, ignored_modes)
     mode_number_vector = mode_number_vector_int.astype(np.float64)
     weight_vector = np.ones_like(mode_number_vector)
 
@@ -144,7 +136,7 @@ def fit_modal_decay(coeffs, dims, n, ignored_modes=1):
             * el_norm_squared[:, np.newaxis])**2
     log_modal_coeffs = np.log(coeffs_squared[:, ignored_modes:] + scaled_baseline)/2
 
-    assert fit_mat.shape[0] == 2 # exponent and log(constant)
+    assert fit_mat.shape[0] == 2  # exponent and log(constant)
 
     fit_values = np.dot(fit_mat, (weight_vector*log_modal_coeffs).T).T
 
@@ -155,27 +147,26 @@ def fit_modal_decay(coeffs, dims, n, ignored_modes=1):
         import matplotlib.pyplot as pt
         pt.plot(log_modal_coeffs.flat, "o-")
 
-        fit = np.log(const[:, np.newaxis] * mode_number_vector**exponent[:, np.newaxis])
+        fit = np.log(const[:, np.newaxis]
+                * mode_number_vector**exponent[:, np.newaxis])
 
         pt.plot(fit.flat)
 
-        #plot_expt = np.zeros_like(log_modal_coeffs)
-        #plot_expt[:] = exponent[:, np.newaxis]
-        #pt.plot(plot_expt.flat)
+        # plot_expt = np.zeros_like(log_modal_coeffs)
+        # plot_expt[:] = exponent[:, np.newaxis]
+        # pt.plot(plot_expt.flat)
 
         pt.show()
 
     return exponent, const
 
 
-
-
 def estimate_relative_expansion_residual(coeffs, dims, n, ignored_modes=1):
     """Use the modal fit to estimate the relative residual of the expansion.
     The arguments to this function exactly match :func:`fit_modal_decay`.
 
     :returns: An array of estimates of the fraction of the :math:`L^2` norm
-        contained in the (unrepresented) tail of 
+        contained in the (unrepresented) tail of
 
     An idea like this is described in this article:
 
diff --git a/modepy/modes.py b/modepy/modes.py
index 0cd5e763c6fb0f501bf8b593c0b9fb1fc68ad1c5..e6dfcc197777e153a5dfe2d1c658db9f21a8f864 100644
--- a/modepy/modes.py
+++ b/modepy/modes.py
@@ -120,12 +120,12 @@ def jacobi(alpha, beta, n, x):
     from modepy.tools import gamma
 
     n = np.int32(n)
-    Nx = len(x)
+    Nx = len(x)  # noqa: N806
     if x.shape[0] > 1:
         x = x.T
 
     # Storage for recursive construction
-    PL = np.zeros((Nx, n+1), np.float64)
+    PL = np.zeros((Nx, n+1), np.float64)  # noqa: N806
 
     # Initial values P_0(x) and P_1(x)
     gamma0 = (2**(alpha+beta+1)
@@ -266,7 +266,7 @@ def grad_pkdo_2d(order, rs):
 # {{{ 3D PKDO
 
 def _rsttoabc(r, s, t):
-    Np = len(r)
+    Np = len(r)  # noqa: N806
     tol = 1e-10
 
     a = np.zeros(Np)
diff --git a/modepy/quadrature/grundmann_moeller.py b/modepy/quadrature/grundmann_moeller.py
index 006abfc0b640c86ae5582e5b5c6435244b73a773..09c2a73c3606faad34562d863ba885c4f09b9776 100644
--- a/modepy/quadrature/grundmann_moeller.py
+++ b/modepy/quadrature/grundmann_moeller.py
@@ -42,12 +42,12 @@ def _extended_euclidean(q, r):
         p, a, b = _extended_euclidean(r, q)
         return p, b, a
 
-    Q = 1, 0
-    R = 0, 1
+    Q = 1, 0  # noqa: N806
+    R = 0, 1  # noqa: N806
 
     while r:
         quot, t = divmod(q, r)
-        T = Q[0] - quot*R[0], Q[1] - quot*R[1]
+        T = Q[0] - quot*R[0], Q[1] - quot*R[1]  # noqa: N806
         q, r = r, t
         Q, R = R, T
 
diff --git a/modepy/quadrature/vr_quad_data_tet.py b/modepy/quadrature/vr_quad_data_tet.py
index 5a1311e31d10e687255ceace60122ae14602cf2e..df8ba46ce5e11c677fdba68b51f6303b28c8fd52 100644
--- a/modepy/quadrature/vr_quad_data_tet.py
+++ b/modepy/quadrature/vr_quad_data_tet.py
@@ -9,6 +9,7 @@ points = "points"
 weights = "weights"
 quad_degree = "quad_degree"
 
+
 def process_rule(table):
     result = {}
     for order, kv in six.iteritems(table):
@@ -22,7 +23,8 @@ def process_rule(table):
 
     return result
 
-tetrahedron_data = process_rule({0: {points: [[0.00000000000000000000000000000000e+00],
+tetrahedron_data = process_rule({  # noqa
+ 0: {points: [[0.00000000000000000000000000000000e+00],
                 [0.00000000000000000000000000000000e+00],
                 [0.00000000000000000000000000000000e+00]],
      quad_degree: 1,
diff --git a/modepy/quadrature/vr_quad_data_tri.py b/modepy/quadrature/vr_quad_data_tri.py
index b83b503dc29eed347df07979d1d13d6d53bcf947..5c58da5013363ad0af4514c19f5b0cdf4a776e18 100644
--- a/modepy/quadrature/vr_quad_data_tri.py
+++ b/modepy/quadrature/vr_quad_data_tri.py
@@ -9,6 +9,7 @@ points = "points"
 weights = "weights"
 quad_degree = "quad_degree"
 
+
 def process_rule(table):
     result = {}
     for order, kv in six.iteritems(table):
@@ -22,7 +23,8 @@ def process_rule(table):
 
     return result
 
-triangle_data = process_rule({0: {points: [[0.00000000000000000000000000000000e+00],
+triangle_data = process_rule({  # noqa
+    0: {points: [[0.00000000000000000000000000000000e+00],
                 [0.00000000000000000000000000000000e+00],
                 []],
      quad_degree: 1,
diff --git a/modepy/quadrature/xg_quad_data.py b/modepy/quadrature/xg_quad_data.py
index 20e04123f5e087b348f9ba55f32dc3ecaf77ef5f..35dff94e70e3655b7f9f67412bf25d6016cd51c6 100644
--- a/modepy/quadrature/xg_quad_data.py
+++ b/modepy/quadrature/xg_quad_data.py
@@ -5,7 +5,8 @@ from __future__ import absolute_import
 import numpy
 import six
 
-triangle_table = {1: {'points': [[-5.55111512312578e-17, 2.22044604925031e-16]],
+triangle_table = {  # noqa
+ 1: {'points': [[-5.55111512312578e-17, 2.22044604925031e-16]],
      'weights': [1.73205080756888]},
  2: {'points': [[-0.5, -0.288675134594813],
                 [0.5, -0.288675134594813],
@@ -16114,7 +16115,8 @@ triangle_table = dict(
     for order, rule in six.iteritems(triangle_table))
 
 
-tetrahedron_table = {1: {'points': [[0.0, 0.0, 0.0]], 'weights': [0.942809041582063]},
+tetrahedron_table = {  # noqa
+ 1: {'points': [[0.0, 0.0, 0.0]], 'weights': [0.942809041582063]},
  2: {'points': [[-0.267702858591007, 0.0593901700619949, -0.396942694114215],
                 [0.151093184153314, 0.435408773247631, 0.215171925430692],
                 [-0.136769919523739, -0.328034111559041, 0.295084651993713],
@@ -18786,4 +18788,3 @@ tetrahedron_table = {1: {'points': [[0.0, 0.0, 0.0]], 'weights': [0.942809041582
 tetrahedron_table = dict(
     (order, dict((name, numpy.array(ary)) for name, ary in six.iteritems(rule)))
     for order, rule in six.iteritems(tetrahedron_table))
-
diff --git a/modepy/quadrature/xiao_gimbutas.py b/modepy/quadrature/xiao_gimbutas.py
index 4db6ac0e362cf849e3d803fe191c9668b318ce91..ac73eefd90e0c231ea885c16d62a368588756c53 100644
--- a/modepy/quadrature/xiao_gimbutas.py
+++ b/modepy/quadrature/xiao_gimbutas.py
@@ -23,13 +23,9 @@ THE SOFTWARE.
 """
 
 
-
-
 from modepy.quadrature import Quadrature, QuadratureRuleUnavailable
 
 
-
-
 class XiaoGimbutasSimplexQuadrature(Quadrature):
     """A (nearly) Gaussian simplicial quadrature with very few quadrature nodes,
     available for low-to-moderate orders.
diff --git a/test/test_nodes.py b/test/test_nodes.py
index 79412d42f12c9c130be6a493521924ef3e58a3a6..7dd2fb9f5cb165feef087f177da348fbc2df3802 100644
--- a/test/test_nodes.py
+++ b/test/test_nodes.py
@@ -23,15 +23,11 @@ THE SOFTWARE.
 """
 
 
-
-
 import numpy as np
 import numpy.linalg as la
 import pytest
 
 
-
-
 @pytest.mark.parametrize("dims", [1, 2, 3])
 def test_barycentric_coordinate_map(dims):
     from random import Random
@@ -60,8 +56,6 @@ def test_barycentric_coordinate_map(dims):
     assert la.norm(unit-unit3) < 1e-14
 
 
-
-
 def test_warp():
     """Check some assumptions on the node warp factor calculator"""
     n = 17
@@ -81,8 +75,6 @@ def test_warp():
     assert abs(lgq(wfc)) < 6e-14
 
 
-
-
 def test_tri_face_node_distribution():
     """Test whether the nodes on the faces of the triangle are distributed
     according to the same proportions on each face.
@@ -114,14 +106,12 @@ def test_tri_face_node_distribution():
         pfp = np.array([np.dot(direction, unodes[:, i]-start) for i in face_i])
         projected_face_points.append(pfp)
 
-    first_points =  projected_face_points[0]
+    first_points = projected_face_points[0]
     for points in projected_face_points[1:]:
         error = la.norm(points-first_points, np.Inf)
         assert error < 1e-15
 
 
-
-
 @pytest.mark.parametrize("dims", [1, 2, 3])
 @pytest.mark.parametrize("n", [1, 3, 6])
 def test_simp_nodes(dims, n):
@@ -135,14 +125,11 @@ def test_simp_nodes(dims, n):
     assert (np.sum(unodes) <= eps).all()
 
 
-
-
-
 def test_affine_map():
     """Check that our cheapo geometry-targeted linear algebra actually works."""
     from modepy.tools import AffineMap
+    # for d in [3]:
     for d in range(1, 5):
-    #for d in [3]:
         for i in range(100):
             a = np.random.randn(d, d)+10*np.eye(d)
             b = np.random.randn(d)
@@ -153,8 +140,6 @@ def test_affine_map():
             assert la.norm(x-m.inverse(m(x))) < 1e-10
 
 
-
-
 # You can test individual routines by typing
 # $ python test_nodes.py 'test_routine()'