From 407b0e0b5d3dc0d7ae9a7a1eb99b14b688e0a1a8 Mon Sep 17 00:00:00 2001 From: Tommy Hofmann Date: Wed, 16 Aug 2023 10:56:54 +0200 Subject: [PATCH] Add is_invertible_with_inverse for arb matrices (#1525) --- src/arb/ComplexMat.jl | 14 ++++++++++---- src/arb/RealMat.jl | 10 ++++++++-- src/arb/acb_mat.jl | 10 ++++++++-- src/arb/arb_mat.jl | 10 ++++++++-- test/arb/ComplexMat-test.jl | 8 ++++++++ test/arb/RealMat-test.jl | 8 ++++++++ test/arb/acb_mat-test.jl | 8 ++++++++ test/arb/arb_mat-test.jl | 8 ++++++++ 8 files changed, 66 insertions(+), 10 deletions(-) diff --git a/src/arb/ComplexMat.jl b/src/arb/ComplexMat.jl index caa63647e..4522d529a 100644 --- a/src/arb/ComplexMat.jl +++ b/src/arb/ComplexMat.jl @@ -440,12 +440,18 @@ Given a $n\times n$ matrix of type `acb_mat`, return an $n\times n$ matrix $X$ such that $AX$ contains the identity matrix. If $A$ cannot be inverted numerically an exception is raised. """ -function inv(x::ComplexMat, prec::Int = precision(Balls)) - ncols(x) != nrows(x) && error("Matrix must be square") +function inv(x::ComplexMat) + fl, z = is_invertible_with_inverse(x) + fl && return z + error("Matrix singular or cannot be inverted numerically") +end + +function is_invertible_with_inverse(x::ComplexMat) + ncols(x) != nrows(x) && return false, x z = similar(x) r = ccall((:acb_mat_inv, libarb), Cint, - (Ref{ComplexMat}, Ref{ComplexMat}, Int), z, x, prec) - Bool(r) ? (return z) : error("Matrix cannot be inverted numerically") + (Ref{ComplexMat}, Ref{ComplexMat}, Int), z, x, precision(Balls)) + return Bool(r), z end ############################################################################### diff --git a/src/arb/RealMat.jl b/src/arb/RealMat.jl index f69fcfe96..ea08f7a5b 100644 --- a/src/arb/RealMat.jl +++ b/src/arb/RealMat.jl @@ -398,11 +398,17 @@ $n\times n$ matrix $X$ such that $AX$ contains the identity matrix. If $A$ cannot be inverted numerically an exception is raised. """ function inv(x::RealMat) - ncols(x) != nrows(x) && error("Matrix must be square") + fl, z = is_invertible_with_inverse(x) + fl && return z + error("Matrix singular or cannot be inverted numerically") +end + +function is_invertible_with_inverse(x::RealMat) + ncols(x) != nrows(x) && return false, x z = similar(x) r = ccall((:arb_mat_inv, libarb), Cint, (Ref{RealMat}, Ref{RealMat}, Int), z, x, precision(Balls)) - Bool(r) ? (return z) : error("Matrix cannot be inverted numerically") + return Bool(r), z end ############################################################################### diff --git a/src/arb/acb_mat.jl b/src/arb/acb_mat.jl index 0bb18e5fc..d0557baff 100644 --- a/src/arb/acb_mat.jl +++ b/src/arb/acb_mat.jl @@ -444,11 +444,17 @@ $n\times n$ matrix $X$ such that $AX$ contains the identity matrix. If $A$ cannot be inverted numerically an exception is raised. """ function inv(x::acb_mat) - ncols(x) != nrows(x) && error("Matrix must be square") + fl, z = is_invertible_with_inverse(x) + fl && return z + error("Matrix singular or cannot be inverted numerically") +end + +function is_invertible_with_inverse(x::acb_mat) + ncols(x) != nrows(x) && return false, x z = similar(x) r = ccall((:acb_mat_inv, libarb), Cint, (Ref{acb_mat}, Ref{acb_mat}, Int), z, x, precision(base_ring(x))) - Bool(r) ? (return z) : error("Matrix cannot be inverted numerically") + return Bool(r), z end ############################################################################### diff --git a/src/arb/arb_mat.jl b/src/arb/arb_mat.jl index 413986742..eab79a27b 100644 --- a/src/arb/arb_mat.jl +++ b/src/arb/arb_mat.jl @@ -402,11 +402,17 @@ $n\times n$ matrix $X$ such that $AX$ contains the identity matrix. If $A$ cannot be inverted numerically an exception is raised. """ function inv(x::arb_mat) - ncols(x) != nrows(x) && error("Matrix must be square") + fl, z = is_invertible_with_inverse(x) + fl && return z + error("Matrix singular or cannot be inverted numerically") +end + +function is_invertible_with_inverse(x::arb_mat) + ncols(x) != nrows(x) && return false, x z = similar(x) r = ccall((:arb_mat_inv, libarb), Cint, (Ref{arb_mat}, Ref{arb_mat}, Int), z, x, precision(base_ring(x))) - Bool(r) ? (return z) : error("Matrix cannot be inverted numerically") + return Bool(r), z end ############################################################################### diff --git a/test/arb/ComplexMat-test.jl b/test/arb/ComplexMat-test.jl index 2e9922e93..1a854cfe8 100644 --- a/test/arb/ComplexMat-test.jl +++ b/test/arb/ComplexMat-test.jl @@ -325,6 +325,14 @@ end @test overlaps(A*C, one(S)) @test contains(C, B) + + fl, C = is_invertible_with_inverse(A) + @test fl && contains(C, B) + + A = CC[1 1; 1 1] + fl, C = is_invertible_with_inverse(A) + @test !fl + @test_throws ErrorException inv(A) end @testset "ComplexMat.divexact" begin diff --git a/test/arb/RealMat-test.jl b/test/arb/RealMat-test.jl index 6a1ed4ecd..b07eb779b 100644 --- a/test/arb/RealMat-test.jl +++ b/test/arb/RealMat-test.jl @@ -294,6 +294,14 @@ end @test overlaps(A*C, one(S)) @test contains(C, B) + + fl, C = is_invertible_with_inverse(A) + @test fl && contains(C, B) + + A = RR[1 1; 1 1] + fl, C = is_invertible_with_inverse(A) + @test !fl + @test_throws ErrorException inv(A) end @testset "RealMat.divexact" begin diff --git a/test/arb/acb_mat-test.jl b/test/arb/acb_mat-test.jl index 33d268a74..e7fcad142 100644 --- a/test/arb/acb_mat-test.jl +++ b/test/arb/acb_mat-test.jl @@ -325,6 +325,14 @@ end @test overlaps(A*C, one(S)) @test contains(C, B) + + fl, C = is_invertible_with_inverse(A) + @test fl && contains(C, B) + + A = CC[1 1; 1 1] + fl, C = is_invertible_with_inverse(A) + @test !fl + @test_throws ErrorException inv(A) end @testset "acb_mat.divexact" begin diff --git a/test/arb/arb_mat-test.jl b/test/arb/arb_mat-test.jl index e1b1e2f6e..60d882f79 100644 --- a/test/arb/arb_mat-test.jl +++ b/test/arb/arb_mat-test.jl @@ -294,6 +294,14 @@ end @test overlaps(A*C, one(S)) @test contains(C, B) + + fl, C = is_invertible_with_inverse(A) + @test fl && contains(C, B) + + A = RR[1 1; 1 1] + fl, C = is_invertible_with_inverse(A) + @test !fl + @test_throws ErrorException inv(A) end @testset "arb_mat.divexact" begin