Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use more flint nullspaces + linear solving clean-up #1681

Merged
merged 14 commits into from
Feb 25, 2024
Merged

Use more flint nullspaces + linear solving clean-up #1681

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
4128e24
Use `fmpz_mod_mat_nullspace` in flint
joschmitt Feb 21, 2024
688cae9
Use `fq_default_mat_nullspace` in flint
joschmitt Feb 21, 2024
948c254
Use `fq_mat_nullspace` in flint
joschmitt Feb 21, 2024
869701c
Use `fq_nmod_mat_nullspace` in flint
joschmitt Feb 21, 2024
2db55ac
Use `fmpz_mod_mat_nullspace` in flint
joschmitt Feb 21, 2024
12a285e
Use `nmod_mat_nullspace` in flint
joschmitt Feb 21, 2024
c456066
Adapt `zzModMatrix` nullspace
joschmitt Feb 21, 2024
a15bf26
Remove outdated `kernel` methods
joschmitt Feb 21, 2024
9379a46
Add generic `kernel` for flint types
joschmitt Feb 21, 2024
7355a4e
Import `AbstractAlgebra.Solve`
joschmitt Feb 23, 2024
c3a3d46
Avoid lazy transpose of C types in solve context
joschmitt Feb 23, 2024
42fd431
Remove outdated linear solving methods
joschmitt Feb 23, 2024
351fc6d
Move `eigenspaces` etc to new `matrix.jl` file
joschmitt Feb 23, 2024
05dd769
Clean up `_cansolve`
joschmitt Feb 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
193 changes: 0 additions & 193 deletions src/HeckeMiscMatrix.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -175,104 +175,6 @@ function lift(a::Generic.Mat{ZZModRingElem})
return z
end

################################################################################
#
# Kernel function
#
################################################################################

function _right_kernel(x::fpMatrix)
z = zero_matrix(base_ring(x), ncols(x), max(nrows(x), ncols(x)))
n = ccall((:nmod_mat_nullspace, libflint), Int, (Ref{fpMatrix}, Ref{fpMatrix}), z, x)
return n, z
end

function _left_kernel(x::fpMatrix)
n, M = _right_kernel(transpose(x))
return n, transpose(M)
end

@doc raw"""
_left_kernel(a::ZZMatrix) -> Int, ZZMatrix

It returns a tuple $(n, M)$ where $M$ is a matrix whose rows generate
the kernel of $a$ and $n$ is the rank of the kernel.
"""
function _left_kernel(x::ZZMatrix)
if nrows(x) == 0
return 0, zero(x, 0, 0)
end
x1 = transpose(hnf(transpose(x)))
H, U = hnf_with_transform(x1)
i = 1
for outer i in 1:nrows(H)
if is_zero_row(H, i)
break
end
end
if is_zero_row(H, i)
return nrows(U) - i + 1, view(U, i:nrows(U), 1:ncols(U))
else
return 0, zero_matrix(FlintZZ, 0, ncols(U))
end
end

function _right_kernel(x::ZZMatrix)
n, M = _left_kernel(transpose(x))
return n, transpose(M)
end

function _right_kernel(M::zzModMatrix)
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 n, k
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 = sub(H, 1:nr, 1:nrows(M))
for i = 1:nrows(h)
if is_zero_row(h, i)
k = sub(H, i:nrows(h), nrows(M)+1:ncols(H))
return nrows(k), transpose(k)
end
end
return 0, zero_matrix(R, nrows(M), 0)
end

function _right_kernel(M::ZZModMatrix)
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 = sub(H, 1:nr, 1:nrows(M))
for i = 1:nrows(h)
if is_zero_row(h, i)
k = sub(H, i:nrows(h), nrows(M)+1:ncols(H))
return nrows(k), transpose(k)
end
end
return 0, zero_matrix(R, nrows(M), 0)
end

################################################################################
#
# Reduce the entries of a matrix modulo p
Expand Down Expand Up @@ -712,98 +614,3 @@ function map_entries(R::ZZModRing, M::ZZMatrix)
end
return N
end

