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

More linear algebra #1821

Merged
merged 10 commits into from
Jul 19, 2024
Merged

More linear algebra #1821

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
eeb19af
Start using traits
joschmitt Jul 3, 2024
c1281f7
Deduplicate code from AbstractAlgebra
joschmitt Jul 10, 2024
fbfa48e
Linear solving over Z/nZ
joschmitt Jul 11, 2024
f89d1d8
Adjust to new type parameter of solve context
joschmitt Jul 17, 2024
6b6c25c
More tests
joschmitt Jul 17, 2024
1e4c222
Bump AbstractAlgebra compat to 0.42
lgoettgens Jul 19, 2024
6819ba3
Bump version to 0.46.0-DEV
lgoettgens Jul 19, 2024
19aa37f
Merge branch 'master' into js/traits
lgoettgens Jul 19, 2024
12441be
More tests
joschmitt Jul 19, 2024
85747e9
Even more tests
joschmitt Jul 19, 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
4 changes: 2 additions & 2 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "Nemo"
uuid = "2edaba10-b0f1-5616-af89-8c11ac63239a"
version = "0.45.7"
version = "0.46.0-DEV"

[deps]
AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
Expand All @@ -13,7 +13,7 @@ RandomExtensions = "fb686558-2515-59ef-acaa-46db3789a887"
SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[compat]
AbstractAlgebra = "0.41.9"
AbstractAlgebra = "0.42.0"
FLINT_jll = "^300.100.100"
Libdl = "1.6"
LinearAlgebra = "1.6"
Expand Down
1 change: 1 addition & 0 deletions src/Nemo.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -178,6 +178,7 @@ import AbstractAlgebra: force_op
import AbstractAlgebra: get_attribute
import AbstractAlgebra: get_cached!
import AbstractAlgebra: Group
import AbstractAlgebra: howell_form!
import AbstractAlgebra: Indent
import AbstractAlgebra: is_terse
import AbstractAlgebra: Lowercase
Expand Down
14 changes: 8 additions & 6 deletions src/arb/ComplexMat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -565,10 +565,13 @@ function _solve_lu_precomp(P::Perm, LU::ComplexMatrix, y::ComplexMatrix)
return z
end

function Solve._can_solve_internal_no_check(A::ComplexMatrix, b::ComplexMatrix, task::Symbol; side::Symbol = :left)
Solve.matrix_normal_form_type(::ComplexField) = Solve.LUTrait()
Solve.matrix_normal_form_type(::ComplexMatrix) = Solve.LUTrait()

function Solve._can_solve_internal_no_check(::Solve.LUTrait, A::ComplexMatrix, b::ComplexMatrix, task::Symbol; side::Symbol = :left)
nrows(A) != ncols(A) && error("Only implemented for square matrices")
if side === :left
fl, sol, K = Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
fl, sol, K = Solve._can_solve_internal_no_check(Solve.LUTrait(), transpose(A), transpose(b), task, side = :right)
return fl, transpose(sol), transpose(K)
end

Expand All @@ -590,9 +593,8 @@ end
#
################################################################################

AbstractAlgebra.solve_context_type(::Type{ComplexFieldElem}) = Solve.SolveCtx{ComplexFieldElem, ComplexMatrix, ComplexMatrix, ComplexMatrix}

function Solve._init_reduce(C::Solve.SolveCtx{ComplexFieldElem})
function Solve._init_reduce(C::Solve.SolveCtx{ComplexFieldElem, Solve.LUTrait})
if isdefined(C, :red) && isdefined(C, :lu_perm)
return nothing
end
Expand All @@ -615,7 +617,7 @@ function Solve._init_reduce(C::Solve.SolveCtx{ComplexFieldElem})
return nothing
end

function Solve._init_reduce_transpose(C::Solve.SolveCtx{ComplexFieldElem})
function Solve._init_reduce_transpose(C::Solve.SolveCtx{ComplexFieldElem, Solve.LUTrait})
if isdefined(C, :red_transp) && isdefined(C, :lu_perm_transp)
return nothing
end
Expand All @@ -638,7 +640,7 @@ function Solve._init_reduce_transpose(C::Solve.SolveCtx{ComplexFieldElem})
return nothing
end

