diff --git a/loopy/cse.py b/loopy/cse.py
index d9fa5656a94651392e455ca1679f4c99d9834303..45579e532110ef6dc20a26123e75faf17b45be0e 100644
--- a/loopy/cse.py
+++ b/loopy/cse.py
@@ -78,140 +78,6 @@ def to_parameters_or_project_out(param_inames, set_inames, set):
 
 
 
-def gaussian_elimination(mat, rhs):
-    m, n = mat.shape
-    i = 0
-    j = 0
-
-    from pymbolic.algorithm import lcm, gcd_many
-
-    while i < m and j < n:
-        # {{{ find pivot in column j, starting in row i
-
-        nonz_row = None
-        for k in range(i, m):
-            if mat[k,j]:
-                nonz_row = k
-                break
-
-        # }}}
-
-        if nonz_row is not None:
-            # swap rows i and nonz
-            mat[i], mat[nonz_row] = \
-                    (mat[nonz_row].copy(), mat[i].copy())
-            rhs[i], rhs[nonz_row] = \
-                    (rhs[nonz_row].copy(), rhs[i].copy())
-
-            for u in range(0, m):
-                if u == i:
-                    continue
-                if not mat[u, j]:
-                    # already 0
-                    continue
-
-                l = lcm(mat[u, j], mat[i, j])
-                u_fac = l//mat[u, j]
-                i_fac = l//mat[i, j]
-
-                mat[u] = u_fac*mat[u] - i_fac*mat[i]
-                rhs[u] = u_fac*rhs[u] - i_fac*rhs[i]
-
-                assert mat[u, j] == 0
-
-            i += 1
-
-        j += 1
-
-    for i in range(m):
-        g = gcd_many(*(
-            [a for a in mat[i] if a]
-            +
-            [a for a in rhs[i] if a]))
-
-        mat[i] //= g
-        rhs[i] //= g
-
-    return mat, rhs
-
-
-
-
-def solve_affine_equations_for(targets, equations):
-    # fix an order for targets
-    targets_list = list(targets)
-    target_idx_lut = dict((tgt_name, idx)
-            for idx, tgt_name in enumerate(targets_list))
-
-    # Find non-target variables, fix order for them
-    # Last non-target is constant.
-    nontargets = set()
-    for lhs, rhs in equations:
-        nontargets.update(get_dependencies(lhs) - targets)
-        nontargets.update(get_dependencies(rhs) - targets)
-    nontargets_list = list(nontargets)
-    nontarget_idx_lut = dict((var_name, idx)
-            for idx, var_name in enumerate(nontargets_list))
-
-    from loopy.symbolic import CoefficientCollector
-    coeff_coll = CoefficientCollector()
-
-    # {{{ build matrix and rhs
-
-    mat = np.zeros((len(equations), len(targets)), dtype=object)
-    rhs_mat = np.zeros((len(equations), len(nontargets)+1), dtype=object)
-
-    for i_eqn, (lhs, rhs) in enumerate(equations):
-        for lhs_factor, coeffs in [(1, coeff_coll(lhs)), (-1, coeff_coll(rhs))]:
-            for key, coeff in coeffs.iteritems():
-                if key in targets:
-                    mat[i_eqn, target_idx_lut[key]] = lhs_factor*coeff
-                elif key in nontargets:
-                    rhs_mat[i_eqn, nontarget_idx_lut[key]] = -lhs_factor*coeff
-                elif key == 1:
-                    rhs_mat[i_eqn, -1] = -lhs_factor*coeff
-                else:
-                    raise ValueError("key '%s' not understood" % key)
-
-    # }}}
-
-
-    mat, rhs_mat = gaussian_elimination(mat, rhs_mat)
-
-    # No need to check for overdetermined system: The case where the
-    # map is empty is treated sufficiently by isl.
-
-    result = {}
-    for j, target in enumerate(targets_list):
-        nonz_row = np.where(mat[:, j])
-        if len(nonz_row) != 1:
-            raise RuntimeError("cannot uniquely solve for '%s'" % target)
-
-        (nonz_row,), = nonz_row
-
-        if abs(mat[nonz_row, j]) != 1:
-            raise RuntimeError("division with remainder in linear solve for '%s'"
-                    % target)
-        div = mat[nonz_row, j]
-
-        target_val = int(rhs_mat[nonz_row, -1]) // div
-        for nontarget, coeff in zip(nontargets_list, rhs_mat[nonz_row]):
-            target_val += (int(coeff) // div) * var(nontarget)
-
-        result[target] = target_val
-
-    if 0:
-        for lhs, rhs in equations:
-            print lhs, '=', rhs
-        print "-------------------"
-        for lhs, rhs in result.iteritems():
-            print lhs, '=', rhs
-
-    return result
-
-
-
-
 def process_cses(kernel, uni_template,
         independent_inames, matching_vars, cse_descriptors):
     if not independent_inames:
@@ -327,19 +193,6 @@ def process_cses(kernel, uni_template,
 
         # }}}
 
-        if 0:
-            # {{{ solve for lead indices
-
-            solve_targets = set()
-            for lhs, rhs in unifier.equations:
-                solve_targets.update(get_dependencies(rhs)
-                        - kernel.non_iname_variable_names())
-
-            csed.from_unified_exprs = solve_affine_equations_for(
-                    solve_targets, unifier.equations)
-
-            # }}}
-
     return footprint, matching_var_values,