Skip to content

Commit

Permalink
Add more special Solve.* methods (#1668)
Browse files Browse the repository at this point in the history 
* Add `Solve.kernel(::ZZMatrix)`

`hnf(lazy_transpose(A))` would not use the correct HNF method, so we
have to really transpose in this case

* Add `nullspace(::LazyTransposeMatElem{QQFieldElem})`

* Add `Solve.kernel` for some more types (using HeckeMiscMatrix)

* Make `side = :left` the default
  • Loading branch information
@joschmitt
joschmitt authored Feb 10, 2024
1 parent 2107727 commit ae0c3b1
Show file tree
Hide file tree
Showing 26 changed files with 420 additions and 184 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ RandomExtensions = "fb686558-2515-59ef-acaa-46db3789a887"
SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[compat]
AbstractAlgebra = "0.37.0"
AbstractAlgebra = "0.37.6"
Antic_jll = "~0.201.500"
Arb_jll = "~200.2300.000"
Calcium_jll = "~0.401.100"
Expand Down
8 changes: 4 additions & 4 deletions src/HeckeMiscMatrix.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -770,15 +770,15 @@ function eigenvalues_with_multiplicities(L::Field, M::MatElem{T}) where T <: Rin
end

@doc raw"""
eigenspace(M::MatElem{T}, lambda::T; side::Symbol = :right)
eigenspace(M::MatElem{T}, lambda::T; side::Symbol = :left)
where T <: FieldElem -> Vector{MatElem{T}}
Return a basis of the eigenspace of $M$ with respect to the eigenvalue $\lambda$.
If side is `:right`, the right eigenspace is computed, i.e. vectors $v$ such that
$Mv = \lambda v$. If side is `:left`, the left eigenspace is computed, i.e. vectors
$v$ such that $vM = \lambda v$.
"""
function eigenspace(M::MatElem{T}, lambda::T; side::Symbol = :right) where T <: FieldElem
function eigenspace(M::MatElem{T}, lambda::T; side::Symbol = :left) where T <: FieldElem
@assert is_square(M)
N = deepcopy(M)
for i = 1:ncols(N)
Expand All @@ -788,7 +788,7 @@ function eigenspace(M::MatElem{T}, lambda::T; side::Symbol = :right) where T <:
end

@doc raw"""
eigenspaces(M::MatElem{T}; side::Symbol = :right)
eigenspaces(M::MatElem{T}; side::Symbol = :left)
where T <: FieldElem -> Dict{T, MatElem{T}}
Return a dictionary containing the eigenvalues of $M$ as keys and bases for the
Expand All @@ -798,7 +798,7 @@ left eigenspaces are computed.
See also `eigenspace`.
"""
function eigenspaces(M::MatElem{T}; side::Symbol = :right) where T<:FieldElem
function eigenspaces(M::MatElem{T}; side::Symbol = :left) where T<:FieldElem

S = eigenvalues(M)
L = Dict{elem_type(base_ring(M)), typeof(M)}()
Expand Down
4 changes: 2 additions & 2 deletions src/arb/ComplexMat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -597,7 +597,7 @@ function solve_lu_precomp(P::Generic.Perm, LU::ComplexMat, y::ComplexMat)
return z
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::ComplexMat, b::ComplexMat, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::ComplexMat, b::ComplexMat, task::Symbol; side::Symbol = :left)
nrows(A) != ncols(A) && error("Only implemented for square matrices")
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
Expand Down Expand Up @@ -672,7 +672,7 @@ function AbstractAlgebra.Solve._init_reduce_transpose(C::AbstractAlgebra.Solve.S
return nothing
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(C::AbstractAlgebra.Solve.SolveCtx{ComplexFieldElem}, b::ComplexMat, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(C::AbstractAlgebra.Solve.SolveCtx{ComplexFieldElem}, b::ComplexMat, task::Symbol; side::Symbol = :left)
if side === :right
LU = AbstractAlgebra.Solve.reduced_matrix(C)
p = AbstractAlgebra.Solve.lu_permutation(C)
Expand Down
4 changes: 2 additions & 2 deletions src/arb/RealMat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -539,7 +539,7 @@ function solve_lu_precomp(P::Generic.Perm, LU::RealMat, y::RealMat)
return z
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::RealMat, b::RealMat, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::RealMat, b::RealMat, task::Symbol; side::Symbol = :left)
nrows(A) != ncols(A) && error("Only implemented for square matrices")
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
Expand Down Expand Up @@ -614,7 +614,7 @@ function AbstractAlgebra.Solve._init_reduce_transpose(C::AbstractAlgebra.Solve.S
return nothing
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(C::AbstractAlgebra.Solve.SolveCtx{RealFieldElem}, b::RealMat, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(C::AbstractAlgebra.Solve.SolveCtx{RealFieldElem}, b::RealMat, task::Symbol; side::Symbol = :left)
if side === :right
LU = AbstractAlgebra.Solve.reduced_matrix(C)
p = AbstractAlgebra.Solve.lu_permutation(C)
Expand Down
4 changes: 2 additions & 2 deletions src/arb/acb_mat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -600,7 +600,7 @@ function solve_lu_precomp(P::Generic.Perm, LU::AcbMatrix, y::AcbMatrix)
return z
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::AcbMatrix, b::AcbMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::AcbMatrix, b::AcbMatrix, task::Symbol; side::Symbol = :left)
nrows(A) != ncols(A) && error("Only implemented for square matrices")
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
Expand Down Expand Up @@ -675,7 +675,7 @@ function AbstractAlgebra.Solve._init_reduce_transpose(C::AbstractAlgebra.Solve.S
return nothing
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(C::AbstractAlgebra.Solve.SolveCtx{AcbFieldElem}, b::AcbMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(C::AbstractAlgebra.Solve.SolveCtx{AcbFieldElem}, b::AcbMatrix, task::Symbol; side::Symbol = :left)
if side === :right
LU = AbstractAlgebra.Solve.reduced_matrix(C)
p = AbstractAlgebra.Solve.lu_permutation(C)
Expand Down
4 changes: 2 additions & 2 deletions src/arb/arb_mat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -566,7 +566,7 @@ function solve_cholesky_precomp(cho::ArbMatrix, y::ArbMatrix)
return z
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::ArbMatrix, b::ArbMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::ArbMatrix, b::ArbMatrix, task::Symbol; side::Symbol = :left)
nrows(A) != ncols(A) && error("Only implemented for square matrices")
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
Expand Down Expand Up @@ -641,7 +641,7 @@ function AbstractAlgebra.Solve._init_reduce_transpose(C::AbstractAlgebra.Solve.S
return nothing
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(C::AbstractAlgebra.Solve.SolveCtx{ArbFieldElem}, b::ArbMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(C::AbstractAlgebra.Solve.SolveCtx{ArbFieldElem}, b::ArbMatrix, task::Symbol; side::Symbol = :left)
if side === :right
LU = AbstractAlgebra.Solve.reduced_matrix(C)
p = AbstractAlgebra.Solve.lu_permutation(C)
Expand Down
13 changes: 12 additions & 1 deletion src/flint/fmpq_mat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -704,7 +704,7 @@ function can_solve(a::QQMatrix, b::QQMatrix; side::Symbol = :right)
return fl
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::QQMatrix, b::QQMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::QQMatrix, b::QQMatrix, task::Symbol; side::Symbol = :left)
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
return fl, transpose(sol), transpose(K)
Expand Down Expand Up @@ -1070,3 +1070,14 @@ function nullspace(A::QQMatrix)

return nullity, NQQ
end

# For compatibility with the generic `kernel` method in AbstractAlgebra:
# Otherwise the generic nullspace would be used for a "lazy transposed QQMatrix"
# instead of flint.
function nullspace(A::AbstractAlgebra.Solve.LazyTransposeMatElem{QQFieldElem, QQMatrix})
B = transpose(AbstractAlgebra.Solve.data(A))
n, N = nullspace(B)
# Looks worse than it is: we have to transpose here, the lazy_transpose is
# just for type stability
return n, AbstractAlgebra.Solve.lazy_transpose(transpose(N))
end
25 changes: 23 additions & 2 deletions src/flint/fmpz_mat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1110,6 +1110,27 @@ function nullspace(x::ZZMatrix)
return ncols(x), identity_matrix(x, ncols(x))
end

function AbstractAlgebra.Solve.kernel(A::ZZMatrix; side::Symbol = :left)
AbstractAlgebra.Solve.check_option(side, [:right, :left], "side")

if side === :left
K = AbstractAlgebra.Solve.kernel(transpose(A), side = :right)
return transpose(K)
end

A = transpose(hnf(A))
H, U = hnf_with_transform(A)
r = nrows(H)
while r > 0 && is_zero_row(H, r)
r -= 1
end
if is_zero(r)
return transpose(U)
else
return transpose(view(U, r + 1:nrows(U), 1:ncols(U)))
end
end

@doc raw"""
nullspace_right_rational(x::ZZMatrix)
Expand Down Expand Up @@ -1293,7 +1314,7 @@ function cansolve(a::ZZMatrix, b::ZZMatrix)
return true, transpose(z*T)
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::ZZMatrix, b::ZZMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::ZZMatrix, b::ZZMatrix, task::Symbol; side::Symbol = :left)
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
return fl, transpose(sol), transpose(K)
Expand All @@ -1303,7 +1324,7 @@ function AbstractAlgebra.Solve._can_solve_internal_no_check(A::ZZMatrix, b::ZZMa
if task === :only_check || task === :with_solution
return fl, sol, zero(A, 0, 0)
end
return fl, sol, AbstractAlgebra.Solve.kernel(A)
return fl, sol, AbstractAlgebra.Solve.kernel(A, side = :right)
end

Base.reduce(::typeof(hcat), A::AbstractVector{ZZMatrix}) = AbstractAlgebra._hcat(A)
Expand Down
39 changes: 37 additions & 2 deletions src/flint/fmpz_mod_mat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -500,7 +500,7 @@ end
=#

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::ZZModMatrix, b::ZZModMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::ZZModMatrix, b::ZZModMatrix, task::Symbol; side::Symbol = :left)
check_parent(A, b)
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
Expand All @@ -513,7 +513,7 @@ function AbstractAlgebra.Solve._can_solve_internal_no_check(A::ZZModMatrix, b::Z
if task === :only_check || task === :with_solution
return Bool(fl), x, zero(A, 0, 0)
end
return Bool(fl), x, AbstractAlgebra.Solve.kernel(A)
return Bool(fl), x, AbstractAlgebra.Solve.kernel(A, side = :right)
end

################################################################################
Expand Down Expand Up @@ -916,3 +916,38 @@ function matrix_space(R::ZZModRing, r::Int, c::Int; cached::Bool = true)
ZZModMatrixSpace(R, r, c)
end

################################################################################
#
# Kernel
#
################################################################################

function AbstractAlgebra.Solve.kernel(M::ZZModMatrix; side::Symbol = :left)
AbstractAlgebra.Solve.check_option(side, [:right, :left], "side")

if side === :left
K = AbstractAlgebra.Solve.kernel(transpose(M), side = :right)
return transpose(K)
end

R = base_ring(M)
N = hcat(transpose(M), identity_matrix(R, ncols(M)))
if nrows(N) < ncols(N)
N = vcat(N, zero_matrix(R, ncols(N) - nrows(N), ncols(N)))
end
howell_form!(N)
H = N
nr = 1
while nr <= nrows(H) && !is_zero_row(H, nr)
nr += 1
end
nr -= 1
h = view(H, 1:nr, 1:nrows(M))
for i = 1:nrows(h)
if is_zero_row(h, i)
k = view(H, i:nrows(h), nrows(M) + 1:ncols(H))
return transpose(k)
end
end
return zero_matrix(R, ncols(M), 0)
end
4 changes: 2 additions & 2 deletions src/flint/fq_default_mat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -450,7 +450,7 @@ function can_solve(a::FqMatrix, b::FqMatrix; side::Symbol = :right)
return fl
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::FqMatrix, b::FqMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::FqMatrix, b::FqMatrix, task::Symbol; side::Symbol = :left)
check_parent(A, b)
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
Expand All @@ -464,7 +464,7 @@ function AbstractAlgebra.Solve._can_solve_internal_no_check(A::FqMatrix, b::FqMa
if task === :only_check || task === :with_solution
return Bool(fl), x, zero(A, 0, 0)
end
return Bool(fl), x, AbstractAlgebra.Solve.kernel(A)
return Bool(fl), x, AbstractAlgebra.Solve.kernel(A, side = :right)
end

################################################################################
Expand Down
4 changes: 2 additions & 2 deletions src/flint/fq_mat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -457,7 +457,7 @@ function can_solve(a::FqPolyRepMatrix, b::FqPolyRepMatrix; side::Symbol = :right
return fl
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::FqPolyRepMatrix, b::FqPolyRepMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::FqPolyRepMatrix, b::FqPolyRepMatrix, task::Symbol; side::Symbol = :left)
check_parent(A, b)
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
Expand All @@ -471,7 +471,7 @@ function AbstractAlgebra.Solve._can_solve_internal_no_check(A::FqPolyRepMatrix,
if task === :only_check || task === :with_solution
return Bool(fl), x, zero(A, 0, 0)
end
return Bool(fl), x, AbstractAlgebra.Solve.kernel(A)
return Bool(fl), x, AbstractAlgebra.Solve.kernel(A, side = :right)
end

################################################################################
Expand Down
4 changes: 2 additions & 2 deletions src/flint/fq_nmod_mat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -445,7 +445,7 @@ function can_solve(a::fqPolyRepMatrix, b::fqPolyRepMatrix; side::Symbol = :right
return fl
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::fqPolyRepMatrix, b::fqPolyRepMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::fqPolyRepMatrix, b::fqPolyRepMatrix, task::Symbol; side::Symbol = :left)
check_parent(A, b)
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
Expand All @@ -459,7 +459,7 @@ function AbstractAlgebra.Solve._can_solve_internal_no_check(A::fqPolyRepMatrix,
if task === :only_check || task === :with_solution
return Bool(fl), x, zero(A, 0, 0)
end
return Bool(fl), x, AbstractAlgebra.Solve.kernel(A)
return Bool(fl), x, AbstractAlgebra.Solve.kernel(A, side = :right)
end

################################################################################
Expand Down
23 changes: 21 additions & 2 deletions src/flint/gfp_mat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -312,7 +312,7 @@ function can_solve(a::fpMatrix, b::fpMatrix; side::Symbol = :right)
return fl
end

function AbstractAlgebra.Solve._can_solve_internal_no_check(A::fpMatrix, b::fpMatrix, task::Symbol; side::Symbol = :right)
function AbstractAlgebra.Solve._can_solve_internal_no_check(A::fpMatrix, b::fpMatrix, task::Symbol; side::Symbol = :left)
check_parent(A, b)
if side === :left
fl, sol, K = AbstractAlgebra.Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
Expand All @@ -326,7 +326,7 @@ function AbstractAlgebra.Solve._can_solve_internal_no_check(A::fpMatrix, b::fpMa
if task === :only_check || task === :with_solution
return Bool(fl), x, zero(A, 0, 0)
end
return Bool(fl), x, AbstractAlgebra.Solve.kernel(A)
return Bool(fl), x, AbstractAlgebra.Solve.kernel(A, side = :right)
end

################################################################################
Expand Down Expand Up @@ -508,3 +508,22 @@ function matrix_space(R::fpField, r::Int, c::Int; cached::Bool = true)
# TODO/FIXME: `cached` is ignored and only exists for backwards compatibility
fpMatrixSpace(R, r, c)
end

################################################################################
#
# Kernel
#
################################################################################

function AbstractAlgebra.Solve.kernel(A::fpMatrix; side::Symbol = :left)
AbstractAlgebra.Solve.check_option(side, [:right, :left], "side")

if side === :left
K = AbstractAlgebra.Solve.kernel(transpose(A), side = :right)
return transpose(K)
end

K = zero_matrix(base_ring(A), ncols(A), max(nrows(A), ncols(A)))
n = ccall((:nmod_mat_nullspace, libflint), Int, (Ref{fpMatrix}, Ref{fpMatrix}), K, A)
return view(K, 1:nrows(K), 1:n)
end
40 changes: 40 additions & 0 deletions src/flint/nmod_mat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -919,3 +919,43 @@ function matrix_space(R::zzModRing, r::Int, c::Int; cached::Bool = true)
zzModMatrixSpace(R, r, c)
end

################################################################################
#
# Kernel
#
################################################################################

function AbstractAlgebra.Solve.kernel(M::zzModMatrix; side::Symbol = :left)
AbstractAlgebra.Solve.check_option(side, [:right, :left], "side")

if side === :left
K = AbstractAlgebra.Solve.kernel(transpose(M), side = :right)
return transpose(K)
end

R = base_ring(M)
if is_prime(modulus(R))
k = zero_matrix(R, ncols(M), ncols(M))
n = ccall((:nmod_mat_nullspace, libflint), Int, (Ref{zzModMatrix}, Ref{zzModMatrix}), k, M)
return view(k, 1:nrows(k), 1:n)
end

H = hcat(transpose(M), identity_matrix(R, ncols(M)))
if nrows(H) < ncols(H)
H = vcat(H, zero_matrix(R, ncols(H) - nrows(H), ncols(H)))
end
howell_form!(H)
nr = 1
while nr <= nrows(H) && !is_zero_row(H, nr)
nr += 1
end
nr -= 1
h = view(H, 1:nr, 1:nrows(M))
for i = 1:nrows(h)
if is_zero_row(h, i)
k = view(H, i:nrows(h), nrows(M) + 1:ncols(H))
return transpose(k)
end
end
return zero_matrix(R, ncols(M), 0)
end
Loading

2 comments on commit ae0c3b1

@thofma
Copy link
Member

@thofma thofma commented on ae0c3b1 Feb 10, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/100614

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.41.5 -m "<description of version>" ae0c3b1b27e0bb0f0536c55c01304a562b10564a
git push origin v0.41.5

Please sign in to comment.