From 088d0c64daffbd1fe5b062e47880b50059e1307d Mon Sep 17 00:00:00 2001 From: Johannes Schmitt Date: Wed, 28 Feb 2024 13:50:39 +0100 Subject: [PATCH 01/10] Remove outdated `_can_solve_with_kernel` --- src/Matrix.jl | 141 ------------------------------------ test/generic/Matrix-test.jl | 10 +-- 2 files changed, 5 insertions(+), 146 deletions(-) diff --git a/src/Matrix.jl b/src/Matrix.jl index bd76400e08..bd5d7ed4d5 100644 --- a/src/Matrix.jl +++ b/src/Matrix.jl @@ -3204,147 +3204,6 @@ function find_pivot(A::MatElem{T}) where T <: RingElement return p end -############################################################################### -# -# Solving with kernel -# -############################################################################### - -function _can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) where T <: FieldElement - @assert base_ring(A) == base_ring(B) - if side === :right - @assert nrows(A) == nrows(B) - return __can_solve_with_kernel(A, B) - elseif side === :left - b, C, K = __can_solve_with_kernel(transpose(A), transpose(B)) - @assert ncols(A) == ncols(B) - if b - return b, transpose(C), transpose(K) - else - return b, C, K - end - else - error("Unsupported argument :$side for side: Must be :left or :right") - end -end - -function __can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}) where T <: FieldElement - R = base_ring(A) - mu = [A B] - rk, mu = rref(mu) - p = find_pivot(mu) - if any(i -> i > ncols(A), p) - return false, B, B - end - sol = zero_matrix(R, ncols(A), ncols(B)) - for i = 1:length(p) - for j = 1:ncols(B) - sol[p[i], j] = mu[i, ncols(A) + j] - end - end - nullity = ncols(A) - length(p) - X = zero(A, ncols(A), nullity) - pivots = zeros(Int, max(nrows(A), ncols(A))) - np = rk - j = k = 1 - for i = 1:rk - while is_zero_entry(mu, i, j) - pivots[np + k] = j - j += 1 - k += 1 - end - pivots[i] = j - j += 1 - end - while k <= nullity - pivots[np + k] = j - j += 1 - k += 1 - end - for i = 1:nullity - for j = 1:rk - X[pivots[j], i] = -mu[j, pivots[np + i]] - end - X[pivots[np + i], i] = one(R) - end - return true, sol, X -end - -@doc raw""" - can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}) where T <: RingElement - -If $Ax = B$ is soluble, returns `true, S, K` where `S` is a particular solution -and $K$ is the kernel. Otherwise returns `false, S, K` where $S$ and $K$ are -undefined, though of the right type for type stability. -Tries to solve $Ax = B$ for $x$ if `side = :right` or $xA = B$ if `side = :left`. -""" -function _can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) where T <: RingElement - @assert base_ring(A) == base_ring(B) - if side === :right - @assert nrows(A) == nrows(B) - b, c, K = __can_solve_with_kernel(transpose(A), B) - return b, transpose(c), K - elseif side === :left - b, C, K = __can_solve_with_kernel(A, transpose(B)) - @assert ncols(A) == ncols(B) - if b - return b, C, transpose(K) - else - return b, C, K - end - else - error("Unsupported argument :$side for side: Must be :left or :right") - end -end - -# Note that _a_ must be supplied transposed and the solution is transposed -function __can_solve_with_kernel(a::MatElem{S}, b::MatElem{S}) where S <: RingElement - H, T = hnf_with_transform(a) - z = zero_matrix(base_ring(a), ncols(b), nrows(a)) - l = min(nrows(a), ncols(a)) - b = deepcopy(b) - for i = 1:ncols(b) - for j = 1:l - k = 1 - while k <= ncols(H) && is_zero_entry(H, j, k) - k += 1 - end - if k > ncols(H) - continue - end - q, r = divrem(b[k, i], H[j, k]) - if !iszero(r) - return false, b, b - end - for h = k:ncols(H) - b[h, i] -= q*H[j, h] - end - z[i, j] = q - end - end - if !iszero(b) - return false, b, b - end - for i = nrows(H):-1:1 - for j = 1:ncols(H) - if !is_zero_entry(H, i, j) - N = zero_matrix(base_ring(a), nrows(a), nrows(H) - i) - for k = 1:nrows(N) - for l = 1:ncols(N) - N[k, l] = T[nrows(T) - l + 1, k] - end - end - return true, z*T, N - end - end - end - N = zero(a, nrows(a), nrows(H)) - for i = 1:min(nrows(N), ncols(N)) - N[i, i] = 1 - end - return true, z*T, N -end - ############################################################################### # # Upper triangular solving diff --git a/test/generic/Matrix-test.jl b/test/generic/Matrix-test.jl index bc2d6cf2ca..3c45347264 100644 --- a/test/generic/Matrix-test.jl +++ b/test/generic/Matrix-test.jl @@ -2982,7 +2982,7 @@ end end end -@testset "Generic.Mat.can_solve_with_kernel" begin +@testset "Generic.Mat.can_solve_with_solution_and_kernel" begin # Exact Ring, Exact field for F in [ZZ, QQ] @@ -3002,7 +3002,7 @@ end B = M*X - flag, Sol, K = AbstractAlgebra._can_solve_with_kernel(M, B) + flag, Sol, K = can_solve_with_solution_and_kernel(M, B, side = :right) @test iszero(M*K) @test flag && M*Sol == B @@ -3010,9 +3010,9 @@ end # fully random B = rand(S, -20:20) - flag, Sol, K = AbstractAlgebra._can_solve_with_kernel(M, B) + flag, Sol, K = can_solve_with_solution_and_kernel(M, B, side = :right) flag2, Sol2 = AbstractAlgebra._can_solve_with_solution(M, B) - + @test flag == flag2 if flag @test M*Sol == B @@ -3030,7 +3030,7 @@ end B = X*M - flag, Sol, K = AbstractAlgebra._can_solve_with_kernel(M, B; side=:left) + flag, Sol, K = can_solve_with_solution_and_kernel(M, B) @test iszero(K*M) @test flag && Sol*M == B From ddf433cfa38db42df7e10b1eb095cd5b56812a6c Mon Sep 17 00:00:00 2001 From: Johannes Schmitt Date: Wed, 28 Feb 2024 14:00:35 +0100 Subject: [PATCH 02/10] Remove outdated generic solving methods * `_can_solve_with_solution` * `_can_solve` * `_solve` * `_solve_left` --- src/Matrix.jl | 127 ------------------------------ src/ModuleHomomorphism.jl | 2 +- src/generic/ModuleHomomorphism.jl | 2 +- test/generic/Matrix-test.jl | 74 ++++++++--------- 4 files changed, 39 insertions(+), 166 deletions(-) diff --git a/src/Matrix.jl b/src/Matrix.jl index bd5d7ed4d5..bcb9abb7cb 100644 --- a/src/Matrix.jl +++ b/src/Matrix.jl @@ -3154,36 +3154,6 @@ function _solve_rational(M::MatElem{T}, b::MatElem{T}) where {T <: PolyRingElem} end end -# @doc raw""" -# solve(a::MatElem{S}, b::MatElem{S}) where {S <: RingElement} -# -# Given an $m\times r$ matrix $a$ over a ring and an $m\times n$ matrix $b$ -# over the same ring, return an $r\times n$ matrix $x$ such that $ax = b$. If -# no such matrix exists, an exception is raised. -# See also [`_solve_left`](@ref). -# """ -function _solve(a::MatElem{S}, b::MatElem{S} - ) where S <: RingElement - can, X = _can_solve_with_solution(a, b, side = :right) - can || throw(ArgumentError("Unable to solve linear system")) - return X -end - -@doc raw""" - _solve_left(a::MatElem{S}, b::MatElem{S}) where S <: RingElement - -Given an $r\times n$ matrix $a$ over a ring and an $m\times n$ matrix $b$ -over the same ring, return an $m\times r$ matrix $x$ such that $xa = b$. If -no such matrix exists, an exception is raised. -See also [`solve`](@ref). -""" -function _solve_left(a::MatElem{S}, b::MatElem{S} - ) where S <: RingElement - (flag, x) = _can_solve_with_solution(a, b; side = :left) - flag || error("Unable to solve linear system") - return x -end - # Find the pivot columns of an rref matrix function find_pivot(A::MatElem{T}) where T <: RingElement p = Int[] @@ -3586,103 +3556,6 @@ function _can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = : error("Unsupported argument :$side for side: Must be :left or :right.") end end - -@doc raw""" - _can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement - -Given two matrices $a$ and $b$ over the same ring, try to solve $ax = b$ -if `side` is `:right` or $xa = b$ if `side` is `:left`. In either case, -return a tuple `(flag, x)`. If a solution exists, `flag` is set to true and -`x` is a solution. If no solution exists, `flag` is set to false and `x` -is arbitrary. If the dimensions of $a$ and $b$ are incompatible, an exception -is raised. -""" -function _can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement - if side == :right - (f, x) = _can_solve_with_solution(transpose(a), transpose(b); side=:left) - return (f, transpose(x)) - elseif side == :left - @assert ncols(a) == ncols(b) - H, T = hnf_with_transform(a) - b = deepcopy(b) - z = zero(a, nrows(b), nrows(a)) - l = min(ncols(a), nrows(a)) - t = base_ring(a)() - for i = 1:nrows(b) - for j = 1:l - k = 1 - while k <= ncols(H) && is_zero_entry(H, j, k) - k += 1 - end - if k > ncols(H) - continue - end - q, r = divrem(b[i, k], H[j, k]) - if r != 0 - return (false, zero(a, 0, 0)) - end - z[i, j] = q - q = -q - for h = k:ncols(H) - t = mul!(t, q, H[j, h]) - b[i, h] = addeq!(b[i, h], t) - end - end - end - if b != 0 - return (false, zero(a, 0, 0)) - end - return (true, z*T) - else - error("Unsupported argument :$side for side: Must be :left or :right.") - end -end - -function _can_solve_with_solution(A::MatElem{T}, B::MatElem{T}; side::Symbol = :right) where T <: FieldElement - if side == :right - (f, x) = _can_solve_with_solution(transpose(A), transpose(B), side = :left) - return (f, transpose(x)) - elseif side == :left - R = base_ring(A) - ncols(A) != ncols(B) && error("Incompatible matrices") - mu = zero_matrix(R, ncols(A), nrows(A) + nrows(B)) - for i = 1:ncols(A) - for j = 1:nrows(A) - mu[i, j] = A[j, i] - end - for j = 1:nrows(B) - mu[i, nrows(A) + j] = B[j, i] - end - end - rk, mu = rref(mu) - p = find_pivot(mu) - if any(i -> i > nrows(A), p) - return (false, zero(A, 0, 0)) - end - sol = zero_matrix(R, nrows(B), nrows(A)) - for i = 1:length(p) - for j = 1:nrows(B) - sol[j, p[i]] = mu[i, nrows(A) + j] - end - end - return (true, sol) - else - error("Unsupported argument :$side for side: Must be :left or :right.") - end -end - -@doc raw""" - _can_solve(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement - -Given two matrices $a$ and $b$ over the same ring, check the solubility -of $ax = b$ if `side` is `:right` or $xa = b$ if `side` is `:left`. -Return true if a solution exists, false otherwise. If the dimensions -of $a$ and $b$ are incompatible, an exception is raised. If a solution -should be computed as well, use `can_solve_with_solution` instead. -""" -function _can_solve(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement - return _can_solve_with_solution(a, b; side=side)[1] -end ############################################################################### # diff --git a/src/ModuleHomomorphism.jl b/src/ModuleHomomorphism.jl index 385d48d12b..00ce4da1cf 100644 --- a/src/ModuleHomomorphism.jl +++ b/src/ModuleHomomorphism.jl @@ -151,7 +151,7 @@ function preimage(f::Map(FPModuleHomomorphism), v::FPModuleElem{T}) where end end # Find left inverse of mat - x = AbstractAlgebra._solve_left(matr, Generic._matrix(v)) + x = solve(matr, Generic._matrix(v)) if q != 0 x = matrix(R, 1, m, T[x[1, i] for i in 1:m]) end diff --git a/src/generic/ModuleHomomorphism.jl b/src/generic/ModuleHomomorphism.jl index b77057a85a..136e53ba8f 100644 --- a/src/generic/ModuleHomomorphism.jl +++ b/src/generic/ModuleHomomorphism.jl @@ -133,7 +133,7 @@ function ModuleIsomorphism(M1::AbstractAlgebra.FPModule{T}, end # Find left inverse of mat N = identity_matrix(R, n) - X = AbstractAlgebra._solve_left(mat, N) + X = solve(mat, N) # Construct matrix of inverse homomorphism from first m columns of X M_inv = X[:, 1:m] end diff --git a/test/generic/Matrix-test.jl b/test/generic/Matrix-test.jl index 3c45347264..fd7853df5d 100644 --- a/test/generic/Matrix-test.jl +++ b/test/generic/Matrix-test.jl @@ -2074,7 +2074,7 @@ end flag1, X, d = Generic._can_solve_with_solution_fflu(A, B) A2 = change_base_ring(QQ, A) B2 = change_base_ring(QQ, B) - flag2, X2 = AbstractAlgebra._can_solve_with_solution(A2, B2) + flag2, X2 = can_solve_with_solution(A2, B2, side = :right) @test flag1 == flag2 end end @@ -2275,7 +2275,7 @@ end M = randmat_with_rank(S, dim, -100:100) b = rand(U, -100:100) - x = AbstractAlgebra._solve(M, b) + x = solve(M, b, side = :right) @test M*x == b end @@ -2319,7 +2319,7 @@ end @test M*x == d*b end -@testset "Generic.Mat.solve" begin +@testset "Generic.Mat.solve_right" begin for R in [ZZ, QQ] for iter = 1:40 for dim = 0:5 @@ -2334,7 +2334,7 @@ end M = rand(S, -20:20) B = M*X1 - X = AbstractAlgebra._solve(M, B) + X = solve(M, B, side = :right) @test M*X == B end @@ -2355,14 +2355,14 @@ end M = rand(S, -1:2, -10:10) B = M*X1 - X = AbstractAlgebra._solve(M, B) + X = solve(M, B, side = :right) @test M*X == B end end end -@testset "Generic.Mat._solve_left" begin +@testset "Generic.Mat.solve_left" begin for R in [ZZ, QQ] for iter = 1:40 for dim = 0:5 @@ -2377,7 +2377,7 @@ end M = rand(U, -20:20) B = X1*M - X = AbstractAlgebra._solve_left(M, X1*M) + X = solve(M, X1*M) @test X*M == B end @@ -2399,7 +2399,7 @@ end M = rand(U, -1:2, -10:10) B = X1*M - X = AbstractAlgebra._solve_left(M, X1*M) + X = solve(M, X1*M) @test X*M == B end @@ -2420,13 +2420,13 @@ end X2 = rand(U, -1:2, -10:10) b = M*X2 - flag, X = AbstractAlgebra._can_solve_with_solution(M, b) + flag, X = can_solve_with_solution(M, b, side = :right) @test flag && M*X == b b = X2*M - flag, X = AbstractAlgebra._can_solve_with_solution(M, b; side=:left) + flag, X = can_solve_with_solution(M, b) @test flag && X*M == b end @@ -2446,7 +2446,7 @@ end M = change_base_ring(S, M) b = change_base_ring(S, b) - flag, X = AbstractAlgebra._can_solve_with_solution(M, b) + flag, X = can_solve_with_solution(M, b, side = :right) @test flag && M*X == b @@ -2457,7 +2457,7 @@ end M = change_base_ring(S, M) b = change_base_ring(S, b) - flag, X = AbstractAlgebra._can_solve_with_solution(M, b; side=:left) + flag, X = can_solve_with_solution(M, b) @test flag && X*M == b end @@ -2477,9 +2477,9 @@ end M = rand(U, -20:20) B = M*X1 - (flag, X) = AbstractAlgebra._can_solve_with_solution(M, M*X1) + (flag, X) = can_solve_with_solution(M, M*X1, side = :right) - @test AbstractAlgebra._can_solve(M, B) + @test can_solve(M, B, side = :right) @test flag && M*X == B end @@ -2491,9 +2491,9 @@ end M = rand(U, -20:20) B = X1*M - (flag, X) = AbstractAlgebra._can_solve_with_solution(M, X1*M; side = :left) + (flag, X) = can_solve_with_solution(M, X1*M) - @test AbstractAlgebra._can_solve(M, B; side = :left) + @test can_solve(M, B) @test flag && X*M == B end end @@ -2516,14 +2516,14 @@ end M = rand(U, -1:2, -10:10) B = M*X1 - (flag, X) = AbstractAlgebra._can_solve_with_solution(M, M*X1) + (flag, X) = can_solve_with_solution(M, M*X1, side = :right) - @test AbstractAlgebra._can_solve(M, B) + @test can_solve(M, B, side = :right) @test flag && M*X == B - (flag, X) = AbstractAlgebra._can_solve_with_solution(M, M*X1; side = :right) + (flag, X) = can_solve_with_solution(M, M*X1, side = :right) - @test AbstractAlgebra._can_solve(M, B; side = :right) + @test can_solve(M, B, side = :right) @test flag && M*X == B end @@ -2535,9 +2535,9 @@ end M = rand(U, -1:2, -10:10) B = X1*M - (flag, X) = AbstractAlgebra._can_solve_with_solution(M, X1*M; side = :left) + (flag, X) = can_solve_with_solution(M, X1*M) - @test AbstractAlgebra._can_solve(M, B; side = :left) + @test can_solve(M, B) @test flag && X*M == B end end @@ -2547,12 +2547,12 @@ end M = matrix(R, 1, 1, [x]) X = matrix(R, 1, 1, [1]) - @assert !AbstractAlgebra._can_solve(M, X) - (flag, _) = AbstractAlgebra._can_solve_with_solution(M, X) + @assert !can_solve(M, X, side = :right) + (flag, _) = can_solve_with_solution(M, X, side = :right) @assert !flag - @assert !AbstractAlgebra._can_solve(M, X; side = :left) - (flag, _) = AbstractAlgebra._can_solve_with_solution(M, X; side = :left) + @assert !can_solve(M, X) + (flag, _) = can_solve_with_solution(M, X) @assert !flag end @@ -2560,12 +2560,12 @@ end M = matrix(ZZ, 2, 2, [1, 1, 1, 1]) X = matrix(ZZ, 2, 1, [1, 0]) - @assert !AbstractAlgebra._can_solve(M, X) - (flag, _) = AbstractAlgebra._can_solve_with_solution(M, X) + @assert !can_solve(M, X, side = :right) + (flag, _) = can_solve_with_solution(M, X, side = :right) @assert !flag - @assert !AbstractAlgebra._can_solve(M, transpose(X); side = :left) - (flag, _) = AbstractAlgebra._can_solve_with_solution(M, transpose(X); side = :left) + @assert !can_solve(M, transpose(X)) + (flag, _) = can_solve_with_solution(M, transpose(X)) @assert !flag end @@ -2573,18 +2573,18 @@ end M = matrix(ZZ, 3, 3, [2, 0, 0, 0, 1, 0, 0, 0, 1]) X1 = matrix(ZZ, 3, 1, [1, 0, 0]) - @assert !AbstractAlgebra._can_solve(M, X1) - (flag, X) = AbstractAlgebra._can_solve_with_solution(M, X1) + @assert !can_solve(M, X1, side = :right) + (flag, X) = can_solve_with_solution(M, X1, side = :right) @assert !flag X2 = matrix(ZZ, 2, 3, [1, 0, 0, 0, 1, 0]) - @assert !AbstractAlgebra._can_solve(M, X2; side = :left) - (flag, _) = AbstractAlgebra._can_solve_with_solution(M, X2; side = :left) + @assert !can_solve(M, X2) + (flag, _) = can_solve_with_solution(M, X2) @assert !flag end - @test_throws Exception AbstractAlgebra._can_solve_with_solution(matrix(ZZ, 2, 2, [1, 0, 0, 1]), matrix(ZZ, 2, 1, [2, 3]), side = :aaa) - @test_throws TypeError AbstractAlgebra._can_solve_with_solution(matrix(ZZ, 2, 2, [1, 0, 0, 1]), matrix(ZZ, 2, 1, [2, 3]), side = "right") + @test_throws Exception can_solve_with_solution(matrix(ZZ, 2, 2, [1, 0, 0, 1]), matrix(ZZ, 2, 1, [2, 3]), side = :aaa) + @test_throws TypeError can_solve_with_solution(matrix(ZZ, 2, 2, [1, 0, 0, 1]), matrix(ZZ, 2, 1, [2, 3]), side = "right") end @testset "Generic.Mat.can_solve_with_solution_interpolation" begin @@ -3011,7 +3011,7 @@ end B = rand(S, -20:20) flag, Sol, K = can_solve_with_solution_and_kernel(M, B, side = :right) - flag2, Sol2 = AbstractAlgebra._can_solve_with_solution(M, B) + flag2, Sol2 = can_solve_with_solution(M, B, side = :right) @test flag == flag2 if flag From e09b860efa97ce5fc158a7cda6908af443bbe2b5 Mon Sep 17 00:00:00 2001 From: Johannes Schmitt Date: Wed, 28 Feb 2024 15:52:16 +0100 Subject: [PATCH 03/10] Specialized linear solving for `FracElem` --- src/Matrix.jl | 62 -------------------------------------- src/Solve.jl | 82 +++++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 80 insertions(+), 64 deletions(-) diff --git a/src/Matrix.jl b/src/Matrix.jl index bcb9abb7cb..93f21f6564 100644 --- a/src/Matrix.jl +++ b/src/Matrix.jl @@ -3495,68 +3495,6 @@ function _can_solve_left_reduced_triu(r::MatElem{T}, return true, x end -function _can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: FracElem{T} where T <: PolyRingElem - if side == :left - (f, x) = _can_solve_with_solution(transpose(a), transpose(b); side=:right) - return (f, transpose(x)) - elseif side == :right - d = numerator(one(base_ring(a))) - for i = 1:nrows(a) - for j = 1:ncols(a) - d = lcm(d, denominator(a[i, j])) - end - end - for i = 1:nrows(b) - for j = 1:ncols(b) - d = lcm(d, denominator(b[i, j])) - end - end - A = matrix(parent(d), nrows(a), ncols(a), [numerator(a[i, j]*d) for i in 1:nrows(a) for j in 1:ncols(a)]) - B = matrix(parent(d), nrows(b), ncols(b), [numerator(b[i, j]*d) for i in 1:nrows(b) for j in 1:ncols(b)]) - flag = false - x = similar(A, ncols(A), nrows(B)) - den = one(base_ring(a)) - try - flag, x, den = _can_solve_with_solution_interpolation(A, B) - catch - flag, x, den = _can_solve_with_solution_fflu(A, B) - end - X = change_base_ring(base_ring(a), x) - X = divexact(X, base_ring(a)(den)) - return flag, X - else - error("Unsupported argument :$side for side: Must be :left or :right.") - end -end - -# The fflu approach is the fastest over a fraction field (see benchmarks on PR 661) -function _can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: Union{FracElem, Rational{BigInt}} - if side == :left - (f, x) = _can_solve_with_solution(transpose(a), transpose(b); side=:right) - return (f, transpose(x)) - elseif side == :right - d = numerator(one(base_ring(a))) - for i = 1:nrows(a) - for j = 1:ncols(a) - d = lcm(d, denominator(a[i, j])) - end - end - for i = 1:nrows(b) - for j = 1:ncols(b) - d = lcm(d, denominator(b[i, j])) - end - end - A = matrix(parent(d), nrows(a), ncols(a), [numerator(a[i, j]*d) for i in 1:nrows(a) for j in 1:ncols(a)]) - B = matrix(parent(d), nrows(b), ncols(b), [numerator(b[i, j]*d) for i in 1:nrows(b) for j in 1:ncols(b)]) - flag, x, den = _can_solve_with_solution_fflu(A, B) - X = change_base_ring(base_ring(a), x) - X = divexact(X, base_ring(a)(den)) - return flag, X - else - error("Unsupported argument :$side for side: Must be :left or :right.") - end -end - ############################################################################### # # Inverse diff --git a/src/Solve.jl b/src/Solve.jl index ed307059d6..3d3a5166a4 100644 --- a/src/Solve.jl +++ b/src/Solve.jl @@ -2,7 +2,8 @@ module Solve using AbstractAlgebra -import AbstractAlgebra: base_ring, nrows, ncols, matrix, rank, Generic, kernel +import AbstractAlgebra: base_ring, nrows, ncols, matrix, rank, Generic, kernel, + _can_solve_with_solution_fflu, _can_solve_with_solution_interpolation ################################################################################ # @@ -353,7 +354,7 @@ end ################################################################################ # -# Internal functionality +# Generic internal functionality # ################################################################################ @@ -535,6 +536,83 @@ function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbo return true, sol, kernel(C, side = side) end +################################################################################ +# +# Special internal functonality +# +################################################################################ + +# _can_solve_internal_no_check methods which only work or are efficient for +# certain types + +function _common_denominator(A::MatElem{T}) where T <: Union{FracElem, Rational{BigInt}} + d = numerator(one(base_ring(A))) + @inbounds for j in 1:ncols(A) + for i in 1:nrows(A) + d = lcm(d, denominator(A[i, j])) + end + end + return d +end + +# The fflu approach is the fastest over a fraction field (see benchmarks on PR 661) +function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: Union{FracElem, Rational{BigInt}} + + if side === :left + fl, _sol, _K = _can_solve_internal_no_check(lazy_transpose(A), lazy_transpose(b), task, side = :right) + # This does not return LazyTransposedMat for sol because the matrices are made integral + return fl, transpose(_sol), data(_K) + end + + d = lcm(_common_denominator(A), _common_denominator(b)) + + Aint = matrix(parent(d), nrows(A), ncols(A), [numerator(A[i, j]*d) for i in 1:nrows(A) for j in 1:ncols(A)]) + bint = matrix(parent(d), nrows(b), ncols(b), [numerator(b[i, j]*d) for i in 1:nrows(b) for j in 1:ncols(b)]) + flag, sol_int, den = _can_solve_with_solution_fflu(Aint, bint) + sol = change_base_ring(base_ring(A), sol_int) + sol = divexact(sol, base_ring(A)(den)) + + if task === :with_kernel + # I don't know how to compute the kernel using an (ff)lu factoring + return flag, sol, kernel(A, side = :right) + else + return flag, sol, zero(A, 0, 0) + end +end + +function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FracElem{TT} where TT <: PolyRingElem + + if side === :left + fl, _sol, _K = _can_solve_internal_no_check(lazy_transpose(A), lazy_transpose(b), task, side = :right) + # This does not return a LazyTransposedMat for sol because the matrices are made integral + return fl, transpose(_sol), data(_K) + end + + d = lcm(_common_denominator(A), _common_denominator(b)) + + Aint = matrix(parent(d), nrows(A), ncols(A), [numerator(A[i, j]*d) for i in 1:nrows(A) for j in 1:ncols(A)]) + bint = matrix(parent(d), nrows(b), ncols(b), [numerator(b[i, j]*d) for i in 1:nrows(b) for j in 1:ncols(b)]) + + flag = false + sol_int = similar(Aint, ncols(Aint), nrows(bint)) + den = one(base_ring(A)) + try + flag, sol_int, den = _can_solve_with_solution_interpolation(Aint, bint) + catch + flag, sol_int, den = _can_solve_with_solution_fflu(Aint, bint) + end + + sol = change_base_ring(base_ring(A), sol_int) + sol = divexact(sol, base_ring(A)(den)) + + if task === :with_kernel + # I don't know how to compute the kernel using an (ff)lu factoring + return flag, sol, kernel(A, side = :right) + else + return flag, sol, zero(A, 0, 0) + end +end + ################################################################################ # # Internals for solving of row reduced matrices From 84601110d58d8849baf46650ca35bd4ba8638b87 Mon Sep 17 00:00:00 2001 From: Johannes Schmitt Date: Fri, 1 Mar 2024 14:51:43 +0100 Subject: [PATCH 04/10] Use LU factoring in solve context over fields --- src/Solve.jl | 75 +++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 57 insertions(+), 18 deletions(-) diff --git a/src/Solve.jl b/src/Solve.jl index 3d3a5166a4..56d4df79e2 100644 --- a/src/Solve.jl +++ b/src/Solve.jl @@ -61,6 +61,9 @@ mutable struct SolveCtx{T, MatT, TranspMatT} pivots::Vector{Int} # pivot and non-pivot columns of red pivots_transp::Vector{Int} # pivot and non-pivot columns of red_transp + kernel_left::MatT + kernel_right::MatT + function SolveCtx{T, MatT, TranspMatT}(A::MatT) where {T, MatT <: MatElem{T}, TranspMatT <: MatElem{T}} z = new{T, MatT, TranspMatT}() z.A = A @@ -86,14 +89,16 @@ end matrix(C::SolveCtx) = C.A function _init_reduce(C::SolveCtx{<:FieldElement}) - if isdefined(C, :red) && isdefined(C, :trafo) + if isdefined(C, :red) && isdefined(C, :lu_perm) return nothing end - r, R, U = _rref_with_transformation(matrix(C)) + B = deepcopy(matrix(C)) + p = one(SymmetricGroup(nrows(matrix(C)))) + r = lu!(p, B) set_rank!(C, r) - C.red = R - C.trafo = U + C.red = B + C.lu_perm = p return nothing end @@ -109,14 +114,16 @@ function _init_reduce(C::SolveCtx{<:RingElement}) end function _init_reduce_transpose(C::SolveCtx{<:FieldElement}) - if isdefined(C, :red_transp) && isdefined(C, :trafo_transp) + if isdefined(C, :red_transp) && isdefined(C, :lu_perm_transp) return nothing end - r, R, U = _rref_with_transformation(lazy_transpose(matrix(C))) + B = lazy_transpose(deepcopy(matrix(C))) + p = one(SymmetricGroup(ncols(matrix(C)))) + r = lu!(p, B) set_rank!(C, r) - C.red_transp = R - C.trafo_transp = U + C.red_transp = B + C.lu_perm_transp = p return nothing end @@ -331,12 +338,18 @@ end function kernel(C::SolveCtx{<:FieldElement}; side::Symbol = :left) check_option(side, [:right, :left], "side") + # I don't know how to compute the kernel using a LU factoring, so we call the + # "usual" method if side === :right - return _kernel_of_rref(reduced_matrix(C), rank(C), pivot_and_non_pivot_cols(C))[2] + if !isdefined(C, :kernel_right) + C.kernel_right = kernel(matrix(C), side = :right) + end + return C.kernel_right else - nullity, X = _kernel_of_rref(reduced_matrix_of_transpose(C), rank(C), pivot_and_non_pivot_cols_of_transpose(C)) - # X is of type LazyTransposeMatElem - return data(X) + if !isdefined(C, :kernel_left) + C.kernel_left = kernel(matrix(C), side = :left) + end + return C.kernel_left end end @@ -344,11 +357,16 @@ function kernel(C::SolveCtx{<:RingElement}; side::Symbol = :left) check_option(side, [:right, :left], "side") if side === :right - return _kernel_of_hnf(matrix(C), reduced_matrix_of_transpose(C), transformation_matrix_of_transpose(C))[2] + if !isdefined(C, :kernel_right) + C.kernel_right = _kernel_of_hnf(matrix(C), reduced_matrix_of_transpose(C), transformation_matrix_of_transpose(C))[2] + end + return C.kernel_right else - nullity, X = _kernel_of_hnf(lazy_transpose(matrix(C)), reduced_matrix(C), transformation_matrix(C)) - # X is of type LazyTransposeMatElem - return data(X) + if !isdefined(C, :kernel_left) + nullity, X = _kernel_of_hnf(lazy_transpose(matrix(C)), reduced_matrix(C), transformation_matrix(C)) + C.kernel_left = data(X) + end + return C.kernel_left end end @@ -509,9 +527,9 @@ end # _can_solve_internal_no_check over FIELDS with SOLVE CONTEXT function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FieldElement if side === :right - fl, sol = _can_solve_with_rref(b, transformation_matrix(C), rank(C), pivot_and_non_pivot_cols(C), task) + fl, sol = _can_solve_with_lu(matrix(C), b, reduced_matrix(C), lu_permutation(C), rank(C)) else - fl, _sol = _can_solve_with_rref(lazy_transpose(b), transformation_matrix_of_transpose(C), rank(C), pivot_and_non_pivot_cols_of_transpose(C), task) + fl, _sol = _can_solve_with_lu(lazy_transpose(matrix(C)), lazy_transpose(b), reduced_matrix_of_transpose(C), lu_permutation_of_transpose(C), rank(C)) sol = data(_sol) end if !fl || task !== :with_kernel @@ -643,9 +661,29 @@ function _can_solve_with_rref(b::MatElem{T}, U::MatElem{T}, r::Int, pivots::Vect return true, sol end +# Solve Ax = b with LU and p a LU decomposition of A of rank r. +# Takes same options for `task` as _can_solve_internal but only returns (flag, solution) +# and no kernel. +function _can_solve_with_lu(A::MatElem{T}, b::MatElem{T}, LU::MatElem{T}, p::Generic.Perm, r::Int) where T <: FieldElement + y = AbstractAlgebra._solve_lu_precomp(p, LU, b) + + # _solve_lu_precomp only takes care of the first r rows + # We now need to check whether the remaining rows are fine as well + n = nrows(A) + fl = true + if r < n + b2 = p*b + A2 = p*A + A3 = A2[r + 1:n, :] + fl = A3*y == b2[r + 1:n, :] + end + return fl, y +end + # Compute a matrix N with RN == 0 where the columns of N give a basis for the kernel. # R must be in rref of rank r and pivots must be of length ncols(R) with the pivot # columns in the first r entries and the non-pivot columns in the remaining entries. +# Only used by Nemo function _kernel_of_rref(R::MatElem{T}, r::Int, pivots::Vector{Int}) where T <: FieldElement @assert length(pivots) == ncols(R) nullity = ncols(R) - r @@ -713,6 +751,7 @@ function _kernel_of_hnf(A::MatElem{T}, H::MatElem{T}, U::MatElem{T}) where T <: end # Copied from Hecke, to be replaced with echelon_form_with_transformation eventually +# Only used by Nemo function _rref_with_transformation(M::MatElem{T}) where T <: FieldElement n = hcat(deepcopy(M), identity_matrix(base_ring(M), nrows(M))) rref!(n) From 83c291254c42d1d27ca6231e6500a3decb016537 Mon Sep 17 00:00:00 2001 From: Johannes Schmitt Date: Fri, 1 Mar 2024 16:18:22 +0100 Subject: [PATCH 05/10] `show` for `SolveCtx` --- src/Solve.jl | 29 +++++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/src/Solve.jl b/src/Solve.jl index 56d4df79e2..eec8b4ab1f 100644 --- a/src/Solve.jl +++ b/src/Solve.jl @@ -2,8 +2,19 @@ module Solve using AbstractAlgebra -import AbstractAlgebra: base_ring, nrows, ncols, matrix, rank, Generic, kernel, - _can_solve_with_solution_fflu, _can_solve_with_solution_interpolation +import AbstractAlgebra: + base_ring, + _can_solve_with_solution_fflu, + _can_solve_with_solution_interpolation, + Generic, + kernel, + matrix, + nrows, + ncols, + PrettyPrinting, + rank + +import Base: show ################################################################################ # @@ -76,6 +87,20 @@ mutable struct SolveCtx{T, MatT, TranspMatT} end end +function Base.show(io::IO, ::MIME"text/plain", C::SolveCtx) + io = PrettyPrinting.pretty(io) + println(io, "Linear solving context object of matrix") + print(io, PrettyPrinting.Indent()) + show(io, MIME"text/plain"(), matrix(C)) + print(io, PrettyPrinting.Dedent()) +end + +function Base.show(io::IO, C::SolveCtx) + PrettyPrinting.@show_name(io, C) + + print(io, "Linear solving context object") +end + @doc raw""" solve_init(A::MatElem) From 7184bfc6e8b418d384cb1da5ee01c43dedf992d8 Mon Sep 17 00:00:00 2001 From: Johannes Schmitt Date: Wed, 6 Mar 2024 14:58:09 +0100 Subject: [PATCH 06/10] Don't break Nemo --- src/Solve.jl | 53 ++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 49 insertions(+), 4 deletions(-) diff --git a/src/Solve.jl b/src/Solve.jl index eec8b4ab1f..2a8f5516bb 100644 --- a/src/Solve.jl +++ b/src/Solve.jl @@ -114,6 +114,29 @@ end matrix(C::SolveCtx) = C.A function _init_reduce(C::SolveCtx{<:FieldElement}) + if isdefined(C, :red) && isdefined(C, :trafo) + return nothing + end + + r, R, U = _rref_with_transformation(matrix(C)) + set_rank!(C, r) + C.red = R + C.trafo = U + return nothing +end + +# All matrices over fields in AbstractAlgebra use LU factoring, but Nemo needs +# _init_reduce (and the calling functions) to do an rref. +# Once we make a breaking release of AbstractAlgebra, we can rename +# * _init_reduce_lu -> _init_reduce +# * _init_reduce_transpose_lu -> _init_reduce_transpose +# * reduced_matrix_lu -> reduced_matrix +# * reduced_matrix_of_transpose_lu -> reduced_matrix_of_transpose +# and then Nemo can have its own version of these functions. +# (transformation_matrix and transformation_matrix_of_transpose should probably +# throw an error for matrices over fields then.) + +function _init_reduce_lu(C::SolveCtx{<:FieldElement}) if isdefined(C, :red) && isdefined(C, :lu_perm) return nothing end @@ -139,6 +162,18 @@ function _init_reduce(C::SolveCtx{<:RingElement}) end function _init_reduce_transpose(C::SolveCtx{<:FieldElement}) + if isdefined(C, :red_transp) && isdefined(C, :trafo_transp) + return nothing + end + + r, R, U = _rref_with_transformation(lazy_transpose(matrix(C))) + set_rank!(C, r) + C.red_transp = R + C.trafo_transp = U + return nothing +end + +function _init_reduce_transpose_lu(C::SolveCtx{<:FieldElement}) if isdefined(C, :red_transp) && isdefined(C, :lu_perm_transp) return nothing end @@ -173,13 +208,23 @@ function reduced_matrix_of_transpose(C::SolveCtx) return C.red_transp end +function reduced_matrix_lu(C::SolveCtx) + _init_reduce_lu(C) + return C.red +end + +function reduced_matrix_of_transpose_lu(C::SolveCtx) + _init_reduce_transpose_lu(C) + return C.red_transp +end + function lu_permutation(C::SolveCtx) - _init_reduce(C) + _init_reduce_lu(C) return C.lu_perm end function lu_permutation_of_transpose(C::SolveCtx) - _init_reduce_transpose(C) + _init_reduce_transpose_lu(C) return C.lu_perm_transp end @@ -552,9 +597,9 @@ end # _can_solve_internal_no_check over FIELDS with SOLVE CONTEXT function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FieldElement if side === :right - fl, sol = _can_solve_with_lu(matrix(C), b, reduced_matrix(C), lu_permutation(C), rank(C)) + fl, sol = _can_solve_with_lu(matrix(C), b, reduced_matrix_lu(C), lu_permutation(C), rank(C)) else - fl, _sol = _can_solve_with_lu(lazy_transpose(matrix(C)), lazy_transpose(b), reduced_matrix_of_transpose(C), lu_permutation_of_transpose(C), rank(C)) + fl, _sol = _can_solve_with_lu(lazy_transpose(matrix(C)), lazy_transpose(b), reduced_matrix_of_transpose_lu(C), lu_permutation_of_transpose(C), rank(C)) sol = data(_sol) end if !fl || task !== :with_kernel From ba681eadbdd0b7e6f7e594531f18327b6ebd9a1e Mon Sep 17 00:00:00 2001 From: Johannes Schmitt Date: Thu, 14 Mar 2024 16:04:17 +0100 Subject: [PATCH 07/10] Use 'FFLU' in `SolveCtx` for fraction fields --- src/Matrix.jl | 45 +++++++++++ src/Solve.jl | 201 ++++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 241 insertions(+), 5 deletions(-) diff --git a/src/Matrix.jl b/src/Matrix.jl index 93f21f6564..5091d9c6a7 100644 --- a/src/Matrix.jl +++ b/src/Matrix.jl @@ -1578,6 +1578,51 @@ function *(P::Perm, x::MatrixElem{T}) where T <: NCRingElement return z end +@doc raw""" + *(x::MatrixElem{T}, P::Perm) where T <: NCRingElement + +Apply the pemutation $P$ to the columns of the matrix $x$ and return the result. + +# Examples + +```jldoctest; setup = :(using AbstractAlgebra) +julia> R, t = polynomial_ring(QQ, "t") +(Univariate polynomial ring in t over rationals, t) + +julia> S = matrix_space(R, 3, 3) +Matrix space of 3 rows and 3 columns + over univariate polynomial ring in t over rationals + +julia> G = SymmetricGroup(3) +Full symmetric group over 3 elements + +julia> A = S([t + 1 t R(1); t^2 t t; R(-2) t + 2 t^2 + t + 1]) +[t + 1 t 1] +[ t^2 t t] +[ -2 t + 2 t^2 + t + 1] + +julia> P = G([1, 3, 2]) +(2,3) + +julia> B = A*P +[t + 1 1 t] +[ t^2 t t] +[ -2 t^2 + t + 1 t + 2] + +``` +""" +function *(x::MatrixElem{T}, P::Perm) where T <: NCRingElement + z = similar(x) + m = nrows(x) + n = ncols(x) + for i = 1:m + for j = 1:n + z[i, P[j]] = x[i, j] + end + end + return z +end + ############################################################################### # # LU factorisation diff --git a/src/Solve.jl b/src/Solve.jl index 2a8f5516bb..7f768ce52a 100644 --- a/src/Solve.jl +++ b/src/Solve.jl @@ -3,16 +3,22 @@ module Solve using AbstractAlgebra import AbstractAlgebra: - base_ring, + @attributes, _can_solve_with_solution_fflu, _can_solve_with_solution_interpolation, + _solve_fflu_precomp, + base_ring, + fflu!, Generic, + get_attribute, + has_attribute, kernel, matrix, nrows, ncols, PrettyPrinting, - rank + rank, + set_attribute! import Base: show @@ -58,15 +64,18 @@ Base.similar(M::LazyTransposeMatElem, i::Int, j::Int) = lazy_transpose(similar(d # ################################################################################ -mutable struct SolveCtx{T, MatT, TranspMatT} +# attribute storage is enable so that the fraction field types can store integral +# matrices. These don't go in `red` or `red_transp` because they don't have the +# correct type. Changing the type parameters would be a breaking change. +@attributes mutable struct SolveCtx{T, MatT, TranspMatT} A::MatT # matrix giving the linear system red::MatT # reduced/canonical form of A (rref, hnf, lu) red_transp::TranspMatT # reduced/canonical form of transpose(A) trafo::MatT # transformation: trafo*A == red (not used for lu) trafo_transp::TranspMatT # transformation: trafo_transp*transpose(A) == red_transp # (not used for lu) - lu_perm::Generic.Perm # permutation used for the lu factorization of A - lu_perm_transp::Generic.Perm # permutation used for the lu factorization of transpose(A) + lu_perm::Generic.Perm{Int} # permutation used for the lu factorization of A + lu_perm_transp::Generic.Perm{Int} # permutation used for the lu factorization of transpose(A) rank::Int # rank of A pivots::Vector{Int} # pivot and non-pivot columns of red @@ -150,6 +159,31 @@ function _init_reduce_lu(C::SolveCtx{<:FieldElement}) return nothing end +function _init_reduce_fflu(C::SolveCtx{<:Union{FracElem, Rational{BigInt}}}) + if has_attribute(C, :reduced_matrix_lu) + return nothing + end + + A = matrix(C) + dA = _common_denominator(A) + Aint = matrix(parent(dA), nrows(A), ncols(A), [numerator(A[i, j]*dA) for i in 1:nrows(A) for j in 1:ncols(A)]) + p = one(SymmetricGroup(nrows(A))) + r, dLU = fflu!(p, Aint) + + set_rank!(C, r) + C.lu_perm = p + d = divexact(dA, base_ring(C)(dLU)) + set_attribute!(C, :reduced_matrix_lu => Aint, :scaling_factor_fflu => d) + if r < nrows(A) + A2 = p*A + A3 = A2[r + 1:nrows(A), :] + set_attribute!(C, :permuted_matrix_fflu => A3) + else + set_attribute!(C, :permuted_matrix_fflu => zero(A, 0, ncols(A))) + end + return nothing +end + function _init_reduce(C::SolveCtx{<:RingElement}) if isdefined(C, :red) && isdefined(C, :trafo) return nothing @@ -187,6 +221,32 @@ function _init_reduce_transpose_lu(C::SolveCtx{<:FieldElement}) return nothing end +function _init_reduce_transpose_fflu(C::SolveCtx{<:Union{FracElem, Rational{BigInt}}}) + if has_attribute(C, :reduced_matrix_of_transpose_lu) + return nothing + end + + A = matrix(C) + dA = _common_denominator(A) + # We transpose A at this point! + Aint = matrix(parent(dA), ncols(A), nrows(A), [numerator(A[i, j]*dA) for j in 1:ncols(A) for i in 1:nrows(A)]) + p = one(SymmetricGroup(nrows(Aint))) + r, dLU = fflu!(p, Aint) + + set_rank!(C, r) + C.lu_perm_transp = p + d = divexact(dA, base_ring(C)(dLU)) + set_attribute!(C, :reduced_matrix_of_transpose_lu => Aint, :scaling_factor_of_transpose_fflu => d) + if r < ncols(A) + A2 = A*p + A3 = A2[:, r + 1:ncols(A)] + set_attribute!(C, :permuted_matrix_of_transpose_fflu => A3) + else + set_attribute!(C, :permuted_matrix_of_transpose_fflu => zero(A, nrows(A), 0)) + end + return nothing +end + function _init_reduce_transpose(C::SolveCtx{<:RingElement}) if isdefined(C, :red_transp) && isdefined(C, :trafo_transp) return nothing @@ -213,21 +273,78 @@ function reduced_matrix_lu(C::SolveCtx) return C.red end +function reduced_matrix_lu(C::SolveCtx{<:FracElem{T}}) where T + _init_reduce_fflu(C) + return get_attribute(C, :reduced_matrix_lu)::dense_matrix_type(T) +end + +function reduced_matrix_lu(C::SolveCtx{<:Rational{BigInt}}) + _init_reduce_fflu(C) + return get_attribute(C, :reduced_matrix_lu)::dense_matrix_type(BigInt) +end + function reduced_matrix_of_transpose_lu(C::SolveCtx) _init_reduce_transpose_lu(C) return C.red_transp end +function reduced_matrix_of_transpose_lu(C::SolveCtx{<:FracElem{T}}) where T + _init_reduce_transpose_fflu(C) + return get_attribute(C, :reduced_matrix_of_transpose_lu)::dense_matrix_type(T) +end + +function reduced_matrix_of_transpose_lu(C::SolveCtx{<:Rational{BigInt}}) + _init_reduce_transpose_fflu(C) + return get_attribute(C, :reduced_matrix_of_transpose_lu)::dense_matrix_type(BigInt) +end + +# Factor by which any solution needs to be multiplied. +# This is the chosen denominator of matrix(C) divided by the denominator returned +# by fflu!(matrix(C)). +function scaling_factor_fflu(C::SolveCtx{T}) where {T <: Union{FracElem, Rational{BigInt}}} + _init_reduce_fflu(C) + return get_attribute(C, :scaling_factor_fflu)::T +end + +function scaling_factor_of_transpose_fflu(C::SolveCtx{T}) where {T <: Union{FracElem, Rational{BigInt}}} + _init_reduce_transpose_fflu(C) + return get_attribute(C, :scaling_factor_of_transpose_fflu)::T +end + +# Let A = matrix(C). +# Return the matrix (p*A)[rank(A) + 1:nrows(A), :] where p is lu_permutation(C). +function permuted_matrix_fflu(C::SolveCtx{FldT, MatT}) where {FldT <: Union{FracElem, Rational{BigInt}}, MatT} + _init_reduce_fflu(C) + return get_attribute(C, :permuted_matrix_fflu)::MatT +end + +# Let A = matrix(C). +# Return the matrix (A*p)[:, rank(A) + 1:ncols(A)] where p is lu_permutation_of_transpose(C). +function permuted_matrix_of_transpose_fflu(C::SolveCtx{FldT, MatT}) where {FldT <: Union{FracElem, Rational{BigInt}}, MatT} + _init_reduce_transpose_fflu(C) + return get_attribute(C, :permuted_matrix_of_transpose_fflu)::MatT +end + function lu_permutation(C::SolveCtx) _init_reduce_lu(C) return C.lu_perm end +function lu_permutation(C::SolveCtx{<:Union{FracElem, Rational{BigInt}}}) + _init_reduce_fflu(C) + return C.lu_perm +end + function lu_permutation_of_transpose(C::SolveCtx) _init_reduce_transpose_lu(C) return C.lu_perm_transp end +function lu_permutation_of_transpose(C::SolveCtx{<:Union{FracElem, Rational{BigInt}}}) + _init_reduce_transpose_fflu(C) + return C.lu_perm_transp +end + function transformation_matrix(C::SolveCtx) _init_reduce(C) return C.trafo @@ -248,11 +365,20 @@ end function AbstractAlgebra.rank(C::SolveCtx{<:FieldElement}) if C.rank < 0 + # We should be calling _init_reduce_lu here, but can't because it has to stay + # compatible with Nemo. _init_reduce(C) end return C.rank end +function AbstractAlgebra.rank(C::SolveCtx{<:Union{FracElem, Rational{BigInt}}}) + if C.rank < 0 + _init_reduce_fflu(C) + end + return C.rank +end + AbstractAlgebra.nrows(C::SolveCtx) = nrows(matrix(C)) AbstractAlgebra.ncols(C::SolveCtx) = ncols(matrix(C)) AbstractAlgebra.base_ring(C::SolveCtx) = base_ring(matrix(C)) @@ -668,6 +794,71 @@ function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol end end +function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: Union{FracElem, Rational{BigInt}} + # Split up in separate functions to make the compiler happy + if side === :right + return _can_solve_internal_no_check_right(C, b, task) + else + return _can_solve_internal_no_check_left(C, b, task) + end +end + +function _can_solve_internal_no_check_right(C::SolveCtx{T}, b::MatElem{T}, task::Symbol) where T <: Union{FracElem, Rational{BigInt}} + K = base_ring(C) + db = _common_denominator(b) + bint = matrix(parent(db), nrows(b), ncols(b), [numerator(b[i, j]*db) for i in 1:nrows(b) for j in 1:ncols(b)]) + fl, y_int = _solve_fflu_precomp(lu_permutation(C), reduced_matrix_lu(C), bint) + if !fl + return fl, zero(b, 0, 0), zero(b, 0, 0) + end + # We have fl == true, but we still have to check whether this really is a solution + d = scaling_factor_fflu(C) + y = change_base_ring(K, y_int) + y = y*divexact(d, K(db)) + if rank(C) < nrows(C) + # We have to check whether y is also a solution for the "lower part" + # of the system + pb = lu_permutation(C)*b + pA = permuted_matrix_fflu(C) + fl = pA*y == pb[rank(C) + 1:nrows(C), :] + end + if task === :with_kernel + # I don't know how to compute the kernel using an (ff)lu factoring + return fl, y, kernel(C, side = :right) + else + return fl, y, zero(b, 0, 0) + end +end + +function _can_solve_internal_no_check_left(C::SolveCtx{T}, b::MatElem{T}, task::Symbol) where T <: Union{FracElem, Rational{BigInt}} + K = base_ring(C) + db = _common_denominator(b) + # bint == transpose(b)*db + bint = matrix(parent(db), ncols(b), nrows(b), [numerator(b[i, j]*db) for j in 1:ncols(b) for i in 1:nrows(b)]) + fl, y_int = _solve_fflu_precomp(lu_permutation_of_transpose(C), reduced_matrix_of_transpose_lu(C), bint) + if !fl + return fl, zero(b, 0, 0), zero(b, 0, 0) + end + # We have fl == true, but we still have to check whether this really is a solution + d = scaling_factor_of_transpose_fflu(C) + # transpose y_int + y = matrix(K, ncols(y_int), nrows(y_int), [ K(y_int[i, j]) for j in 1:ncols(y_int) for i in 1:nrows(y_int)]) + y = y*divexact(d, K(db)) + if rank(C) < ncols(C) + # We have to check whether y is also a solution for the "right hand part" + # of the system + pb = b*lu_permutation_of_transpose(C) + pA = permuted_matrix_of_transpose_fflu(C) + fl = y*pA == pb[:, rank(C) + 1:ncols(C)] + end + if task === :with_kernel + # I don't know how to compute the kernel using an (ff)lu factoring + return fl, y, kernel(C, side = :left) + else + return fl, y, zero(b, 0, 0) + end +end + function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FracElem{TT} where TT <: PolyRingElem if side === :left From 3eb60456d757e4e9e0ecd8628c1ca915e721dca5 Mon Sep 17 00:00:00 2001 From: Johannes Schmitt Date: Thu, 14 Mar 2024 16:40:09 +0100 Subject: [PATCH 08/10] Another workaround to not break Nemo --- src/Solve.jl | 34 +++++++++++++++++++++++----------- 1 file changed, 23 insertions(+), 11 deletions(-) diff --git a/src/Solve.jl b/src/Solve.jl index 7f768ce52a..746b413ac3 100644 --- a/src/Solve.jl +++ b/src/Solve.jl @@ -653,14 +653,26 @@ function _can_solve_internal(A::Union{MatElem{T}, SolveCtx{T}}, b::MatElem{T}, t return _can_solve_internal_no_check(A, b, task, side = side) end +# The internal AbstractAlgebra functions are called _can_solve_internal_no_check_dr +# (dr = done right), to allow other packages to overwrite +# _can_solve_internal_no_check without creating any ambiguity. +# Example: AbstractAlgebra defines +# _can_solve_internal_no_check(::SolveCtx{T}, ::MatElem{T}, ...) where {T <: FracElem} +# and Nemo defines +# _can_solve_internal_no_check(::SolveCtx{T, MatT}, ::MatT, ...) where {MatT <: _FieldMatTypes} +# This is ambiguous for T == QQFieldElem and MatT == QQMatrix. +# Once we do a breaking release, we can remove any _dr (and resolve method ambiguities +# in our own packages). +_can_solve_internal_no_check(args...; kwargs...) = _can_solve_internal_no_check_dr(args...; kwargs...) + # _can_solve_internal_no_check over FIELDS -function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FieldElement +function _can_solve_internal_no_check_dr(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FieldElement R = base_ring(A) if side === :left # For side == :left, we pretend that A and b are transposed - fl, _sol, _K = _can_solve_internal_no_check(lazy_transpose(A), lazy_transpose(b), task, side = :right) + fl, _sol, _K = _can_solve_internal_no_check_dr(lazy_transpose(A), lazy_transpose(b), task, side = :right) return fl, data(_sol), data(_K) end @@ -700,13 +712,13 @@ function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol end # _can_solve_internal_no_check over RINGS -function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: RingElement +function _can_solve_internal_no_check_dr(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: RingElement R = base_ring(A) if side === :left # For side == :left, we pretend that A and b are transposed - fl, _sol, _K = _can_solve_internal_no_check(lazy_transpose(A), lazy_transpose(b), task, side = :right) + fl, _sol, _K = _can_solve_internal_no_check_dr(lazy_transpose(A), lazy_transpose(b), task, side = :right) return fl, data(_sol), data(_K) end @@ -721,7 +733,7 @@ function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol end # _can_solve_internal_no_check over FIELDS with SOLVE CONTEXT -function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FieldElement +function _can_solve_internal_no_check_dr(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FieldElement if side === :right fl, sol = _can_solve_with_lu(matrix(C), b, reduced_matrix_lu(C), lu_permutation(C), rank(C)) else @@ -736,7 +748,7 @@ function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbo end # _can_solve_internal_no_check over RINGS with SOLVE CONTEXT -function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: RingElement +function _can_solve_internal_no_check_dr(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: RingElement if side === :right fl, sol = _can_solve_with_hnf(b, reduced_matrix_of_transpose(C), transformation_matrix_of_transpose(C), task) else @@ -770,10 +782,10 @@ function _common_denominator(A::MatElem{T}) where T <: Union{FracElem, Rational{ end # The fflu approach is the fastest over a fraction field (see benchmarks on PR 661) -function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: Union{FracElem, Rational{BigInt}} +function _can_solve_internal_no_check_dr(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: Union{FracElem, Rational{BigInt}} if side === :left - fl, _sol, _K = _can_solve_internal_no_check(lazy_transpose(A), lazy_transpose(b), task, side = :right) + fl, _sol, _K = _can_solve_internal_no_check_dr(lazy_transpose(A), lazy_transpose(b), task, side = :right) # This does not return LazyTransposedMat for sol because the matrices are made integral return fl, transpose(_sol), data(_K) end @@ -794,7 +806,7 @@ function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol end end -function _can_solve_internal_no_check(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: Union{FracElem, Rational{BigInt}} +function _can_solve_internal_no_check_dr(C::SolveCtx{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: Union{FracElem, Rational{BigInt}} # Split up in separate functions to make the compiler happy if side === :right return _can_solve_internal_no_check_right(C, b, task) @@ -859,10 +871,10 @@ function _can_solve_internal_no_check_left(C::SolveCtx{T}, b::MatElem{T}, task:: end end -function _can_solve_internal_no_check(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FracElem{TT} where TT <: PolyRingElem +function _can_solve_internal_no_check_dr(A::MatElem{T}, b::MatElem{T}, task::Symbol; side::Symbol = :left) where T <: FracElem{TT} where TT <: PolyRingElem if side === :left - fl, _sol, _K = _can_solve_internal_no_check(lazy_transpose(A), lazy_transpose(b), task, side = :right) + fl, _sol, _K = _can_solve_internal_no_check_dr(lazy_transpose(A), lazy_transpose(b), task, side = :right) # This does not return a LazyTransposedMat for sol because the matrices are made integral return fl, transpose(_sol), data(_K) end From 5c6fcceee3ad8f16bccd39910cde7c7ac5a12e32 Mon Sep 17 00:00:00 2001 From: Johannes Schmitt Date: Fri, 15 Mar 2024 16:30:30 +0100 Subject: [PATCH 09/10] Increase test coverage (?) --- src/Solve.jl | 4 ++-- test/Solve-test.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Solve.jl b/src/Solve.jl index 746b413ac3..10fdbccf3b 100644 --- a/src/Solve.jl +++ b/src/Solve.jl @@ -98,7 +98,7 @@ end function Base.show(io::IO, ::MIME"text/plain", C::SolveCtx) io = PrettyPrinting.pretty(io) - println(io, "Linear solving context object of matrix") + println(io, "Linear solving context of matrix") print(io, PrettyPrinting.Indent()) show(io, MIME"text/plain"(), matrix(C)) print(io, PrettyPrinting.Dedent()) @@ -107,7 +107,7 @@ end function Base.show(io::IO, C::SolveCtx) PrettyPrinting.@show_name(io, C) - print(io, "Linear solving context object") + print(io, "Linear solving context") end @doc raw""" diff --git a/test/Solve-test.jl b/test/Solve-test.jl index 7279719293..2eeb2ca62a 100644 --- a/test/Solve-test.jl +++ b/test/Solve-test.jl @@ -1,5 +1,5 @@ @testset "Linear solving" begin - for R in [ QQ, ZZ ] + for R in [ QQ, ZZ, GF(101) ] M = matrix(R, [1 2 3 4 5; 0 0 8 9 10; 0 0 0 14 15]) @test_throws ErrorException AbstractAlgebra.Solve.solve(M, [ R(1) ]) @@ -81,7 +81,7 @@ end @testset "Linear solving context" begin - for R in [ QQ, ZZ ] + for R in [ QQ, ZZ, GF(101) ] M = matrix(R, [1 2 3 4 5; 0 0 8 9 10; 0 0 0 14 15]) C = AbstractAlgebra.Solve.solve_init(M) From fb2c344f0281fb0889733e439821ee31fbab47a9 Mon Sep 17 00:00:00 2001 From: Johannes Schmitt Date: Fri, 15 Mar 2024 17:05:56 +0100 Subject: [PATCH 10/10] More test coverage --- src/Solve.jl | 24 ++++++++++++++++++++++++ test/Solve-test.jl | 3 ++- 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/Solve.jl b/src/Solve.jl index 10fdbccf3b..1e8a361ba4 100644 --- a/src/Solve.jl +++ b/src/Solve.jl @@ -115,6 +115,30 @@ end Return a context object `C` that allows to efficiently solve linear systems $Ax = b$ or $xA = b$ for different $b$. + +# Example + +```jldoctest; setup = :(using AbstractAlgebra) +julia> A = QQ[1 2 3; 0 3 0; 5 0 0]; + +julia> C = solve_init(A) +Linear solving context of matrix + [1//1 2//1 3//1] + [0//1 3//1 0//1] + [5//1 0//1 0//1] + +julia> solve(C, [QQ(1), QQ(1), QQ(1)], side = :left) +3-element Vector{Rational{BigInt}}: + 1//3 + 1//9 + 2//15 + +julia> solve(C, [QQ(1), QQ(1), QQ(1)], side = :right) +3-element Vector{Rational{BigInt}}: + 1//5 + 1//3 + 2//45 +``` """ function solve_init(A::MatElem) return SolveCtx(A) diff --git a/test/Solve-test.jl b/test/Solve-test.jl index 2eeb2ca62a..b600300856 100644 --- a/test/Solve-test.jl +++ b/test/Solve-test.jl @@ -81,7 +81,8 @@ end @testset "Linear solving context" begin - for R in [ QQ, ZZ, GF(101) ] + # QQ, ring, field, fraction field different from QQ to test different implementations + for R in [ QQ, ZZ, GF(101), fraction_field(QQ["x"][1]) ] M = matrix(R, [1 2 3 4 5; 0 0 8 9 10; 0 0 0 14 15]) C = AbstractAlgebra.Solve.solve_init(M)