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

Jump to bottom

feat: Change solving interface #1601

Merged
merged 1 commit into from
Feb 12, 2024
Merged

feat: Change solving interface #1601

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
16 changes: 5 additions & 11 deletions docs/src/linear_solving.md
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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 @@
# 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 @@
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 @@
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 @@
return r1 + r2
end

function solve_triu(T::MatElem, b::MatElem; cutoff::Int = cutoff)
function _solve_triu(T::MatElem, b::MatElem; cutoff::Int = cutoff)

Check warning on line 274 in src/Matrix-Strassen.jl

View check run for this annotation

Codecov / codecov/patch

src/Matrix-Strassen.jl#L274 

Added line #L274 was not covered by tests
#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)

Check warning on line 278 in src/Matrix-Strassen.jl

View check run for this annotation

Codecov / codecov/patch

src/Matrix-Strassen.jl#L278 

Added line #L278 was not covered by tests
return R
end

Expand All @@ -292,16 +292,16 @@
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)

Check warning on line 296 in src/Matrix-Strassen.jl

View check run for this annotation

Codecov / codecov/patch

src/Matrix-Strassen.jl#L295-L296 

Added lines #L295 - L296 were not covered by tests

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

Check warning on line 300 in src/Matrix-Strassen.jl

View check run for this annotation

Codecov / codecov/patch

src/Matrix-Strassen.jl#L300 

Added line #L300 was not covered by tests

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

Check warning on line 304 in src/Matrix-Strassen.jl

View check run for this annotation

Codecov / codecov/patch

src/Matrix-Strassen.jl#L304 

Added line #L304 was not covered by tests

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