From 557d1615da7735049d67c0ad2b610ec12bcfb744 Mon Sep 17 00:00:00 2001 From: Morgan Rodgers Date: Wed, 26 Jul 2023 11:31:37 +0200 Subject: [PATCH 1/6] Fix `lu!` requiring real matrices to be square. --- src/arb/RealMat.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/arb/RealMat.jl b/src/arb/RealMat.jl index f69fcfe96..c3b959763 100644 --- a/src/arb/RealMat.jl +++ b/src/arb/RealMat.jl @@ -496,7 +496,6 @@ end ############################################################################### function lu!(P::Generic.Perm, x::RealMat) - ncols(x) != nrows(x) && error("Matrix must be square") parent(P).n != nrows(x) && error("Permutation does not match matrix") P.d .-= 1 r = ccall((:arb_mat_lu, libarb), Cint, From 5110bfa8ec36f2c403d3e6fa10bc5db01549414b Mon Sep 17 00:00:00 2001 From: Morgan Rodgers Date: Mon, 31 Jul 2023 14:28:47 +0200 Subject: [PATCH 2/6] add test, fix `lu()` method --- src/arb/RealMat.jl | 11 +++++++---- test/arb/RealMat-test.jl | 12 ++++++++++++ 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/arb/RealMat.jl b/src/arb/RealMat.jl index c3b959763..42e271325 100644 --- a/src/arb/RealMat.jl +++ b/src/arb/RealMat.jl @@ -508,20 +508,23 @@ function lu!(P::Generic.Perm, x::RealMat) end function lu(x::RealMat, P = SymmetricGroup(nrows(x))) + m = nrows(x) + n = ncols(x) + P.n != m && error("Permutation does not match matrix") p = one(P) R = base_ring(x) - L = similar(x) + L = similar(x, R, m, m) U = deepcopy(x) - n = ncols(x) + r = lu!(p, U) - for i = 1:n + for i = 1:m for j = 1:n if i > j L[i, j] = U[i, j] U[i, j] = R() elseif i == j L[i, j] = one(R) - else + elseif j <= m L[i, j] = R() end end diff --git a/test/arb/RealMat-test.jl b/test/arb/RealMat-test.jl index 6a1ed4ecd..035eef789 100644 --- a/test/arb/RealMat-test.jl +++ b/test/arb/RealMat-test.jl @@ -361,6 +361,18 @@ end @test overlaps(B, C) end +@testset "RealMat.lu_nonsquare" begin + S = matrix_space(RR, 2, 3) + + A = S(["1.0 +/- 0.01" "1.0 +/- 0.01" "1.0 +/- 0.01"; + "1.0 +/- 0.01" "0.0 +/- 0.01" "-1.0 +/- 0.01"]) + + r, p, L, U = lu(A) + + @test overlaps(L*U, p*A) + @test r == 2 +end + @testset "RealMat.linear_solving" begin S = matrix_space(RR, 3, 3) T = matrix_space(ZZ, 3, 3) From eb93d9ac84215a7353af43ab20606278fda94108 Mon Sep 17 00:00:00 2001 From: Morgan Rodgers Date: Thu, 3 Aug 2023 14:35:41 +0200 Subject: [PATCH 3/6] fix LU algorithms for other arb types --- src/arb/ComplexMat.jl | 13 +++++++------ src/arb/RealMat.jl | 1 - src/arb/acb_mat.jl | 13 +++++++------ src/arb/arb_mat.jl | 11 ++++++----- 4 files changed, 20 insertions(+), 18 deletions(-) diff --git a/src/arb/ComplexMat.jl b/src/arb/ComplexMat.jl index caa63647e..1a0deeab6 100644 --- a/src/arb/ComplexMat.jl +++ b/src/arb/ComplexMat.jl @@ -566,21 +566,22 @@ function lu!(P::Generic.Perm, x::ComplexMat) end function lu(P::Generic.Perm, x::ComplexMat) - ncols(x) != nrows(x) && error("Matrix must be square") - parent(P).n != nrows(x) && error("Permutation does not match matrix") + m = nrows(x) + n = ncols(x) + P.n != m && error("Permutation does not match matrix") + p = one(P) R = base_ring(x) - L = similar(x) + L = similar(x, R, m, m) U = deepcopy(x) - n = ncols(x) lu!(P, U) - for i = 1:n + for i = 1:m for j = 1:n if i > j L[i, j] = U[i, j] U[i, j] = R() elseif i == j L[i, j] = one(R) - else + elseif j <= m L[i, j] = R() end end diff --git a/src/arb/RealMat.jl b/src/arb/RealMat.jl index 42e271325..3c1f50f6a 100644 --- a/src/arb/RealMat.jl +++ b/src/arb/RealMat.jl @@ -515,7 +515,6 @@ function lu(x::RealMat, P = SymmetricGroup(nrows(x))) R = base_ring(x) L = similar(x, R, m, m) U = deepcopy(x) - r = lu!(p, U) for i = 1:m for j = 1:n diff --git a/src/arb/acb_mat.jl b/src/arb/acb_mat.jl index 0bb18e5fc..92f17dcc4 100644 --- a/src/arb/acb_mat.jl +++ b/src/arb/acb_mat.jl @@ -569,21 +569,22 @@ function lu!(P::Generic.Perm, x::acb_mat) end function lu(P::Generic.Perm, x::acb_mat) - ncols(x) != nrows(x) && error("Matrix must be square") - parent(P).n != nrows(x) && error("Permutation does not match matrix") + m = nrows(x) + n = ncols(x) + P.n != m && error("Permutation does not match matrix") + p = one(P) R = base_ring(x) - L = similar(x) + L = similar(x, R, m, m) U = deepcopy(x) - n = ncols(x) lu!(P, U) - for i = 1:n + for i = 1:m for j = 1:n if i > j L[i, j] = U[i, j] U[i, j] = R() elseif i == j L[i, j] = one(R) - else + elseif j <= m L[i, j] = R() end end diff --git a/src/arb/arb_mat.jl b/src/arb/arb_mat.jl index 413986742..35c9a9604 100644 --- a/src/arb/arb_mat.jl +++ b/src/arb/arb_mat.jl @@ -500,7 +500,6 @@ end ############################################################################### function lu!(P::Generic.Perm, x::arb_mat) - ncols(x) != nrows(x) && error("Matrix must be square") parent(P).n != nrows(x) && error("Permutation does not match matrix") P.d .-= 1 r = ccall((:arb_mat_lu, libarb), Cint, @@ -513,20 +512,22 @@ function lu!(P::Generic.Perm, x::arb_mat) end function lu(x::arb_mat, P = SymmetricGroup(nrows(x))) + m = nrows(x) + n = ncols(x) + P.n != m && error("Permutation does not match matrix") p = one(P) R = base_ring(x) - L = similar(x) + L = similar(x, R, m, m) U = deepcopy(x) - n = ncols(x) r = lu!(p, U) - for i = 1:n + for i = 1:m for j = 1:n if i > j L[i, j] = U[i, j] U[i, j] = R() elseif i == j L[i, j] = one(R) - else + elseif j <= m L[i, j] = R() end end From f8153e6e2e65d28db062f266b282d1975eebc609 Mon Sep 17 00:00:00 2001 From: Morgan Rodgers Date: Thu, 3 Aug 2023 14:50:15 +0200 Subject: [PATCH 4/6] add test for arb_mat --- test/arb/arb_mat-test.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/arb/arb_mat-test.jl b/test/arb/arb_mat-test.jl index e1b1e2f6e..0779b3a54 100644 --- a/test/arb/arb_mat-test.jl +++ b/test/arb/arb_mat-test.jl @@ -361,6 +361,18 @@ end @test overlaps(B, C) end +@testset "RealMat.lu_nonsquare" begin + S = matrix_space(RR, 2, 3) + + A = S(["1.0 +/- 0.01" "1.0 +/- 0.01" "1.0 +/- 0.01"; + "1.0 +/- 0.01" "0.0 +/- 0.01" "-1.0 +/- 0.01"]) + + r, p, L, U = lu(A) + + @test overlaps(L*U, p*A) + @test r == 2 +end + @testset "arb_mat.linear_solving" begin S = matrix_space(RR, 3, 3) T = matrix_space(ZZ, 3, 3) From fbd7c27feca3b01d132012c35b75413ccfe06a2c Mon Sep 17 00:00:00 2001 From: Morgan Rodgers Date: Fri, 4 Aug 2023 10:35:17 +0200 Subject: [PATCH 5/6] tests and using AbstractAlgebra generic `lu()` method --- src/arb/ComplexMat.jl | 24 ------------------------ src/arb/RealMat.jl | 24 ------------------------ src/arb/acb_mat.jl | 24 ------------------------ src/arb/arb_mat.jl | 24 ------------------------ test/arb/ComplexMat-test.jl | 12 ++++++++++++ test/arb/acb_mat-test.jl | 12 ++++++++++++ test/arb/arb_mat-test.jl | 2 +- 7 files changed, 25 insertions(+), 97 deletions(-) diff --git a/src/arb/ComplexMat.jl b/src/arb/ComplexMat.jl index 1a0deeab6..511dcd323 100644 --- a/src/arb/ComplexMat.jl +++ b/src/arb/ComplexMat.jl @@ -565,30 +565,6 @@ function lu!(P::Generic.Perm, x::ComplexMat) return nrows(x) end -function lu(P::Generic.Perm, x::ComplexMat) - m = nrows(x) - n = ncols(x) - P.n != m && error("Permutation does not match matrix") - p = one(P) - R = base_ring(x) - L = similar(x, R, m, m) - U = deepcopy(x) - lu!(P, U) - for i = 1:m - for j = 1:n - if i > j - L[i, j] = U[i, j] - U[i, j] = R() - elseif i == j - L[i, j] = one(R) - elseif j <= m - L[i, j] = R() - end - end - end - return L, U -end - function solve!(z::ComplexMat, x::ComplexMat, y::ComplexMat) r = ccall((:acb_mat_solve, libarb), Cint, (Ref{ComplexMat}, Ref{ComplexMat}, Ref{ComplexMat}, Int), diff --git a/src/arb/RealMat.jl b/src/arb/RealMat.jl index 3c1f50f6a..a51cf2583 100644 --- a/src/arb/RealMat.jl +++ b/src/arb/RealMat.jl @@ -507,30 +507,6 @@ function lu!(P::Generic.Perm, x::RealMat) return nrows(x) end -function lu(x::RealMat, P = SymmetricGroup(nrows(x))) - m = nrows(x) - n = ncols(x) - P.n != m && error("Permutation does not match matrix") - p = one(P) - R = base_ring(x) - L = similar(x, R, m, m) - U = deepcopy(x) - r = lu!(p, U) - for i = 1:m - for j = 1:n - if i > j - L[i, j] = U[i, j] - U[i, j] = R() - elseif i == j - L[i, j] = one(R) - elseif j <= m - L[i, j] = R() - end - end - end - return r, p, L, U -end - function solve!(z::RealMat, x::RealMat, y::RealMat) r = ccall((:arb_mat_solve, libarb), Cint, (Ref{RealMat}, Ref{RealMat}, Ref{RealMat}, Int), diff --git a/src/arb/acb_mat.jl b/src/arb/acb_mat.jl index 92f17dcc4..a606a4b23 100644 --- a/src/arb/acb_mat.jl +++ b/src/arb/acb_mat.jl @@ -568,30 +568,6 @@ function lu!(P::Generic.Perm, x::acb_mat) return nrows(x) end -function lu(P::Generic.Perm, x::acb_mat) - m = nrows(x) - n = ncols(x) - P.n != m && error("Permutation does not match matrix") - p = one(P) - R = base_ring(x) - L = similar(x, R, m, m) - U = deepcopy(x) - lu!(P, U) - for i = 1:m - for j = 1:n - if i > j - L[i, j] = U[i, j] - U[i, j] = R() - elseif i == j - L[i, j] = one(R) - elseif j <= m - L[i, j] = R() - end - end - end - return L, U -end - function solve!(z::acb_mat, x::acb_mat, y::acb_mat) r = ccall((:acb_mat_solve, libarb), Cint, (Ref{acb_mat}, Ref{acb_mat}, Ref{acb_mat}, Int), diff --git a/src/arb/arb_mat.jl b/src/arb/arb_mat.jl index 35c9a9604..dc7ef135d 100644 --- a/src/arb/arb_mat.jl +++ b/src/arb/arb_mat.jl @@ -511,30 +511,6 @@ function lu!(P::Generic.Perm, x::arb_mat) return nrows(x) end -function lu(x::arb_mat, P = SymmetricGroup(nrows(x))) - m = nrows(x) - n = ncols(x) - P.n != m && error("Permutation does not match matrix") - p = one(P) - R = base_ring(x) - L = similar(x, R, m, m) - U = deepcopy(x) - r = lu!(p, U) - for i = 1:m - for j = 1:n - if i > j - L[i, j] = U[i, j] - U[i, j] = R() - elseif i == j - L[i, j] = one(R) - elseif j <= m - L[i, j] = R() - end - end - end - return r, p, L, U -end - function solve!(z::arb_mat, x::arb_mat, y::arb_mat) r = ccall((:arb_mat_solve, libarb), Cint, (Ref{arb_mat}, Ref{arb_mat}, Ref{arb_mat}, Int), diff --git a/test/arb/ComplexMat-test.jl b/test/arb/ComplexMat-test.jl index 2e9922e93..4f20a8148 100644 --- a/test/arb/ComplexMat-test.jl +++ b/test/arb/ComplexMat-test.jl @@ -395,6 +395,18 @@ end @test overlaps(B, C) end +@testset "ComplexMat.lu_nonsquare" begin + S = matrix_space(CC, 2, 3) + + A = S(["1.0 +/- 0.01" "1.0 +/- 0.01" "1.0 +/- 0.01"; + "1.0 +/- 0.01" "0.0 +/- 0.01" "-1.0 +/- 0.01"]) + + r, p, L, U = lu(A) + + @test overlaps(L*U, p*A) + @test r == 2 +end + @testset "ComplexMat.linear_solving" begin S = matrix_space(CC, 3, 3) T = matrix_space(ZZ, 3, 3) diff --git a/test/arb/acb_mat-test.jl b/test/arb/acb_mat-test.jl index 33d268a74..16382210f 100644 --- a/test/arb/acb_mat-test.jl +++ b/test/arb/acb_mat-test.jl @@ -395,6 +395,18 @@ end @test overlaps(B, C) end +@testset "acb_mat.lu_nonsquare" begin + S = matrix_space(CC, 2, 3) + + A = S(["1.0 +/- 0.01" "1.0 +/- 0.01" "1.0 +/- 0.01"; + "1.0 +/- 0.01" "0.0 +/- 0.01" "-1.0 +/- 0.01"]) + + r, p, L, U = lu(A) + + @test overlaps(L*U, p*A) + @test r == 2 +end + @testset "acb_mat.linear_solving" begin S = matrix_space(CC, 3, 3) T = matrix_space(ZZ, 3, 3) diff --git a/test/arb/arb_mat-test.jl b/test/arb/arb_mat-test.jl index 0779b3a54..1afbc2f04 100644 --- a/test/arb/arb_mat-test.jl +++ b/test/arb/arb_mat-test.jl @@ -361,7 +361,7 @@ end @test overlaps(B, C) end -@testset "RealMat.lu_nonsquare" begin +@testset "arb_mat.lu_nonsquare" begin S = matrix_space(RR, 2, 3) A = S(["1.0 +/- 0.01" "1.0 +/- 0.01" "1.0 +/- 0.01"; From da57044abb9ac62b01cee03f386551911d280f8d Mon Sep 17 00:00:00 2001 From: Tommy Hofmann Date: Tue, 15 Aug 2023 19:04:37 +0200 Subject: [PATCH 6/6] Monkey patch --- src/arb/ComplexMat.jl | 2 +- src/arb/RealMat.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/arb/ComplexMat.jl b/src/arb/ComplexMat.jl index 511dcd323..62cbb8f2d 100644 --- a/src/arb/ComplexMat.jl +++ b/src/arb/ComplexMat.jl @@ -562,7 +562,7 @@ function lu!(P::Generic.Perm, x::ComplexMat) r == 0 && error("Could not find $(nrows(x)) invertible pivot elements") P.d .+= 1 inv!(P) - return nrows(x) + return min(nrows(x), ncols(x)) end function solve!(z::ComplexMat, x::ComplexMat, y::ComplexMat) diff --git a/src/arb/RealMat.jl b/src/arb/RealMat.jl index a51cf2583..8517a0a97 100644 --- a/src/arb/RealMat.jl +++ b/src/arb/RealMat.jl @@ -504,7 +504,7 @@ function lu!(P::Generic.Perm, x::RealMat) r == 0 && error("Could not find $(nrows(x)) invertible pivot elements") P.d .+= 1 inv!(P) - return nrows(x) + return min(nrows(x), ncols(x)) end function solve!(z::RealMat, x::RealMat, y::RealMat)