Don't ccall arb_mat_lu

Nemocas · Oct 23, 2024 · 34f3c41 · 34f3c41
1 parent c05ecf2
commit 34f3c41
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 32 deletions.
diff --git a/src/arb/RealMat.jl b/src/arb/RealMat.jl
@@ -467,18 +467,22 @@ end
 #
 ###############################################################################
 
-function lu!(P::Perm, x::RealMatrix)
+function lu!(P::Perm, z::RealMatrix, x::RealMatrix)
   parent(P).n != nrows(x) && error("Permutation does not match matrix")
   P.d .-= 1
   r = ccall((:arb_mat_lu, libflint), Cint,
             (Ptr{Int}, Ref{RealMatrix}, Ref{RealMatrix}, Int),
-            P.d, x, x, precision(Balls))
+            P.d, z, x, precision(Balls))
   r == 0 && error("Could not find $(nrows(x)) invertible pivot elements")
   P.d .+= 1
   inv!(P)
   return min(nrows(x), ncols(x))
 end
 
+function lu!(P::Perm, x::RealMatrix)
+  return lu!(P, x, x)
+end
+
 function _solve!(z::RealMatrix, x::RealMatrix, y::RealMatrix)
   r = ccall((:arb_mat_solve, libflint), Cint,
             (Ref{RealMatrix}, Ref{RealMatrix}, Ref{RealMatrix}, Int),
@@ -537,13 +541,7 @@ function Solve._init_reduce(C::Solve.SolveCtx{RealFieldElem, Solve.LUTrait})
   A = matrix(C)
   P = Perm(nrows(C))
   x = similar(A, nrows(A), ncols(A))
-  P.d .-= 1
-  fl = ccall((:arb_mat_lu, libflint), Cint,
-             (Ptr{Int}, Ref{RealMatrix}, Ref{RealMatrix}, Int),
-             P.d, x, A, precision(Balls))
-  fl == 0 && error("Could not find $(nrows(x)) invertible pivot elements")
-  P.d .+= 1
-  inv!(P)
+  lu!(P, x, A)
 
   C.red = x
   C.lu_perm = P
@@ -560,13 +558,7 @@ function Solve._init_reduce_transpose(C::Solve.SolveCtx{RealFieldElem, Solve.LUT
   A = transpose(matrix(C))
   P = Perm(nrows(C))
   x = similar(A, nrows(A), ncols(A))
-  P.d .-= 1
-  fl = ccall((:arb_mat_lu, libflint), Cint,
-             (Ptr{Int}, Ref{RealMatrix}, Ref{RealMatrix}, Int),
-             P.d, x, A, precision(Balls))
-  fl == 0 && error("Could not find $(nrows(x)) invertible pivot elements")
-  P.d .+= 1
-  inv!(P)
+  lu!(P, x, A)
 
   C.red_transp = x
   C.lu_perm_transp = P

diff --git a/src/arb/arb_mat.jl b/src/arb/arb_mat.jl
@@ -480,18 +480,22 @@ function cholesky(x::ArbMatrix)
   return z
 end
 
-function lu!(P::Perm, x::ArbMatrix)
+function lu!(P::Perm, z::ArbMatrix, x::ArbMatrix)
   parent(P).n != nrows(x) && error("Permutation does not match matrix")
   P.d .-= 1
   r = ccall((:arb_mat_lu, libflint), Cint,
             (Ptr{Int}, Ref{ArbMatrix}, Ref{ArbMatrix}, Int),
-            P.d, x, x, precision(base_ring(x)))
+            P.d, z, x, precision(base_ring(x)))
   r == 0 && error("Could not find $(nrows(x)) invertible pivot elements")
   P.d .+= 1
   inv!(P)
   return nrows(x)
 end
 
+function lu!(P::Perm, x::ArbMatrix)
+  return lu!(P, x, x)
+end
+
 function _solve!(z::ArbMatrix, x::ArbMatrix, y::ArbMatrix)
   r = ccall((:arb_mat_solve, libflint), Cint,
             (Ref{ArbMatrix}, Ref{ArbMatrix}, Ref{ArbMatrix}, Int),
@@ -564,13 +568,7 @@ function Solve._init_reduce(C::Solve.SolveCtx{ArbFieldElem, Solve.LUTrait})
   A = matrix(C)
   P = Perm(nrows(C))
   x = similar(A, nrows(A), ncols(A))
-  P.d .-= 1
-  fl = ccall((:arb_mat_lu, libflint), Cint,
-             (Ptr{Int}, Ref{ArbMatrix}, Ref{ArbMatrix}, Int),
-             P.d, x, A, precision(base_ring(A)))
-  fl == 0 && error("Could not find $(nrows(x)) invertible pivot elements")
-  P.d .+= 1
-  inv!(P)
+  lu!(P, x, A)
 
   C.red = x
   C.lu_perm = P
@@ -587,13 +585,7 @@ function Solve._init_reduce_transpose(C::Solve.SolveCtx{ArbFieldElem, Solve.LUTr
   A = transpose(matrix(C))
   P = Perm(nrows(C))
   x = similar(A, nrows(A), ncols(A))
-  P.d .-= 1
-  fl = ccall((:arb_mat_lu, libflint), Cint,
-             (Ptr{Int}, Ref{ArbMatrix}, Ref{ArbMatrix}, Int),
-             P.d, x, A, precision(base_ring(A)))
-  fl == 0 && error("Could not find $(nrows(x)) invertible pivot elements")
-  P.d .+= 1
-  inv!(P)
+  lu!(P, x, A)
 
   C.red_transp = x
   C.lu_perm_transp = P