Skip to content

Commit d6add4a

Browse files
mfasimfasi
committed
A few minor improvements, high order functions removed, added type checking to ensure the precision to be double
1 parent 05eb38b commit d6add4a

File tree

4 files changed

+40
-34
lines changed

4 files changed

+40
-34
lines changed

‎base/linalg/dense.jl

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -281,14 +281,18 @@ end
281281
function logm{T}(A::StridedMatrix{T})
282282
# If possible, use diagonalization
283283
if ishermitian(A)
284-
return funm(log, Hermitian(A))
284+
return logm(Hermitian(A))
285+
end
286+
287+
if T != Float64 && T != Complex{Float64}
288+
warn("The code was originally intended for double precision (Float64, Complex{Float64}) only")
285289
end
286290

287291
# Use Schur decomposition
288292
n = chksquare(A)
289293
S,Q,d = schur(complex(A))
290-
R = full(logm(UpperTriangular(S)))
291-
retmat = Q*R*Q'
294+
R = logm(UpperTriangular(S))
295+
retmat = Q * R * Q'
292296

293297
# Check whether the matrix has nonpositive real eigs
294298
np_real_eigs = false
@@ -313,7 +317,7 @@ logm(a::Complex) = log(a)
313317

314318
function sqrtm{T<:Real}(A::StridedMatrix{T})
315319
if issym(A)
316-
return funm(sqrt, Symmetric(A))
320+
return sqrtm(Symmetric(A))
317321
end
318322
n = chksquare(A)
319323
SchurF = schurfact(complex(A))
@@ -323,7 +327,7 @@ function sqrtm{T<:Real}(A::StridedMatrix{T})
323327
end
324328
function sqrtm{T<:Complex}(A::StridedMatrix{T})
325329
if ishermitian(A)
326-
return funm(sqrt, Hermitian(A))
330+
return sqrtm(Hermitian(A))
327331
end
328332
n = chksquare(A)
329333
SchurF = schurfact(A)

‎base/linalg/symmetric.jl

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -112,8 +112,13 @@ end
112112

