Skip to content

Commit

Permalink
feat: Change solving interface (Nemocas#1601)
Browse files Browse the repository at this point in the history 
- Use the `Solve.*` functionality
- Hide all previous `*solve*` and `*kernel*` functions using
  underscores

[breaking]
  • Loading branch information
@thofma
thofma authored Feb 12, 2024
1 parent cf05e04 commit f8e2306
Show file tree
Hide file tree
Showing 15 changed files with 218 additions and 393 deletions.
16 changes: 5 additions & 11 deletions docs/src/linear_solving.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,7 @@
CurrentModule = AbstractAlgebra.Solve
```

# Linear solving

!!! note

This functionality is experimental and subject to change.
# [Linear solving](@id solving_chapter)

## Overview of the functionality

Expand All @@ -17,23 +13,21 @@ The module `AbstractAlgebra.Solve` provides the following four functions for sol
* `can_solve_with_solution_and_kernel`

All of these take the same set of arguments, namely:
* a matrix $A$ of type `MatElem{T}`;
* a vector or matrix $B$ of type `Vector{T}` or `MatElem{T}`;
* a matrix $A$ of type `MatElem`;
* a vector or matrix $B$ of type `Vector` or `MatElem`;
* a keyword argument `side` which can be either `:left` (default) or `:right`.

If `side` is `:left`, the system $xA = B$ is solved, otherwise the system $Ax = B$ is solved.
For matrices defined over a field, the functions internally rely on `rref`.
If the matrices are defined over a ring, the function `hnf_with_transform` is required internally.

The functionality of the functions can be summarized as follows.
* `solve`: return a solution, if it exists, otherwise throw an error.
* `can_solve`: return `true`, if a solution exists, `false` otherwise.
* `can_solve_with_solution`: return `true` and a solution, if this exists, and `false` and an empty vector or matrix otherwise.
* `can_solve_with_solution_and_kernel`: like `can_solve_with_solution` and additionally return a matrix whose columns (respectively rows) give a basis of the kernel of $A$.
* `can_solve_with_solution_and_kernel`: like `can_solve_with_solution` and additionally return a matrix whose rows (respectively columns) give a basis of the kernel of $A$.

## Solving with several right hand sides

Systems $Ax = b_1,\dots, Ax = b_k$ with the same matrix $A$, but several right hand sides $b_i$ can be solved more efficiently, by first initializing a "context object" `C` of type `SolveCtx`.
Systems $xA = b_1,\dots, xA = b_k$ with the same matrix $A$, but several right hand sides $b_i$ can be solved more efficiently, by first initializing a "context object" `C`.
```@docs
solve_init
```
Expand Down
170 changes: 1 addition & 169 deletions docs/src/matrix.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -959,143 +959,7 @@ julia> pfaffians(M, 2)

### Linear solving

```@docs
solve{T <: FieldElem}(::MatElem{T}, ::MatElem{T})
```

```@docs
solve_rational{T <: RingElem}(::MatElem{T}, ::MatElem{T})
```

```@docs
can_solve_with_solution{T <: RingElement}(::MatElem{T}, ::MatElem{T})
```

```@docs
can_solve{T <: RingElement}(::MatElem{T}, ::MatElem{T})
```

```@docs
solve_left{T <: RingElem}(::MatElem{T}, ::MatElem{T})
```

```@docs
solve_triu{T <: FieldElem}(::MatElem{T}, ::MatElem{T}, ::Bool)
```

```@docs
can_solve_left_reduced_triu{T <: RingElement}(::MatElem{T}, ::MatElem{T})
```

**Examples**

```jldoctest
julia> R, x = polynomial_ring(QQ, "x")
(Univariate polynomial ring in x over rationals, x)
julia> K, = residue_field(R, x^3 + 3x + 1); a = K(x);
julia> S = matrix_space(K, 3, 3)
Matrix space of 3 rows and 3 columns
over residue field of univariate polynomial ring modulo x^3 + 3*x + 1
julia> U = matrix_space(K, 3, 1)
Matrix space of 3 rows and 1 column
over residue field of univariate polynomial ring modulo x^3 + 3*x + 1
julia> A = S([K(0) 2a + 3 a^2 + 1; a^2 - 2 a - 1 2a; a^2 + 3a + 1 2a K(1)])
[ 0 2*x + 3 x^2 + 1]
[ x^2 - 2 x - 1 2*x]
[x^2 + 3*x + 1 2*x 1]
julia> b = U([2a, a + 1, (-a - 1)])
[ 2*x]
[ x + 1]
[-x - 1]
julia> x = solve(A, b)
[ 1984//7817*x^2 + 1573//7817*x - 937//7817]
[ -2085//7817*x^2 + 1692//7817*x + 965//7817]
[-3198//7817*x^2 + 3540//7817*x - 3525//7817]
julia> A = matrix(ZZ, 2, 2, [1, 2, 0, 2])
[1 2]
[0 2]
julia> b = matrix(ZZ, 2, 1, [2, 1])
[2]
[1]
julia> can_solve(A, b, side = :right)
false
julia> A = matrix(QQ, 2, 2, [3, 4, 5, 6])
[3//1 4//1]
[5//1 6//1]
julia> b = matrix(QQ, 1, 2, [2, 1])
[2//1 1//1]
julia> can_solve_with_solution(A, b; side = :left)
(true, [-7//2 5//2])
julia> A = S([a + 1 2a + 3 a^2 + 1; K(0) a^2 - 1 2a; K(0) K(0) a])
[x + 1 2*x + 3 x^2 + 1]
[ 0 x^2 - 1 2*x]
[ 0 0 x]
julia> bb = U([2a, a + 1, (-a - 1)])
[ 2*x]
[ x + 1]
[-x - 1]
julia> x = solve_triu(A, bb, false)
[ 1//3*x^2 + 8//3*x + 13//3]
[-3//5*x^2 - 3//5*x - 12//5]
[ x^2 + 2]
julia> R, x = polynomial_ring(ZZ, "x")
(Univariate polynomial ring in x over integers, x)
julia> S = matrix_space(R, 3, 3)
Matrix space of 3 rows and 3 columns
over univariate polynomial ring in x over integers
julia> U = matrix_space(R, 3, 2)
Matrix space of 3 rows and 2 columns
over univariate polynomial ring in x over integers
julia> A = S([R(0) 2x + 3 x^2 + 1; x^2 - 2 x - 1 2x; x^2 + 3x + 1 2x R(1)])
[ 0 2*x + 3 x^2 + 1]
[ x^2 - 2 x - 1 2*x]
[x^2 + 3*x + 1 2*x 1]
julia> bbb = U(transpose([2x x + 1 (-x - 1); x + 1 (-x) x^2]))
[ 2*x x + 1]
[ x + 1 -x]
[-x - 1 x^2]
julia> x, d = solve_rational(A, bbb)
([3*x^4-10*x^3-8*x^2-11*x-4 -x^5+3*x^4+x^3-2*x^2+3*x-1; -2*x^5-x^4+6*x^3+2*x+1 x^6+x^5+4*x^4+9*x^3+8*x^2+5*x+2; 6*x^4+12*x^3+15*x^2+6*x-3 -2*x^5-4*x^4-6*x^3-9*x^2-4*x+1], x^5 + 2*x^4 + 15*x^3 + 18*x^2 + 8*x + 7)
julia> S = matrix_space(ZZ, 3, 3)
Matrix space of 3 rows and 3 columns
over integers
julia> T = matrix_space(ZZ, 3, 1)
Matrix space of 3 rows and 1 column
over integers
julia> A = S([BigInt(2) 3 5; 1 4 7; 9 2 2])
[2 3 5]
[1 4 7]
[9 2 2]
julia> B = T([BigInt(4), 5, 7])
[4]
[5]
[7]
```
See [Linear Solving & Kernel](@ref solving_chapter)

### Inverse

Expand Down Expand Up @@ -1162,38 +1026,6 @@ julia> X, d = pseudo_inv(A)
nullspace{T <: FieldElem}(::MatElem{T})
```

### Kernel

```@docs
kernel{T <: RingElem}(::MatElem{T})
left_kernel{T <: RingElem}(::MatElem{T})
right_kernel{T <: RingElem}(::MatElem{T})
```

**Examples**

```jldoctest
julia> S = matrix_space(ZZ, 4, 4)
Matrix space of 4 rows and 4 columns
over integers
julia> M = S([1 2 0 4;
2 0 1 1;
0 1 1 -1;
2 -1 0 2])
[1 2 0 4]
[2 0 1 1]
[0 1 1 -1]
[2 -1 0 2]
julia> nr, Nr = kernel(M)
(1, [-8; -6; 11; 5])
julia> nl, Nl = left_kernel(M)
(1, [0 -1 1 1])
```

### Hessenberg form

```@docs
Expand Down
15 changes: 5 additions & 10 deletions src/AbstractAlgebra.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -696,6 +696,11 @@ import .Generic: sturm_sequence

include("Solve.jl")

import ..Solve: solve
import ..Solve: solve_init
import ..Solve: can_solve
import ..Solve: can_solve_with_solution

# Do not export inv, div, divrem, exp, log, sqrt, numerator and denominator as we define our own
export _check_dim
export _checkbounds
Expand Down Expand Up @@ -727,10 +732,7 @@ export basis
export block_diagonal_matrix
export cached
export can_solve
export can_solve_left_reduced_triu
export can_solve_with_kernel
export can_solve_with_solution
export can_solve_with_solution_interpolation
export canonical_unit
export change_base_ring
export change_coefficient_ring
Expand Down Expand Up @@ -920,7 +922,6 @@ export leading_exponent_vector
export leading_exponent_word
export leading_monomial
export leading_term
export left_kernel
export leglength
export length
export lift
Expand Down Expand Up @@ -1057,7 +1058,6 @@ export reverse_cols
export reverse_cols!
export reverse_rows
export reverse_rows!
export right_kernel
export rising_factorial
export rising_factorial2
export rowlength
Expand Down Expand Up @@ -1088,11 +1088,6 @@ export snf_kb_with_transform
export snf_kb!
export snf_with_transform
export solve
export solve_ff
export solve_left
export solve_rational
export solve_triu
export solve_with_det
export sort_terms!
export SparsePolynomialRing
export Strassen
Expand Down
6 changes: 3 additions & 3 deletions src/Generic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,9 @@ import Base: //
import Base: /
import Base: !=

# The type and helper function for the dictionaries for hashing
import ..AbstractAlgebra: _can_solve_with_solution_fflu
import ..AbstractAlgebra: _can_solve_with_solution_lu
import ..AbstractAlgebra: _left_kernel
import ..AbstractAlgebra: @attributes
import ..AbstractAlgebra: @enable_all_show_via_expressify
import ..AbstractAlgebra: add!
Expand All @@ -125,8 +127,6 @@ import ..AbstractAlgebra: addmul!
import ..AbstractAlgebra: base_ring
import ..AbstractAlgebra: base_ring_type
import ..AbstractAlgebra: CacheDictType
import ..AbstractAlgebra: can_solve_with_solution_fflu
import ..AbstractAlgebra: can_solve_with_solution_lu
import ..AbstractAlgebra: canonical_unit
import ..AbstractAlgebra: change_base_ring
import ..AbstractAlgebra: characteristic
Expand Down
28 changes: 14 additions & 14 deletions src/Matrix-Strassen.jl
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
"""
Provides generic asymptotically fast matrix methods:
- mul and mul! using the Strassen scheme
- solve_tril!
- _solve_tril!
- lu!
- solve_triu
- _solve_triu
Just prefix the function by "Strassen." all 4 functions support a keyword
argument "cutoff" to indicate when the base case should be used.
Expand Down Expand Up @@ -171,9 +171,9 @@ end
# B C Y V
# => X = solve(A, U)
# Y = solve(C, V - B*X)
function solve_tril!(A::MatElem{T}, B::MatElem{T}, C::MatElem{T}, f::Int = 0; cutoff::Int = 2) where T
function _solve_tril!(A::MatElem{T}, B::MatElem{T}, C::MatElem{T}, f::Int = 0; cutoff::Int = 2) where T
if nrows(A) < cutoff || ncols(A) < cutoff
return AbstractAlgebra.solve_tril!(A, B, C, f)
return AbstractAlgebra._solve_tril!(A, B, C, f)
end
n = nrows(B)
n2 = div(n, 2)
Expand All @@ -184,10 +184,10 @@ function solve_tril!(A::MatElem{T}, B::MatElem{T}, C::MatElem{T}, f::Int = 0; cu
X2 = view(A, n2+1:n, 1:ncols(A))
C1 = view(C, 1:n2, 1:ncols(A))
C2 = view(C, n2+1:n, 1:ncols(A))
solve_tril!(X1, B11, C1, f; cutoff)
_solve_tril!(X1, B11, C1, f; cutoff)
x = B21 * X1 # strassen...
sub!(X2, C2, x)
solve_tril!(X2, B22, X2, f; cutoff)
_solve_tril!(X2, B22, X2, f; cutoff)
end

function apply!(A::MatElem, P::Perm{Int}; offset::Int = 0)
Expand Down Expand Up @@ -249,8 +249,8 @@ function lu!(P::Perm{Int}, A; cutoff::Int = 300)
if r1 > 0
#Note: A00 is a view of A0 thus a view of A
# A0 is lu!, thus implicitly two triangular matrices giving the
# lu decomosition. solve_tril! looks ONLY at the lower part of A00
solve_tril!(A01, A00, A01, 1)
# lu decomosition. _solve_tril! looks ONLY at the lower part of A00
_solve_tril!(A01, A00, A01, 1)
X = A10 * A01
sub!(A11, A11, X)
end
Expand All @@ -271,11 +271,11 @@ function lu!(P::Perm{Int}, A; cutoff::Int = 300)
return r1 + r2
end

function solve_triu(T::MatElem, b::MatElem; cutoff::Int = cutoff)
function _solve_triu(T::MatElem, b::MatElem; cutoff::Int = cutoff)
#b*inv(T), thus solves Tx = b for T upper triangular
n = ncols(T)
if n <= cutoff
R = AbstractAlgebra.solve_triu(T, b)
R = AbstractAlgebra._solve_triu(T, b)
return R
end

Expand All @@ -292,16 +292,16 @@ function solve_triu(T::MatElem, b::MatElem; cutoff::Int = cutoff)
B = view(T, 1:n2, 1+n2:n)
C = view(T, 1+n2:n, 1+n2:n)

S = solve_triu(A, U; cutoff)
R = solve_triu(A, X; cutoff)
S = _solve_triu(A, U; cutoff)
R = _solve_triu(A, X; cutoff)

SS = mul(S, B; cutoff)
sub!(SS, V, SS)
SS = solve_triu(C, SS; cutoff)
SS = _solve_triu(C, SS; cutoff)

RR = mul(R, B; cutoff)
sub!(RR, Y, RR)
RR = solve_triu(C, RR; cutoff)
RR = _solve_triu(C, RR; cutoff)

return [S SS; R RR]
end
Expand Down
Loading

0 comments on commit f8e2306

Please sign in to comment.