function Solve._can_solve_internal_no_check(C::Solve.SolveCtx{ComplexFieldElem}, b::ComplexMatrix, task::Symbol; side::Symbol = :left)
function Solve._can_solve_internal_no_check(::Solve.LUTrait, C::Solve.SolveCtx{ComplexFieldElem, Solve.LUTrait}, b::ComplexMatrix, task::Symbol; side::Symbol = :left)
if side === :right
LU = Solve.reduced_matrix(C)
p = Solve.lu_permutation(C)
Expand Down
15 changes: 8 additions & 7 deletions src/arb/RealMat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -508,10 +508,13 @@ function _solve_lu_precomp(P::Perm, LU::RealMatrix, y::RealMatrix)
return z
end

function Solve._can_solve_internal_no_check(A::RealMatrix, b::RealMatrix, task::Symbol; side::Symbol = :left)
Solve.matrix_normal_form_type(::RealField) = Solve.LUTrait()
Solve.matrix_normal_form_type(::RealMatrix) = Solve.LUTrait()

function Solve._can_solve_internal_no_check(::Solve.LUTrait, A::RealMatrix, b::RealMatrix, task::Symbol; side::Symbol = :left)
nrows(A) != ncols(A) && error("Only implemented for square matrices")
if side === :left
fl, sol, K = Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
fl, sol, K = Solve._can_solve_internal_no_check(Solve.LUTrait(), transpose(A), transpose(b), task, side = :right)
return fl, transpose(sol), transpose(K)
end

Expand All @@ -533,9 +536,7 @@ end
#
################################################################################

AbstractAlgebra.solve_context_type(::Type{RealFieldElem}) = Solve.SolveCtx{RealFieldElem, RealMatrix, RealMatrix, RealMatrix}

function Solve._init_reduce(C::Solve.SolveCtx{RealFieldElem})
function Solve._init_reduce(C::Solve.SolveCtx{RealFieldElem, Solve.LUTrait})
if isdefined(C, :red) && isdefined(C, :lu_perm)
return nothing
end
Expand All @@ -558,7 +559,7 @@ function Solve._init_reduce(C::Solve.SolveCtx{RealFieldElem})
return nothing
end

function Solve._init_reduce_transpose(C::Solve.SolveCtx{RealFieldElem})
function Solve._init_reduce_transpose(C::Solve.SolveCtx{RealFieldElem, Solve.LUTrait})
if isdefined(C, :red_transp) && isdefined(C, :lu_perm_transp)
return nothing
end
Expand All @@ -581,7 +582,7 @@ function Solve._init_reduce_transpose(C::Solve.SolveCtx{RealFieldElem})
return nothing
end

function Solve._can_solve_internal_no_check(C::Solve.SolveCtx{RealFieldElem}, b::RealMatrix, task::Symbol; side::Symbol = :left)
function Solve._can_solve_internal_no_check(::Solve.LUTrait, C::Solve.SolveCtx{RealFieldElem, Solve.LUTrait}, b::RealMatrix, task::Symbol; side::Symbol = :left)
if side === :right
LU = Solve.reduced_matrix(C)
p = Solve.lu_permutation(C)
Expand Down
15 changes: 8 additions & 7 deletions src/arb/acb_mat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -568,10 +568,13 @@ function _solve_lu_precomp(P::Perm, LU::AcbMatrix, y::AcbMatrix)
return z
end

function Solve._can_solve_internal_no_check(A::AcbMatrix, b::AcbMatrix, task::Symbol; side::Symbol = :left)
Solve.matrix_normal_form_type(::AcbField) = Solve.LUTrait()
Solve.matrix_normal_form_type(::AcbMatrix) = Solve.LUTrait()

