Separate copying of A from log computation

ElOceanografo · Mar 10, 2021 · 23becc5 · 23becc5
1 parent 37d0344
commit 23becc5
Showing 1 changed file with 30 additions and 20 deletions.
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
@@ -1797,6 +1797,34 @@ log(A::LowerTriangular) = copy(transpose(log(copy(transpose(A)))))
 log(A::UnitLowerTriangular) = copy(transpose(log(copy(transpose(A)))))
 
 function log_quasitriu(A0::AbstractMatrix{T}) where T<:BlasFloat
+    n = checksquare(A0)
+
+    # copy A0, making it real if the logarithm will be real
+    if A0 isa UnitUpperTriangular
+        Ap = (T <: Complex && isreal(A0)) ? real(parent(A0)) : copy(parent(A0))
+        @inbounds for i in 1:n
+            Ap[i,i] = 1
+        end
+        A = UpperTriangular(Ap)
+    elseif isreal(A0) && (!istriu(A0) || !any(x -> real(x) < zero(real(T)), diag(A0)))
+        A = T <: Complex ? real(A0) : copy(A0)
+    else
+        A = complex(A0)
+    end
+
+    Y0 = _log_quasitriu(A)
+
+    # return complex result for complex input
+    Y = T <: Complex ? complex(Y0) : Y0
+
+    if A0 isa UpperTriangular || A0 isa UnitUpperTriangular
+        return UpperTriangular(Y)
+    else
+        return Y
+    end
+end
+
+function _log_quasitriu(A0)
     maxsqrt = 100
     theta = [1.586970738772063e-005,
          2.313807884242979e-003,
@@ -1807,21 +1835,10 @@ function log_quasitriu(A0::AbstractMatrix{T}) where T<:BlasFloat
          2.879093714241194e-001]
     tmax = size(theta, 1)
     n = size(A0, 1)
-    # allocate real A if log(A) will be real and complex A otherwise
-    if isreal(A0) && (!istriu(A0) || !any(x -> real(x) < zero(real(T)), diag(A0)))
-        A = eltype(A0) <: Complex ? real(A0) : copy(A0)
-    else
-        A = complex(A0)
-    end
-    if A0 isa UnitUpperTriangular
-        A = UpperTriangular(parent(A))
-        @inbounds for i in 1:n
-            A[i,i] = 1
-        end
-    end
     p = 0
     m = 0
 
+    A = A0
     # Compute repeated roots
     d = complex(diag(A))
     dm1 = d .- 1
@@ -1926,14 +1943,7 @@ function log_quasitriu(A0::AbstractMatrix{T}) where T<:BlasFloat
     # Compute accurate diagonal and superdiagonal of log(A)
     _log_diag_quasitriu!(Y, A0)
 
-    # return complex result for complex input
-    Yc = eltype(A0) <: Complex ? complex(Y) : Y
-
-    if A0 isa UpperTriangular || A0 isa UnitUpperTriangular
-        return UpperTriangular(Yc)
-    else
-        return Yc
-    end
+    return Y
 end
 
 # Auxiliary functions for matrix logarithm and matrix power