Specialize lu instead of inv

Author: lbenet

src/TaylorSeries.jl

@@ -23,10 +23,12 @@ using Markdown
 using Requires
 
 using LinearAlgebra: norm, mul!,
-    lu, istriu, triu!, istril, tril!, UpperTriangular, LowerTriangular,
-    LinearAlgebra.inv!, LinearAlgebra.checksquare
+    lu, lu!, LinearAlgebra.lutype, LinearAlgebra.copy_oftype,
+    LinearAlgebra.issuccess
+    # istriu, triu!, istril, tril!, UpperTriangular, LowerTriangular,
+    # LinearAlgebra.inv!, LinearAlgebra.checksquare
 
-import LinearAlgebra: norm, mul!
+import LinearAlgebra: norm, mul!, lu
 
 import Base: ==, +, -, *, /, ^
 

src/arithmetic.jl

@@ -566,18 +566,34 @@ end
 # licensed under MIT "Expat".
 # Specialize a method of `inv` for Matrix{Taylor1{T}}. Simply, avoid pivoting,
 # since the polynomial field is not an ordered one.
-function Base.inv(A::StridedMatrix{Taylor1{T}}) where T
-    checksquare(A)
-    S = Taylor1{typeof((one(T)*zero(T) + one(T)*zero(T))/one(T))}
-    AA = convert(AbstractArray{S}, A)
-    if istriu(AA)
-        Ai = triu!(parent(inv(UpperTriangular(AA))))
-    elseif istril(AA)
-        Ai = tril!(parent(inv(LowerTriangular(AA))))
+# function Base.inv(A::StridedMatrix{Taylor1{T}}) where T
+#     checksquare(A)
+#     S = Taylor1{typeof((one(T)*zero(T) + one(T)*zero(T))/one(T))}
+#     AA = convert(AbstractArray{S}, A)
+#     if istriu(AA)
+#         Ai = triu!(parent(inv(UpperTriangular(AA))))
+#     elseif istril(AA)
+#         Ai = tril!(parent(inv(LowerTriangular(AA))))
+#     else
+#         # Do not use pivoting !!
+#         Ai = inv!(lu(AA, Val(false)))
+#         Ai = convert(typeof(parent(Ai)), Ai)
+#     end
+#     return Ai
+# end
+
+# Adapted from (Julia v1.2) stdlib/v1.2/LinearAlgebra/src/lu.jl#240-253
+# and (Julia v1.4.0-dev) stdlib/LinearAlgebra/v1.4/src/lu.jl#270-274,
+# licensed under MIT "Expat".
+# Specialize a method of `lu` for Matrix{Taylor1{T}}, which avoids pivoting,
+# since the polynomial field is not an ordered one.
+# We can't assume an ordered field so we first try without pivoting
+function lu(A::AbstractMatrix{Taylor1{T}}; check::Bool = true) where T
+    S = Taylor1{lutype(T)}
+    F = lu!(copy_oftype(A, S), Val(false); check = false)
+    if issuccess(F)
+        return F
     else
-        # Do not use pivoting !!
-        Ai = inv!(lu(AA, Val(false)))
-        Ai = convert(typeof(parent(Ai)), Ai)
+        return lu!(copy_oftype(A, S), Val(true); check = check)
     end
-    return Ai
 end

test/onevariable.jl

@@ -506,22 +506,30 @@ end
 
 @testset "Test `inv` for `Matrix{Taylor1{Float64}}``" begin
     t = Taylor1(5)
-
     a = Diagonal(rand(0:10,3)) + rand(3, 3)
+    ainv = inv(a)
     b = Taylor1.(a, 5)
     binv = inv(b)
-    @test norm(binv - inv(a), Inf) ≤ 1e-11
-    @test norm(b*binv - I, Inf) ≤ 1e-11
-    @test norm(binv*b - I, Inf) ≤ 1e-11
-    @test norm(triu(b)*inv(UpperTriangular(b)) - I, Inf) ≤ 1e-11
-    @test norm(inv(LowerTriangular(b))*tril(b) - I, Inf) ≤ 1e-11
-
-    b .= b .+ t
-    binv = inv(b)
-    @test norm(b*binv - I, Inf) ≤ 1e-11
-    @test norm(binv*b - I, Inf) ≤ 1e-11
-    @test norm(triu(b)*inv(triu(b)) - I, Inf) ≤ 1e-11
-    @test norm(inv(tril(b))*tril(b) - I, Inf) ≤ 1e-11
+    tol = 1.0e-11
+
+    for its = 1:10
+        a .= Diagonal(rand(2:12,3)) + rand(3, 3)
+        ainv .= inv(a)
+        b .= Taylor1.(a, 5)
+        binv .= inv(b)
+        @test norm(binv - ainv, Inf) ≤ tol
+        @test norm(b*binv - I, Inf) ≤ tol
+        @test norm(binv*b - I, Inf) ≤ tol
+        @test norm(triu(b)*inv(UpperTriangular(b)) - I, Inf) ≤ tol
+        @test norm(inv(LowerTriangular(b))*tril(b) - I, Inf) ≤ tol
+
+        b .= b .+ t
+        binv .= inv(b)
+        @test norm(b*binv - I, Inf) ≤ tol
+        @test norm(binv*b - I, Inf) ≤ tol
+        @test norm(triu(b)*inv(triu(b)) - I, Inf) ≤ tol
+        @test norm(inv(tril(b))*tril(b) - I, Inf) ≤ tol
+    end
 end
 
 @testset "Matrix multiplication for Taylor1" begin