Skip to content

Commit

Permalink
Fixed #12294 and added logm for a few classes of special matrices
Browse files Browse the repository at this point in the history
  • Loading branch information
mfasi committed Jul 24, 2015
1 parent 71566da commit 6fdad80
Show file tree
Hide file tree
Showing 4 changed files with 15 additions and 5 deletions.
11 changes: 8 additions & 3 deletions base/linalg/dense.jl
Expand Up @@ -286,9 +286,14 @@ function logm(A::StridedMatrix)

# Use Schur decomposition
n = chksquare(A)
S,Q,d = schur(complex(A))
R = logm(UpperTriangular(S))
retmat = Q * R * Q'
if istriu(A)
retmat = full(logm(UpperTriangular(complex(A))))
d = diag(A)
else
S,Q,d = schur(complex(A))
R = logm(UpperTriangular(S))
retmat = Q * R * Q'
end

# Check whether the matrix has nonpositive real eigs
np_real_eigs = false
Expand Down
1 change: 1 addition & 0 deletions base/linalg/diagonal.jl
Expand Up @@ -119,6 +119,7 @@ end
eye{T}(::Type{Diagonal{T}}, n::Int) = Diagonal(ones(T,n))

expm(D::Diagonal) = Diagonal(exp(D.diag))
logm(D::Diagonal) = Diagonal(log(D.diag))
sqrtm(D::Diagonal) = Diagonal(sqrt(D.diag))

#Linear solver
Expand Down
5 changes: 3 additions & 2 deletions base/linalg/triangular.jl
Expand Up @@ -910,7 +910,7 @@ function logm{T<:Union{Float64,Complex{Float64}}}(A0::UpperTriangular{T})
2.879093714241194e-001]
tmax = size(theta, 1)
n = size(A0, 1)
A = A0
A = copy(A0)
p = 0
m = 0

Expand All @@ -919,7 +919,7 @@ function logm{T<:Union{Float64,Complex{Float64}}}(A0::UpperTriangular{T})
dm1 = Array(T, n)
s = 0
for i = 1:n
dm1[i] = dm1[i] - 1.
dm1[i] = d[i] - 1.
end
while norm(dm1, Inf) > theta[tmax]
for i = 1:n
Expand Down Expand Up @@ -1078,6 +1078,7 @@ function logm{T<:Union{Float64,Complex{Float64}}}(A0::UpperTriangular{T})
return UpperTriangular(Y)

end
logm(A::LowerTriangular) = logm(A.').'

function sqrtm{T}(A::UpperTriangular{T})
n = size(A, 1)
Expand Down
3 changes: 3 additions & 0 deletions test/linalg3.jl
Expand Up @@ -182,6 +182,9 @@ for elty in (Float64, Complex{Float64})
@test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.")
redirect_stderr(OLD_STDERR)

A7 = convert(Matrix{elty}, [1 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1])
@test_approx_eq expm(logm(A7)) A7

end

A8 = 100 * [-1+1im 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1];
Expand Down

0 comments on commit 6fdad80

Please sign in to comment.