Skip to content

Commit

Permalink
fix: linear solving over rings of integers (#1651)
Browse files Browse the repository at this point in the history
- Take into account, that the pseudo-HNF might have zero rows
  • Loading branch information
thofma authored Oct 18, 2024
1 parent 276ac56 commit 6b7d0e0
Show file tree
Hide file tree
Showing 2 changed files with 25 additions and 12 deletions.
9 changes: 7 additions & 2 deletions src/NumFieldOrd/NfOrd/LinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -871,10 +871,15 @@ function _contained_in_span_of_pseudohnf(v::Generic.Mat{T}, P::PMat{T, S}, ::Val
for i = start:step:stop
# find pivot
if shape === :upperright
piv = findfirst(k -> !iszero(P.matrix[i, k]), 1:ncols(P))::Int
_piv = findfirst(k -> !iszero(P.matrix[i, k]), 1:ncols(P))
else
piv = findlast(k -> !iszero(P.matrix[i, k]), 1:ncols(P))::Int
_piv = findlast(k -> !iszero(P.matrix[i, k]), 1:ncols(P))
end
# The pseudo-HNF might be the zero matrix?
if _piv isa Nothing
continue
end
piv = _piv::Int
if !(w[1, piv]//P.matrix[i, piv] in P.coeffs[i])
if with_solution
return false, sol
Expand Down
28 changes: 18 additions & 10 deletions test/LinearAlgebra/SolveIntegral.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,9 @@
for b in (matrix(OK, 1, 1, [alpha + beta]), OK.([alpha + beta]))
b = matrix(OK, 1, 1, [alpha + beta])
fl, X, K = can_solve_with_solution_and_kernel(M, b; side = :left)
@assert fl
@assert X * M == b
@assert is_zero(K * M)
@test fl
@test X * M == b
@test is_zero(K * M)
end
for b in (matrix(OK, 1, 1, [one(OK)]), [OK(1)])
fl = can_solve(M, b; side = :left)
Expand All @@ -29,9 +29,9 @@
@test ncols(K) >= 2
b = matrix(OK, 1, 1, [alpha + beta])
fl, X, K = can_solve_with_solution_and_kernel(M, b; side = :right)
@assert fl
@assert M * X == b
@assert is_zero(M * K)
@test fl
@test M * X == b
@test is_zero(M * K)
b = matrix(OK, 1, 1, [one(OK)])
fl = can_solve(M, b; side = :right)
@test !fl
Expand All @@ -46,13 +46,21 @@
X = rand(matrix_space(OK, k, n), 5)
B = X * A
fl, v = can_solve_with_solution(A, B, side = :left)
@assert fl
@assert v * A == B
@test fl
@test v * A == B

X = rand(matrix_space(OK, m, k), 5)
B = A * X
fl, v = can_solve_with_solution(A, B, side = :right)
@assert fl
@assert A * v == B
@test fl
@test A * v == B
end

# fix some issue with matrices not of full rank

A = zero_matrix(OK, 1, 1)
B = zero_matrix(OK, 2, 1)
fl, v = can_solve_with_solution(A, B, side = :left)
@test fl
@test v * A == B
end

0 comments on commit 6b7d0e0

Please sign in to comment.