################################################################################
#
# Eigenvalues and eigenspaces
#
################################################################################

@doc raw"""
eigenvalues(M::MatElem{T}) where T <: FieldElem

Return the eigenvalues of `M`.
"""
function eigenvalues(M::MatElem{T}) where T <: FieldElem
@assert is_square(M)
K = base_ring(M)
f = charpoly(M)
return roots(f)
end

@doc raw"""
eigenvalues_with_multiplicities(M::MatElem{T}) where T <: FieldElem

Return the eigenvalues of `M` together with their algebraic multiplicities as a
vector of tuples.
"""
function eigenvalues_with_multiplicities(M::MatElem{T}) where T <: FieldElem
@assert is_square(M)
K = base_ring(M)
Kx, x = polynomial_ring(K, "x", cached = false)
f = charpoly(Kx, M)
r = roots(f)
return [ (a, valuation(f, x - a)) for a in r ]
end

@doc raw"""
eigenvalues(L::Field, M::MatElem{T}) where T <: RingElem

Return the eigenvalues of `M` over the field `L`.
"""
function eigenvalues(L::Field, M::MatElem{T}) where T <: RingElem
@assert is_square(M)
M1 = change_base_ring(L, M)
return eigenvalues(M1)
end

@doc raw"""
eigenvalues_with_multiplicities(L::Field, M::MatElem{T}) where T <: RingElem

Return the eigenvalues of `M` over the field `L` together with their algebraic
multiplicities as a vector of tuples.
"""
function eigenvalues_with_multiplicities(L::Field, M::MatElem{T}) where T <: RingElem
@assert is_square(M)
M1 = change_base_ring(L, M)
return eigenvalues_with_multiplicities(M1)
end

@doc raw"""
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 = :left) where T <: FieldElem
@assert is_square(M)
N = deepcopy(M)
for i = 1:ncols(N)
N[i, i] -= lambda
end
return kernel(N, side = side)
end

@doc raw"""
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
corresponding eigenspaces as values.
If side is `:right`, the right eigenspaces are computed, if it is `:left` then the
left eigenspaces are computed.

