diff --git a/src/Matrix.jl b/src/Matrix.jl
index 417bf2e28..b6d2c3a62 100644
--- a/src/Matrix.jl
+++ b/src/Matrix.jl
@@ -3351,8 +3351,7 @@ triangular, and an $m\times n$ matrix $b$ over the same ring, return an
 $m\times n$ matrix $x$ such that $xU = b$. If this is not possible, an error
 will be raised.
 
-See also [`_solve_triu`](@ref) or [`can__solve_left_reduced_triu`](@ref) when
-$U$ is not square or not of full rank.
+See also [`_solve_triu`](@ref) when $U$ is not square or not of full rank.
 """
 function __solve_triu_left(b::MatElem{T}, U::MatElem{T}) where {T <: RingElement}
    n = ncols(U)
@@ -3482,64 +3481,6 @@ function is_diagonal(A::MatrixElem)
     return true
 end
 
-###############################################################################
-#
-#   Can solve
-#
-###############################################################################
-
-@doc raw"""
-    can_solve_left_reduced_triu(r::MatElem{T},
-                          M::MatElem{T}) where T <: RingElement
-Return a tuple `flag, x` where `flag` is set to true if $xM = r$ has a
-solution, where $M$ is an $m\times n$ matrix in (upper triangular) Hermite
-normal form or reduced row echelon form and $r$ and $x$ are row vectors with
-$m$ columns (i.e. $1 \times m$ matrices). If there is no solution, flag is set
-to `false` and $x$ is set to zero.
-"""
-function _can_solve_left_reduced_triu(r::MatElem{T},
-                          M::MatElem{T}) where T <: RingElement
-   ncols(r) != ncols(M) && error("Incompatible matrices")
-   r = deepcopy(r) # do not destroy input
-   m = ncols(r)
-   n = nrows(M)
-   if n == 0
-      return true, r
-   end
-   R = base_ring(r)
-   x = zero_matrix(R, 1, n)
-   j = 1 # row in M
-   k = 1 # column in M
-   t = R()
-   for i = 1:m # column in r
-      if is_zero_entry(r, 1, i)
-         continue
-      end
-      while k <= i && j <= n
-         if is_zero_entry(M, j, k)
-            k += 1
-         elseif k < i
-            j += 1
-         else
-            break
-         end
-      end
-      if k != i
-         return false, x
-      end
-      x[1, j], r[1, i] = divrem(r[1, i], M[j, k])
-      if !is_zero_entry(r, 1, i)
-         return false, x
-      end
-      q = -x[1, j]
-      for l = i + 1:m
-         t = mul!(t, q, M[j, l])
-         r[1, l] = addeq!(r[1, l], t)
-      end
-   end
-   return true, x
-end
-
 ###############################################################################
 #
 #   Inverse
diff --git a/src/Module.jl b/src/Module.jl
index 32fd96c28..3befa97de 100644
--- a/src/Module.jl
+++ b/src/Module.jl
@@ -223,20 +223,17 @@ function ==(M::FPModule{T}, N::FPModule{T}) where T <: RingElement
          mat2[i + r2, j] = prels[i][1, j]
       end
    end
-   # Put the matrices into reduced form
-   mat1 = reduced_form(mat1)
-   mat2 = reduced_form(mat2)
+   sol_ctx1 = solve_init(mat1)
+   sol_ctx2 = solve_init(mat2)
    # Check containment of rewritten gens of M in row space of mat2
    for v in G1
-      flag, r = _can_solve_left_reduced_triu(Generic._matrix(v), mat2)
-      if !flag
+      if !can_solve(sol_ctx2, Generic._matrix(v))
          return false
       end
    end
    # Check containment of rewritten gens of N in row space of mat1
    for v in G2
-      flag, r = _can_solve_left_reduced_triu(Generic._matrix(v), mat1)
-      if !flag
+      if !can_solve(sol_ctx1, Generic._matrix(v))
          return false
       end
    end
diff --git a/test/generic/Matrix-test.jl b/test/generic/Matrix-test.jl
index fd7853df5..024662eb3 100644
--- a/test/generic/Matrix-test.jl
+++ b/test/generic/Matrix-test.jl
@@ -2639,24 +2639,6 @@ end
    end
 end
 
-@testset "Generic.Mat._solve_left_reduced_triu" begin
-   for iter = 1:40
-      n = rand(1:6)
-      m = rand(1:n)
-      S = matrix_space(ZZ, m, n)
-      U = matrix_space(ZZ, 1, n)
-
-      M = randmat_with_rank(S, rand(1:m), -20:20)
-      r = rand(U, -20:20)
-
-      M = hnf(M)
-
-      flag, x = AbstractAlgebra._can_solve_left_reduced_triu(r, M)
-
-      @test flag == false || x*M == r
-   end
-end
-
 @testset "Generic.Mat.rref" begin
    # Non-integral domain