LU-factorize without pivoting for inv of Taylor1 matrices

src/TaylorSeries.jl

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

src/arithmetic.jl

@@ -560,3 +560,40 @@ function mul!(y::Vector{Taylor1{T}},
 
  return y
 end
+
+
+# Adapted from (Julia v1.2) stdlib/v1.2/LinearAlgebra/src/dense.jl#721-734,
+# 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))))
+# 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
+ return lu!(copy_oftype(A, S), Val(true); check = check)
+ end
+end

test/onevariable.jl

@@ -502,14 +502,34 @@ eeuler = Base.MathConstants.e
  @test Taylor1{Float64}(false) == Taylor1([0.0])
  @test Taylor1{Int}(true) == Taylor1([1])
  @test Taylor1{Int}(false) == Taylor1([0])
+end
 
- # Test use of `inv` for computation of matrix inverse of `Matrix{Taylor1{Float64}}`
- a = rand(3, 3)
+@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)
- I_t1_5 = Taylor1.(Matrix{Float64}(I, size(b)), 5) # 5x5 Taylor1{Float64} identity matrix, order 5
- @test norm(b*binv - I_t1_5, Inf) ≤ 1e-12
- @test norm(binv*b - I_t1_5, Inf) ≤ 1e-12
+ 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