Copy solve context + RREFTrait from Nemo

src/Solve.jl

@@ -249,6 +249,22 @@ matrix_normal_form_type(C::SolveCtx{T, NF}) where {T, NF} = NF()
 
 matrix(C::SolveCtx) = C.A
 
+function _init_reduce(C::SolveCtx{T, RREFTrait}) where T
+  if isdefined(C, :red) && isdefined(C, :trafo)
+    return nothing
+  end
+
+  B, U = echelon_form_with_transformation(matrix(C))
+  r = nrows(B)
+  while r > 0 && is_zero_row(B, r)
+    r -= 1
+  end
+  set_rank!(C, r)
+  C.red = B
+  C.trafo = U
+  return nothing
+end
+
 function _init_reduce(C::SolveCtx{T, LUTrait}) where T
   if isdefined(C, :red) && isdefined(C, :lu_perm)
     return nothing
@@ -327,6 +343,22 @@ function _init_reduce(C::SolveCtx{T, HowellFormTrait}) where T
   return nothing
 end
 
+function _init_reduce_transpose(C::SolveCtx{T, RREFTrait}) where T
+  if isdefined(C, :red_transp) && isdefined(C, :trafo_transp)
+    return nothing
+  end
+
+  B, U = echelon_form_with_transformation(lazy_transpose(matrix(C)))
+  r = nrows(B)
+  while r > 0 && is_zero_row(B, r)
+    r -= 1
+  end
+  set_rank!(C, r)
+  C.red_transp = B
+  C.trafo_transp = U
+  return nothing
+end
+
 function _init_reduce_transpose(C::SolveCtx{T, LUTrait}) where T
   if isdefined(C, :red_transp) && isdefined(C, :lu_perm_transp)
     return nothing
@@ -968,7 +1000,26 @@ function _can_solve_internal_no_check(::RREFTrait, A::MatElem{T}, b::MatElem{T},
   return true, sol, X
 end
 
-# RREFTrait with SolveCtx is not implemented (it is in Nemo)
+function _can_solve_internal_no_check(::RREFTrait, C::SolveCtx{T, RREFTrait}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T
+  if side === :right
+    fl, sol = _can_solve_with_rref(b, transformation_matrix(C), rank(C), pivot_and_non_pivot_cols(C), task)
+    if !fl || task !== :with_kernel
+      return fl, sol, zero(b, 0, 0)
+    end
+
+    _, K = _kernel_of_rref(reduced_matrix(C), rank(C), pivot_and_non_pivot_cols(C))
+    return fl, sol, K
+  else# side === :left
+    fl, sol_transp = _can_solve_with_rref(lazy_transpose(b), transformation_matrix_of_transpose(C), rank(C), pivot_and_non_pivot_cols_of_transpose(C), task)
+    sol = lazy_transpose(sol_transp)
+    if !fl || task !== :with_kernel
+      return fl, sol, zero(b, 0, 0)
+    end
+
+    _, K = _kernel_of_rref(reduced_matrix_of_transpose(C), rank(C), pivot_and_non_pivot_cols_of_transpose(C))
+    return fl, sol, lazy_transpose(K)
+  end
+end
 
 ### LUTrait
 

test/Solve-test.jl

@@ -91,18 +91,32 @@
   @test K == identity_matrix(R, 2) || K == swap_cols!(identity_matrix(R, 2), 1, 2)
 end
 
-@testset "Linear solving context over $R" for R in [ QQ, ZZ, GF(101), fraction_field(QQ["x"][1]), residue_ring(ZZ, 16)[1] ]
+@testset "Linear solving context over $R with $NFTrait" for (R, NFTrait, is_default) in [
+    (QQ, AbstractAlgebra.Solve.FFLUTrait, true),
+    (ZZ, AbstractAlgebra.Solve.HermiteFormTrait, true),
+    (GF(101), AbstractAlgebra.Solve.LUTrait, true),
+    (GF(101), AbstractAlgebra.Solve.RREFTrait, false),
+    (fraction_field(QQ["x"][1]), AbstractAlgebra.Solve.FFLUTrait, true),
+    (residue_ring(ZZ, 16)[1], AbstractAlgebra.Solve.HowellFormTrait, true)
+   ]
   M = matrix(R, [1 2 3 4 5; 0 0 8 9 10; 0 0 0 14 15])
-  C = AbstractAlgebra.Solve.solve_init(M)
-
-  @test C isa AbstractAlgebra.solve_context_type(R)
-  @test C isa AbstractAlgebra.solve_context_type(M)
-  @test C isa AbstractAlgebra.solve_context_type(AbstractAlgebra.Solve.matrix_normal_form_type(C), elem_type(R))
-  @test C isa AbstractAlgebra.solve_context_type(AbstractAlgebra.Solve.matrix_normal_form_type(C), R(1))
-  @test C isa AbstractAlgebra.solve_context_type(AbstractAlgebra.Solve.matrix_normal_form_type(C), typeof(R))
-  @test C isa AbstractAlgebra.solve_context_type(AbstractAlgebra.Solve.matrix_normal_form_type(C), R)
-  @test C isa AbstractAlgebra.solve_context_type(AbstractAlgebra.Solve.matrix_normal_form_type(C), typeof(M))
-  @test C isa AbstractAlgebra.solve_context_type(AbstractAlgebra.Solve.matrix_normal_form_type(C), M)
+
+  if is_default
+    C = solve_init(M)
+    @test AbstractAlgebra.Solve.matrix_normal_form_type(C) === NFTrait()
+    @test C isa AbstractAlgebra.solve_context_type(R)
+    @test C isa AbstractAlgebra.solve_context_type(M)
+  end
+
+  C = solve_init(NFTrait(), M)
+
+  @test AbstractAlgebra.Solve.matrix_normal_form_type(C) === NFTrait()
+  @test C isa AbstractAlgebra.solve_context_type(NFTrait(), elem_type(R))
+  @test C isa AbstractAlgebra.solve_context_type(NFTrait(), R(1))
+  @test C isa AbstractAlgebra.solve_context_type(NFTrait(), typeof(R))
+  @test C isa AbstractAlgebra.solve_context_type(NFTrait(), R)
+  @test C isa AbstractAlgebra.solve_context_type(NFTrait(), typeof(M))
+  @test C isa AbstractAlgebra.solve_context_type(NFTrait(), M)
 
   @test_throws ErrorException AbstractAlgebra.Solve.solve(C, [ R(1) ])
   @test_throws ErrorException AbstractAlgebra.Solve.solve(C, [ R(1) ], side = :right)