function Solve._can_solve_internal_no_check(::Solve.LUTrait, 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 = Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
fl, sol, K = Solve._can_solve_internal_no_check(Solve.LUTrait(), transpose(A), transpose(b), task, side = :right)
return fl, transpose(sol), transpose(K)
end

Expand All @@ -593,9 +596,7 @@ end
#
################################################################################

AbstractAlgebra.solve_context_type(::Type{AcbFieldElem}) = Solve.SolveCtx{AcbFieldElem, AcbMatrix, AcbMatrix, AcbMatrix}

function Solve._init_reduce(C::Solve.SolveCtx{AcbFieldElem})
function Solve._init_reduce(C::Solve.SolveCtx{AcbFieldElem, Solve.LUTrait})
if isdefined(C, :red) && isdefined(C, :lu_perm)
return nothing
end
Expand All @@ -618,7 +619,7 @@ function Solve._init_reduce(C::Solve.SolveCtx{AcbFieldElem})
return nothing
end

function Solve._init_reduce_transpose(C::Solve.SolveCtx{AcbFieldElem})
function Solve._init_reduce_transpose(C::Solve.SolveCtx{AcbFieldElem, Solve.LUTrait})
if isdefined(C, :red_transp) && isdefined(C, :lu_perm_transp)
return nothing
end
Expand All @@ -641,7 +642,7 @@ function Solve._init_reduce_transpose(C::Solve.SolveCtx{AcbFieldElem})
return nothing
end

function Solve._can_solve_internal_no_check(C::Solve.SolveCtx{AcbFieldElem}, b::AcbMatrix, task::Symbol; side::Symbol = :left)
function Solve._can_solve_internal_no_check(::Solve.LUTrait, C::Solve.SolveCtx{AcbFieldElem, Solve.LUTrait}, b::AcbMatrix, task::Symbol; side::Symbol = :left)
if side === :right
LU = Solve.reduced_matrix(C)
p = Solve.lu_permutation(C)
Expand Down
15 changes: 8 additions & 7 deletions src/arb/arb_mat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -535,10 +535,13 @@ function _solve_cholesky_precomp(cho::ArbMatrix, y::ArbMatrix)
return z
end

function Solve._can_solve_internal_no_check(A::ArbMatrix, b::ArbMatrix, task::Symbol; side::Symbol = :left)
Solve.matrix_normal_form_type(::ArbField) = Solve.LUTrait()
Solve.matrix_normal_form_type(::ArbMatrix) = Solve.LUTrait()

function Solve._can_solve_internal_no_check(::Solve.LUTrait, 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 = Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
fl, sol, K = Solve._can_solve_internal_no_check(Solve.LUTrait(), transpose(A), transpose(b), task, side = :right)
return fl, transpose(sol), transpose(K)
end

Expand All @@ -560,9 +563,7 @@ end
#
################################################################################

AbstractAlgebra.solve_context_type(::Type{ArbFieldElem}) = Solve.SolveCtx{ArbFieldElem, ArbMatrix, ArbMatrix, ArbMatrix}

function Solve._init_reduce(C::Solve.SolveCtx{ArbFieldElem})
function Solve._init_reduce(C::Solve.SolveCtx{ArbFieldElem, Solve.LUTrait})
if isdefined(C, :red) && isdefined(C, :lu_perm)
return nothing
end
Expand All @@ -585,7 +586,7 @@ function Solve._init_reduce(C::Solve.SolveCtx{ArbFieldElem})
return nothing
end

function Solve._init_reduce_transpose(C::Solve.SolveCtx{ArbFieldElem})
function Solve._init_reduce_transpose(C::Solve.SolveCtx{ArbFieldElem, Solve.LUTrait})
if isdefined(C, :red_transp) && isdefined(C, :lu_perm_transp)
return nothing
end
Expand All @@ -608,7 +609,7 @@ function Solve._init_reduce_transpose(C::Solve.SolveCtx{ArbFieldElem})
return nothing
end

function Solve._can_solve_internal_no_check(C::Solve.SolveCtx{ArbFieldElem}, b::ArbMatrix, task::Symbol; side::Symbol = :left)
function Solve._can_solve_internal_no_check(::Solve.LUTrait, C::Solve.SolveCtx{ArbFieldElem, Solve.LUTrait}, b::ArbMatrix, task::Symbol; side::Symbol = :left)
if side === :right
LU = Solve.reduced_matrix(C)
p = Solve.lu_permutation(C)
Expand Down
98 changes: 13 additions & 85 deletions src/flint/fmpq_mat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -653,9 +653,20 @@ function _solve_dixon(a::QQMatrix, b::QQMatrix)
return z
end

function Solve._can_solve_internal_no_check(A::QQMatrix, b::QQMatrix, task::Symbol; side::Symbol = :left)
# Actually we use a modular algorithm in flint...
Solve.matrix_normal_form_type(::QQField) = Solve.RREFTrait()
Solve.matrix_normal_form_type(::QQMatrix) = Solve.RREFTrait()

# fflu is much slower in some cases, so we do an rref (with transformation)
# here and let flint choose an algorithm, see
# https://github.com/Nemocas/Nemo.jl/issues/1710.
function Solve.solve_context_type(::QQField)
return Solve.solve_context_type(Solve.RREFTrait(), QQFieldElem)
end

function Solve._can_solve_internal_no_check(::Solve.RREFTrait, A::QQMatrix, b::QQMatrix, task::Symbol; side::Symbol = :left)
if side === :left
fl, sol, K = Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
fl, sol, K = Solve._can_solve_internal_no_check(Solve.RREFTrait(), transpose(A), transpose(b), task, side = :right)
return fl, transpose(sol), transpose(K)
end

Expand All @@ -669,89 +680,6 @@ function Solve._can_solve_internal_no_check(A::QQMatrix, b::QQMatrix, task::Symb
return Bool(fl), x, kernel(A, side = :right)
end

function Solve._init_reduce(C::Solve.SolveCtx{QQFieldElem})
if isdefined(C, :red) && isdefined(C, :trafo)
return nothing
end

# fflu is much slower in some cases, so we do an rref (with transformation)
# here and let flint choose an algorithm, see
# https://github.com/Nemocas/Nemo.jl/issues/1710.
A = matrix(C)
B = hcat(deepcopy(A), identity_matrix(base_ring(A), nrows(A)))
rref!(B)
r = nrows(B)
while r > 0 && is_zero(view(B, r:r, 1:ncols(A)))
r -= 1
end
Solve.set_rank!(C, r)
C.red = view(B, 1:nrows(A), 1:ncols(A))
C.trafo = view(B, 1:nrows(A), ncols(A) + 1:ncols(B))
return nothing
end

function Solve.reduced_matrix(C::Solve.SolveCtx{QQFieldElem})
Solve._init_reduce(C)
return C.red
end

function Solve.transformation_matrix(C::Solve.SolveCtx{QQFieldElem})
Solve._init_reduce(C)
return C.trafo
end

function Solve._init_reduce_transpose(C::Solve.SolveCtx{QQFieldElem})
if isdefined(C, :red_transp) && isdefined(C, :trafo_transp)
return nothing
end

# fflu is much slower in some cases, so we do an rref (with transformation)
# here and let flint choose an algorithm, see
# https://github.com/Nemocas/Nemo.jl/issues/1710.
A = matrix(C)
B = hcat(transpose(A), identity_matrix(base_ring(A), ncols(A)))
rref!(B)
r = nrows(B)
while r > 0 && is_zero(view(B, r:r, 1:nrows(A)))
r -= 1
end
Solve.set_rank!(C, r)
C.red_transp = view(B, 1:ncols(A), 1:nrows(A))
C.trafo_transp = view(B, 1:ncols(A), nrows(A) + 1:ncols(B))
return nothing
end

function Solve.reduced_matrix_of_transpose(C::Solve.SolveCtx{QQFieldElem})
Solve._init_reduce_transpose(C)
return C.red_transp
end

function Solve.transformation_matrix_of_transpose(C::Solve.SolveCtx{QQFieldElem})
Solve._init_reduce_transpose(C)
return C.trafo_transp
end

function Solve._can_solve_internal_no_check(C::Solve.SolveCtx{QQFieldElem}, b::QQMatrix, task::Symbol; side::Symbol = :left)
if side === :right
fl, sol = Solve._can_solve_with_rref(b, Solve.transformation_matrix(C), rank(C), Solve.pivot_and_non_pivot_cols(C), task)
if !fl || task !== :with_kernel
return fl, sol, zero(b, 0, 0)
end

_, K = Solve._kernel_of_rref(Solve.reduced_matrix(C), rank(C), Solve.pivot_and_non_pivot_cols(C))
return fl, sol, K
else# side === :left
fl, sol_transp = Solve._can_solve_with_rref(transpose(b), Solve.transformation_matrix_of_transpose(C), rank(C), Solve.pivot_and_non_pivot_cols_of_transpose(C), task)
sol = transpose(sol_transp)
if !fl || task !== :with_kernel
return fl, sol, zero(b, 0, 0)
end

_, K = Solve._kernel_of_rref(Solve.reduced_matrix_of_transpose(C), rank(C), Solve.pivot_and_non_pivot_cols_of_transpose(C))
return fl, sol, transpose(K)
end
end

###############################################################################
#
# Trace
Expand Down
57 changes: 6 additions & 51 deletions src/flint/fmpz_mat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1220,27 +1220,6 @@ function nullspace(x::ZZMatrix)
return ncols(x), identity_matrix(x, ncols(x))
end

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

if side === :left
K = 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 @@ -1476,9 +1455,12 @@ end
#
###############################################################################

function Solve._can_solve_internal_no_check(A::ZZMatrix, b::ZZMatrix, task::Symbol; side::Symbol = :left)
Solve.matrix_normal_form_type(::ZZRing) = Solve.HermiteFormTrait()
Solve.matrix_normal_form_type(::ZZMatrix) = Solve.HermiteFormTrait()

function Solve._can_solve_internal_no_check(::Solve.HermiteFormTrait, A::ZZMatrix, b::ZZMatrix, task::Symbol; side::Symbol = :left)
if side === :left
fl, sol, K = Solve._can_solve_internal_no_check(transpose(A), transpose(b), task, side = :right)
fl, sol, K = Solve._can_solve_internal_no_check(Solve.HermiteFormTrait(), transpose(A), transpose(b), task, side = :right)
return fl, transpose(sol), transpose(K)
end

Expand Down Expand Up @@ -1529,37 +1511,10 @@ function Solve._can_solve_internal_no_check(A::ZZMatrix, b::ZZMatrix, task::Symb
if task === :with_solution
return true, sol, zero(A, 0, 0)
end
_, K = Solve._kernel_of_hnf(A, H, T)
K = transpose(Solve._kernel_of_hnf(H, T))
return fl, sol, K
end

# Overwrite some solve context functionality so that it uses `transpose` and not
# `lazy_transpose`

AbstractAlgebra.solve_context_type(::Type{ZZRingElem}) = Solve.SolveCtx{ZZRingElem, ZZMatrix, ZZMatrix, ZZMatrix}

function Solve._init_reduce_transpose(C::Solve.SolveCtx{ZZRingElem})
if isdefined(C, :red_transp) && isdefined(C, :trafo_transp)
return nothing
end

R, U = hnf_with_transform(transpose(matrix(C)))
C.red_transp = R
C.trafo_transp = U
return nothing
end

function Solve.kernel(C::Solve.SolveCtx{ZZRingElem}; side::Symbol = :left)
Solve.check_option(side, [:right, :left], "side")

if side === :right
return Solve._kernel_of_hnf(matrix(C), Solve.reduced_matrix_of_transpose(C), Solve.transformation_matrix_of_transpose(C))[2]
else
nullity, X = Solve._kernel_of_hnf(matrix(C), Solve.reduced_matrix(C), Solve.transformation_matrix(C))
return transpose(X)
end
end

Base.reduce(::typeof(hcat), A::AbstractVector{ZZMatrix}) = AbstractAlgebra._hcat(A)

Base.reduce(::typeof(vcat), A::AbstractVector{ZZMatrix}) = AbstractAlgebra._vcat(A)
Expand Down
Loading
Loading