113113
#Matrix-valued functions
114114
expm{T<:Real}(A::RealHermSymComplexHerm{T}) = (F = eigfact(A); F.vectors*Diagonal(exp(F.values))*F.vectors')
115-
function funm{T<:Real}(f::Function, A::RealHermSymComplexHerm{T})
115+
function logm{T<:Real}(A::RealHermSymComplexHerm{T})
116116
F = eigfact(A)
117-
isposdef(F) && return F.vectors*Diagonal(f(F.values))*F.vectors'
118-
return F.vectors*Diagonal(f(complex(F.values)))*F.vectors'
117+
isposdef(F) && return F.vectors*Diagonal(log(F.values))*F.vectors'
118+
return F.vectors*Diagonal(log(complex(F.values)))*F.vectors'
119+
end
120+
function sqrtm{T<:Real}(A::RealHermSymComplexHerm{T})
121+
F = eigfact(A)
122+
isposdef(F) && return F.vectors*Diagonal(sqrt(F.values))*F.vectors'
123+
return F.vectors*Diagonal(sqrt(complex(F.values)))*F.vectors'
119124
end

‎base/linalg/triangular.jl

Lines changed: 14 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -895,6 +895,7 @@ end
895895
# Al-Mohy, Higham and Relton, "Computing the Frechet derivative of the matrix
896896
# logarithm and estimating the condition number", SIAM J. Sci. Comput., 35(4),
897897
# (2013), C394-C410.
898+
# http://eprints.ma.man.ac.uk/1687/02/logm.zip
898899
function logm{T}(A0::UpperTriangular{T})
899900
maxsqrt = 100
900901
theta = [1.586970738772063e-005,
@@ -931,9 +932,9 @@ function logm{T}(A0::UpperTriangular{T})
931932
A = sqrtm(A)
932933
end
933934

934-
AmI = A - UpperTriangular(eye(n))
935-
d2 = norm(AmI^2, 1)^(1/2)
936-
d3 = norm(AmI^3, 1)^(1/3)
935+
AmI = UpperTriangular(full(A) - I)
936+
d2 = sqrt(norm(AmI^2, 1))
937+
d3 = cbrt(norm(AmI^3, 1))
937938
alpha2 = max(d2, d3)
938939
foundm = false
939940
if alpha2 <= theta[2]
@@ -944,7 +945,7 @@ function logm{T}(A0::UpperTriangular{T})
944945
while ~foundm
945946
more = false
946947
if s > s0
947-
d3 = norm(AmI^3, 1)^(1/3)
948+
d3 = cbrt(norm(AmI^3, 1))
948949
end
949950
d4 = norm(AmI^4, 1)^(1/4)
950951
alpha3 = max(d3, d4)
@@ -957,7 +958,7 @@ function logm{T}(A0::UpperTriangular{T})
957958
if j <= 6
958959
m = j
959960
break
960-
elseif alpha3/2 <= theta[5] && p < 2
961+
elseif alpha3 / 2 <= theta[5] && p < 2
961962
more = true
962963
p = p + 1
963964
end
@@ -984,12 +985,12 @@ function logm{T}(A0::UpperTriangular{T})
984985
break
985986
end
986987
A = sqrtm(A)
987-
AmI = A - UpperTriangular(eye(n))
988+
AmI = UpperTriangular(full(A) - I)
988989
s = s + 1
989990
end
990991

991992
# Compute accurate superdiagonal of T
992-
p = (1/2^s)
993+
p = 1 / 2^s
993994
for k = 1:n-1
994995
Ak = A0[k,k]
995996
Akp1 = A0[k+1,k+1]
@@ -999,7 +1000,7 @@ function logm{T}(A0::UpperTriangular{T})
9991000
A[k+1,k+1] = Akp1p
10001001
if Ak == Akp1
10011002
A[k,k+1] = p * A0[k,k+1] * Ak^(p-1)
1002-
elseif abs(Ak) < abs(Akp1)/2 || abs(Akp1) < abs(Ak)/2
1003+
elseif 2 * abs(Ak) < abs(Akp1) || 2 * abs(Akp1) < abs(Ak)
10031004
A[k,k+1] = A0[k,k+1] * (Akp1p - Akp) / (Akp1 - Ak)
10041005
else
10051006
logAk = log(Ak)
@@ -1017,7 +1018,7 @@ function logm{T}(A0::UpperTriangular{T})
10171018
r = a - 1
10181019
end
10191020
s0 = s
1020-
if angle(a) >= pi/2
1021+
if angle(a) >= pi / 2
10211022
a = sqrt(a)
10221023
s0 = s - 1
10231024
end
@@ -1034,26 +1035,20 @@ function logm{T}(A0::UpperTriangular{T})
10341035
# Compute the Gauss-Legendre quadrature formula
10351036
R = zeros(Float64, m, m)
10361037
for i = 1:m - 1
1037-
R[i,i+1] = i / sqrt((2 * i).^2-1)
1038+
R[i,i+1] = i / sqrt((2 * i)^2 - 1)
10381039
R[i+1,i] = R[i,i+1]
10391040
end
10401041
x,V = eig(R)
10411042
w = Array(Float64, m)
10421043
for i = 1:m
10431044
x[i] = (x[i] + 1) / 2
1044-
w[i] = (V[1,i])^2
1045+
w[i] = V[1,i]^2
10451046
end
10461047

10471048
# Compute the Padé approximation
1048-
B = zeros(T, n, n)
10491049
Y = zeros(T, n, n)
10501050
for k = 1:m
1051-
for j = 1:n
1052-
B[j,j] = 1. + A[j,j] * x[k]
1053-
for i = 1:j-1
1054-
B[i,j] = A[i,j] * x[k]
1055-
end
1056-
end
1051+
B = scale(x[k], full(A)) + eye(n)
10571052
B = A * inv(B)
10581053
scale!(w[k], B)
10591054
Y = Y + B
@@ -1072,7 +1067,7 @@ function logm{T}(A0::UpperTriangular{T})
10721067
Y[k+1,k+1] = logAkp1
10731068
if Ak == Akp1
10741069
Y[k,k+1] = A0[k,k+1] / Ak
1075-
elseif abs(Ak) < abs(Akp1)/2 || abs(Akp1) < abs(Ak)/2
1070+
elseif 2 * abs(Ak) < abs(Akp1) || 2 * abs(Akp1) < abs(Ak)
10761071
Y[k,k+1] = A0[k,k+1] * (logAkp1 - logAk) / (Akp1 - Ak)
10771072
else
10781073
w = atanh((Akp1 - Ak)/(Akp1 + Ak) + im*pi*(ceil((imag(logAkp1-logAk) - pi)/(2*pi))))

‎test/linalg3.jl

Lines changed: 9 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -154,7 +154,15 @@ for elty in (Float32, Float64, Complex64, Complex128)
154154
0.135335281175235 0.406005843524598 0.541341126763207]')
155155
@test_approx_eq expm(A3) eA3
156156

157-
A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5;
157+
# Hessenberg
158+
@test_approx_eq hessfact(A1)[:H] convert(Matrix{elty},
159+
[4.000000000000000 -1.414213562373094 -1.414213562373095
160+
-1.414213562373095 4.999999999999996 -0.000000000000000
161+
0 -0.000000000000002 3.000000000000000])
162+
end
163+
164+
for elty in (Float64, Complex{Float64})
165+
A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+eps();
158166
1/3 1/4 1/5 1/6;
159167
1/4 1/5 1/6 1/7;
160168
1/5 1/6 1/7 1/8])
@@ -168,12 +176,6 @@ for elty in (Float32, Float64, Complex64, Complex128)
168176

169177
A7 = convert(Matrix{elty}, [1 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1])
170178
@test_approx_eq expm(logm(A7)) A7
171-
172-
# Hessenberg
173-
@test_approx_eq hessfact(A1)[:H] convert(Matrix{elty},
174-
[4.000000000000000 -1.414213562373094 -1.414213562373095
175-
-1.414213562373095 4.999999999999996 -0.000000000000000
176-
0 -0.000000000000002 3.000000000000000])
177179
end
178180

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

0 commit comments

Comments
 (0)