See also `eigenspace`.
"""
function eigenspaces(M::MatElem{T}; side::Symbol = :left) where T<:FieldElem

S = eigenvalues(M)
L = Dict{elem_type(base_ring(M)), typeof(M)}()
for k in S
push!(L, k => vcat(eigenspace(M, k, side = side)))
end
return L
end
3 changes: 2 additions & 1 deletion src/Nemo.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,6 +196,7 @@ import AbstractAlgebra: promote_rule
import AbstractAlgebra: Ring
import AbstractAlgebra: Set
import AbstractAlgebra: set_attribute!
import AbstractAlgebra: Solve

include("Exports.jl")

Expand Down Expand Up @@ -485,7 +486,7 @@ include("embedding/embedding.jl")

include("Rings.jl")

include("view.jl")
include("matrix.jl")

include("HeckeMiscFiniteField.jl")
include("HeckeMiscInfinity.jl")
Expand Down
28 changes: 10 additions & 18 deletions src/arb/ComplexMat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -574,14 +574,6 @@ function _solve!(z::ComplexMat, x::ComplexMat, y::ComplexMat)
nothing
end

function _solve(x::ComplexMat, y::ComplexMat)
ncols(x) != nrows(x) && error("First argument must be square")
ncols(x) != nrows(y) && error("Matrix dimensions are wrong")
z = similar(y)
_solve!(z, x, y)
return z
end

function _solve_lu_precomp!(z::ComplexMat, P::Generic.Perm, LU::ComplexMat, y::ComplexMat)
Q = inv(P)
ccall((:acb_mat_solve_lu_precomp, libarb), Nothing,
Expand All @@ -597,10 +589,10 @@ 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 = :left)
function 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)
fl, sol, K = Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
return fl, transpose(sol), transpose(K)
end

Expand All @@ -623,10 +615,10 @@ end
################################################################################

function solve_init(A::ComplexMat)
return AbstractAlgebra.Solve.SolveCtx{ComplexFieldElem, ComplexMat, ComplexMat}(A)
return Solve.SolveCtx{ComplexFieldElem, ComplexMat, ComplexMat}(A)
end

function AbstractAlgebra.Solve._init_reduce(C::AbstractAlgebra.Solve.SolveCtx{ComplexFieldElem})
function Solve._init_reduce(C::Solve.SolveCtx{ComplexFieldElem})
if isdefined(C, :red) && isdefined(C, :lu_perm)
return nothing
end
Expand All @@ -649,7 +641,7 @@ function AbstractAlgebra.Solve._init_reduce(C::AbstractAlgebra.Solve.SolveCtx{Co
return nothing
end

function AbstractAlgebra.Solve._init_reduce_transpose(C::AbstractAlgebra.Solve.SolveCtx{ComplexFieldElem})
function Solve._init_reduce_transpose(C::Solve.SolveCtx{ComplexFieldElem})
if isdefined(C, :red_transp) && isdefined(C, :lu_perm_transp)
return nothing
end
Expand All @@ -672,13 +664,13 @@ 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 = :left)
function Solve._can_solve_internal_no_check(C::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)
LU = Solve.reduced_matrix(C)
p = Solve.lu_permutation(C)
else
LU = AbstractAlgebra.Solve.reduced_matrix_of_transpose(C)
p = AbstractAlgebra.Solve.lu_permutation_of_transpose(C)
LU = Solve.reduced_matrix_of_transpose(C)
p = Solve.lu_permutation_of_transpose(C)
b = transpose(b)
end

Expand Down
28 changes: 10 additions & 18 deletions src/arb/RealMat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -516,14 +516,6 @@ function _solve!(z::RealMat, x::RealMat, y::RealMat)
nothing
end

function _solve(x::RealMat, y::RealMat)
ncols(x) != nrows(x) && error("First argument must be square")
ncols(x) != nrows(y) && error("Matrix dimensions are wrong")
z = similar(y)
_solve!(z, x, y)
return z
end

function _solve_lu_precomp!(z::RealMat, P::Generic.Perm, LU::RealMat, y::RealMat)
Q = inv(P)
ccall((:arb_mat_solve_lu_precomp, libarb), Nothing,
Expand All @@ -539,10 +531,10 @@ 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 = :left)
function 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)
fl, sol, K = Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
return fl, transpose(sol), transpose(K)
end

Expand All @@ -565,10 +557,10 @@ end
################################################################################

function solve_init(A::RealMat)
return AbstractAlgebra.Solve.SolveCtx{RealFieldElem, RealMat, RealMat}(A)
return Solve.SolveCtx{RealFieldElem, RealMat, RealMat}(A)
end

function AbstractAlgebra.Solve._init_reduce(C::AbstractAlgebra.Solve.SolveCtx{RealFieldElem})
function Solve._init_reduce(C::Solve.SolveCtx{RealFieldElem})
if isdefined(C, :red) && isdefined(C, :lu_perm)
return nothing
end
Expand All @@ -591,7 +583,7 @@ function AbstractAlgebra.Solve._init_reduce(C::AbstractAlgebra.Solve.SolveCtx{Re
return nothing
end

function AbstractAlgebra.Solve._init_reduce_transpose(C::AbstractAlgebra.Solve.SolveCtx{RealFieldElem})
function Solve._init_reduce_transpose(C::Solve.SolveCtx{RealFieldElem})
if isdefined(C, :red_transp) && isdefined(C, :lu_perm_transp)
return nothing
end
Expand All @@ -614,13 +606,13 @@ 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 = :left)
function Solve._can_solve_internal_no_check(C::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)
LU = Solve.reduced_matrix(C)
p = Solve.lu_permutation(C)
else
LU = AbstractAlgebra.Solve.reduced_matrix_of_transpose(C)
p = AbstractAlgebra.Solve.lu_permutation_of_transpose(C)
LU = Solve.reduced_matrix_of_transpose(C)
p = Solve.lu_permutation_of_transpose(C)
b = transpose(b)
end

Expand Down
Loading
Loading