Skip to content

Commit 8a2ea49

Browse files
andyferrisandyferris
authored
Merge pull request #424 from schmrlng/pull-request/7289dd37
Utilize lu for det, inv, solve, expm
2 parents 64cf1cb + da34ded commit 8a2ea49

File tree

11 files changed

+313
-19
lines changed

11 files changed

+313
-19
lines changed

‎perf/bench_matrix_ops.txt

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
Julia Version 0.6.2
2+
Commit d386e40c17* (2017-12-13 18:08 UTC)
3+
Platform Info:
4+
OS: Linux (x86_64-linux-gnu)
5+
CPU: AMD Ryzen 7 1800X Eight-Core Processor
6+
WORD_SIZE: 64
7+
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Barcelona)
8+
LAPACK: libopenblas64_
9+
LIBM: libopenlibm
10+
LLVM: libLLVM-3.9.1 (ORCJIT, generic)
11+
80×6 DataFrames.DataFrame
12+
│ Row │ SIZE │ STAT │ det │ inv │ expm │ \\ │
13+
├─────┼──────┼───────────────────┼─────────────┼─────────────┼─────────────┼─────────────┤
14+
│ 1 │ 1 │ compile time (s) │ 0.0358709 │ 0.00410504 │ 0.029984 │ 0.0932384 │
15+
│ 2 │ │ StaticArrays (μs) │ 0.001813 │ 0.001813 │ 0.0106907 │ 0.004648 │
16+
│ 3 │ │ Base (μs) │ 0.0548293 │ 0.221979 │ 0.77535 │ 0.0877129 │
17+
│ 4 │ │ max error │ 0.0 │ 0.0 │ 0.0 │ 0.0 │
18+
│ 5 │ 2 │ compile time (s) │ 0.00551981 │ 0.00544899 │ 0.0109721 │ 0.109253 │
19+
│ 6 │ │ StaticArrays (μs) │ 0.001813 │ 0.002064 │ 0.0479443 │ 0.0349496 │
20+
│ 7 │ │ Base (μs) │ 0.192087 │ 0.554522 │ 2.60044 │ 1.4488 │
21+
│ 8 │ │ max error │ 2.22045e-16 │ 1.31167e-9 │ 3.17919e-15 │ 2.27485e-13 │
22+
│ 9 │ 3 │ compile time (s) │ 0.0810862 │ 0.0491976 │ 1.11676 │ 0.0138903 │
23+
│ 10 │ │ StaticArrays (μs) │ 0.001853 │ 0.00998799 │ 0.148526 │ 0.0606544 │
24+
│ 11 │ │ Base (μs) │ 0.244762 │ 0.671916 │ 3.19611 │ 1.8104 │
25+
│ 12 │ │ max error │ 2.22045e-16 │ 3.40212e-7 │ 5.55049e-15 │ 8.41688e-11 │
26+
│ 13 │ 4 │ compile time (s) │ 0.0180314 │ 0.846753 │ 1.10441 │ 0.0243759 │
27+
│ 14 │ │ StaticArrays (μs) │ 0.0104805 │ 0.0403582 │ 0.30312 │ 0.113047 │
28+
│ 15 │ │ Base (μs) │ 0.295459 │ 0.92675 │ 5.10971 │ 2.09289 │
29+
│ 16 │ │ max error │ 2.77556e-16 │ 2.62963e-10 │ 1.16735e-14 │ 2.01389e-9 │
30+
│ 17 │ 5 │ compile time (s) │ 2.92353 │ 0.41665 │ 1.58469 │ 0.0376695 │
31+
│ 18 │ │ StaticArrays (μs) │ 0.105874 │ 0.204571 │ 0.519827 │ 0.198071 │
32+
│ 19 │ │ Base (μs) │ 0.373642 │ 1.047 │ 5.89783 │ 2.2522 │
33+
│ 20 │ │ max error │ 1.66533e-16 │ 2.15962e-9 │ 1.84223e-14 │ 5.0253e-10 │
34+
│ 21 │ 6 │ compile time (s) │ 0.594495 │ 0.26922 │ 2.571 │ 0.0567599 │
35+
│ 22 │ │ StaticArrays (μs) │ 0.20246 │ 0.331432 │ 0.83373 │ 0.321775 │
36+
│ 23 │ │ Base (μs) │ 0.442449 │ 1.2143 │ 6.5424 │ 2.3995 │
37+
│ 24 │ │ max error │ 1.66533e-16 │ 1.53235e-8 │ 4.67587e-14 │ 6.4429e-9 │
38+
│ 25 │ 7 │ compile time (s) │ 0.679751 │ 0.367467 │ 4.41147 │ 0.0837663 │
39+
│ 26 │ │ StaticArrays (μs) │ 0.285484 │ 0.479979 │ 1.3115 │ 0.481318 │
40+
│ 27 │ │ Base (μs) │ 0.533771 │ 1.4007 │ 7.57425 │ 2.65178 │
41+
│ 28 │ │ max error │ 1.66533e-16 │ 1.54685e-7 │ 1.13318e-13 │ 3.0299e-9 │
42+
│ 29 │ 8 │ compile time (s) │ 0.798822 │ 0.501215 │ 7.42485 │ 0.125385 │
43+
│ 30 │ │ StaticArrays (μs) │ 0.400655 │ 0.672679 │ 1.7624 │ 0.678955 │
44+
│ 31 │ │ Base (μs) │ 0.613287 │ 1.603 │ 8.035 │ 2.94444 │
45+
│ 32 │ │ max error │ 2.77556e-16 │ 4.73129e-9 │ 1.46772e-13 │ 2.61726e-9 │
46+
│ 33 │ 9 │ compile time (s) │ 0.91446 │ 0.677859 │ 3.29461 │ 0.188982 │
47+
│ 34 │ │ StaticArrays (μs) │ 0.574361 │ 0.9838 │ 2.81078 │ 0.965474 │
48+
│ 35 │ │ Base (μs) │ 0.7649 │ 1.8765 │ 9.128 │ 3.17475 │
49+
│ 36 │ │ max error │ 3.46945e-16 │ 6.68251e-9 │ 3.30702e-13 │ 2.04139e-9 │
50+
│ 37 │ 10 │ compile time (s) │ 1.09413 │ 0.923797 │ 4.16234 │ 0.258198 │
51+
│ 38 │ │ StaticArrays (μs) │ 0.770716 │ 1.3145 │ 4.16643 │ 1.3064 │
52+
│ 39 │ │ Base (μs) │ 2.82311 │ 4.12057 │ 12.835 │ 5.3835 │
53+
│ 40 │ │ max error │ 1.94289e-16 │ 1.7691e-9 │ 2.42841e-13 │ 1.97784e-9 │
54+
│ 41 │ 11 │ compile time (s) │ 1.30448 │ 1.21905 │ 5.35341 │ 0.396076 │
55+
│ 42 │ │ StaticArrays (μs) │ 1.037 │ 1.7583 │ 5.43033 │ 1.7323 │
56+
│ 43 │ │ Base (μs) │ 3.24737 │ 4.82417 │ 16.481 │ 6.83967 │
57+
│ 44 │ │ max error │ 2.42861e-16 │ 1.91576e-9 │ 5.14422e-13 │ 3.63391e-9 │
58+
│ 45 │ 12 │ compile time (s) │ 1.55567 │ 1.66793 │ 6.75529 │ 0.601386 │
59+
│ 46 │ │ StaticArrays (μs) │ 1.3044 │ 2.203 │ 6.8148 │ 2.18522 │
60+
│ 47 │ │ Base (μs) │ 3.3915 │ 5.04783 │ 15.549 │ 6.2778 │
61+
│ 48 │ │ max error │ 2.77556e-16 │ 6.40505e-10 │ 7.41502e-13 │ 9.72575e-10 │
62+
│ 49 │ 13 │ compile time (s) │ 1.93222 │ 2.21769 │ 8.61503 │ 0.847173 │
63+
│ 50 │ │ StaticArrays (μs) │ 1.6922 │ 2.83978 │ 8.623 │ 2.812 │
64+
│ 51 │ │ Base (μs) │ 3.841 │ 5.813 │ 25.147 │ 6.2298 │
65+
│ 52 │ │ max error │ 2.22045e-16 │ 3.45741e-8 │ 1.43535e-12 │ 3.7065e-10 │
66+
│ 53 │ 14 │ compile time (s) │ 2.26075 │ 2.94048 │ 10.9119 │ 1.22455 │
67+
│ 54 │ │ StaticArrays (μs) │ 2.11289 │ 3.49162 │ 10.63 │ 3.43775 │
68+
│ 55 │ │ Base (μs) │ 6.0394 │ 8.1555 │ 29.706 │ 8.3256 │
69+
│ 56 │ │ max error │ 5.55112e-16 │ 2.75214e-8 │ 2.81287e-12 │ 1.23022e-8 │
70+
│ 57 │ 15 │ compile time (s) │ 0.0109467 │ 0.0650982 │ 16.0084 │ 0.059017 │
71+
│ 58 │ │ StaticArrays (μs) │ 6.788 │ 8.787 │ 37.942 │ 11.321 │
72+
│ 59 │ │ Base (μs) │ 6.6 │ 8.43933 │ 24.837 │ 9.819 │
73+
│ 60 │ │ max error │ 0.0 │ 0.0 │ 6.88512e-12 │ 0.0 │
74+
│ 61 │ 16 │ compile time (s) │ 0.0108788 │ 0.0883042 │ 18.8852 │ 0.0656792 │
75+
│ 62 │ │ StaticArrays (μs) │ 7.4815 │ 10.4963 │ 44.585 │ 11.091 │
76+
│ 63 │ │ Base (μs) │ 7.116 │ 9.638 │ 36.8 │ 12.223 │
77+
│ 64 │ │ max error │ 0.0 │ 0.0 │ 8.35981e-12 │ 0.0 │
78+
│ 65 │ 17 │ compile time (s) │ 0.0110253 │ 0.103676 │ 23.3029 │ 0.0761904 │
79+
│ 66 │ │ StaticArrays (μs) │ 8.656 │ 11.402 │ 53.001 │ 10.37 │
80+
│ 67 │ │ Base (μs) │ 8.386 │ 11.592 │ 47.86 │ 12.814 │
81+
│ 68 │ │ max error │ 0.0 │ 0.0 │ 2.26813e-11 │ 0.0 │
82+
│ 69 │ 18 │ compile time (s) │ 0.0128607 │ 0.0897007 │ 27.215 │ 0.0820536 │
83+
│ 70 │ │ StaticArrays (μs) │ 11.061 │ 13.706 │ 62.048 │ 17.152 │
84+
│ 71 │ │ Base (μs) │ 12.584 │ 14.096 │ 55.174 │ 14.878 │
85+
│ 72 │ │ max error │ 0.0 │ 0.0 │ 1.89276e-11 │ 0.0 │
86+
│ 73 │ 19 │ compile time (s) │ 0.0110386 │ 0.113333 │ 32.6636 │ 0.090609 │
87+
│ 74 │ │ StaticArrays (μs) │ 10.48 │ 14.387 │ 71.615 │ 17.143 │
88+
│ 75 │ │ Base (μs) │ 13.325 │ 13.636 │ 60.074 │ 12.513 │
89+
│ 76 │ │ max error │ 0.0 │ 0.0 │ 2.97618e-11 │ 0.0 │
90+
│ 77 │ 20 │ compile time (s) │ 0.0111067 │ 0.110741 │ 39.7384 │ 0.0973694 │
91+
│ 78 │ │ StaticArrays (μs) │ 11.2547 │ 14.768 │ 69.351 │ 13.365 │
92+
│ 79 │ │ Base (μs) │ 10.59 │ 13.887 │ 62.789 │ 17.813 │
93+
│ 80 │ │ max error │ 0.0 │ 0.0 │ 1.64377e-11 │ 0.0 │

‎perf/benchmark_matrix_ops.jl

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
using StaticArrays, BenchmarkTools, DataFrames, DataStructures
2+
3+
Nmax = 20
4+
unary = (det, inv, expm)
5+
binary = (\,)
6+
7+
data = OrderedDict{Symbol,Any}()
8+
data[:SIZE] = vcat(([i, "", "", ""] for i in 1:Nmax)...)
9+
data[:STAT] = [stat for sz in 1:Nmax for stat in ("compile time (s)", "StaticArrays (μs)", "Base (μs)", "max error")]
10+
11+
for f in unary
12+
f_data = Float64[]
13+
for N in 1:Nmax
14+
print("\r$((f,N))")
15+
SA = @SMatrix rand(N,N)
16+
A = Array(SA)
17+
push!(f_data, @elapsed f(SA))
18+
push!(f_data, 1e6*@belapsed $f($SA))
19+
push!(f_data, 1e6*@belapsed $f($A))
20+
push!(f_data, maximum([begin
21+
SA = @SMatrix rand(N,N)
22+
A = Array(SA)
23+
norm(f(A) - f(SA))
24+
end for i in 1:1000]))
25+
end
26+
data[Symbol(f)] = f_data
27+
end
28+
29+
for f in binary
30+
f_data = Float64[]
31+
for N in 1:Nmax
32+
print("\r$((f,N))")
33+
SA = @SMatrix rand(N,N)
34+
A = Array(SA)
35+
SB = @SMatrix rand(N,N)
36+
B = Array(SB)
37+
push!(f_data, @elapsed f(SA,SB))
38+
push!(f_data, 1e6*@belapsed $f($SA,$SB))
39+
push!(f_data, 1e6*@belapsed $f($A,$B))
40+
push!(f_data, maximum([begin
41+
SA = @SMatrix rand(N,N)
42+
A = Array(SA)
43+
SB = @SMatrix rand(N,N)
44+
B = Array(SB)
45+
norm(f(A,B) - f(SA,SB))
46+
end for i in 1:1000]))
47+
end
48+
data[Symbol(f)] = f_data
49+
end
50+
51+
df = DataFrame(data...)
52+
53+
open("bench_matrix_ops.txt", "w") do f
54+
versioninfo(f)
55+
showall(f, df)
56+
end

‎src/det.jl

Lines changed: 34 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -34,10 +34,39 @@ end
3434
A[5] * A[2] * A[11] * A[16] + A[1] * A[6] * A[11] * A[16])
3535
end
3636

37-
@inline function _det(::Size, A::StaticMatrix)
38-
return det(Matrix(A))
37+
@inline function _parity(p) # inefficient compared to computing cycle lengths, but non-allocating
38+
s = 0
39+
for i in 1:length(p), j in i+1:length(p)
40+
s += p[i] > p[j]
41+
end
42+
-2*rem(s, 2) + 1
3943
end
4044

41-
@inline logdet(a::StaticMatrix) = _logdet(Size(a), a)
42-
@inline _logdet(::Union{Size{(1,1)}, Size{(2,2)}, Size{(3,3)}}, a::StaticMatrix) = log(det(a))
43-
@inline _logdet(::Size, a::StaticMatrix) = logdet(drop_sdims(a))
45+
@generated function _det(::Size{S}, A::StaticMatrix) where S
46+
LinearAlgebra.checksquare(A)
47+
if prod(S) ≤ 14*14
48+
quote
49+
@_inline_meta
50+
L, U, p = lu(A)
51+
det(U)*_parity(p)
52+
end
53+
else
54+
:(@_inline_meta; det(Matrix(A)))
55+
end
56+
end
57+
58+
@inline logdet(A::StaticMatrix) = _logdet(Size(A), A)
59+
@inline _logdet(::Union{Size{(1,1)}, Size{(2,2)}, Size{(3,3)}, Size{(4,4)}}, A::StaticMatrix) = log(det(A))
60+
@generated function _logdet(::Size{S}, A::StaticMatrix) where S
61+
LinearAlgebra.checksquare(A)
62+
if prod(S) ≤ 14*14
63+
quote
64+
@_inline_meta
65+
L, U, p = lu(A)
66+
d, s = logabsdet(U)
67+
d + log(s*_parity(p))
68+
end
69+
else
70+
:(@_inline_meta; logdet(drop_sdims(A)))
71+
end
72+
end

‎src/expm.jl

Lines changed: 54 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -66,9 +66,58 @@ end
6666
(newtype)((m11, m21, m12, m22))
6767
end
6868

69-
# TODO add special case for 3x3 matrices
70-
if VERSION < v"0.7-"
71-
@inline _exp(s::Size, A::StaticArray) = s(Base.expm(Array(A)))
72-
else
73-
@inline _exp(s::Size, A::StaticArray) = s(Base.exp(Array(A)))
69+
# Adapted from implementation in Base; algorithm from
70+
# Higham, "Functions of Matrices: Theory and Computation", SIAM, 2008
71+
function _exp(::Size, A::StaticMatrix{<:Any,<:Any,T}) where T
72+
# omitted: matrix balancing, i.e., LAPACK.gebal!
73+
nA = maximum(sum(abs.(A), Val{1})) # marginally more performant than norm(A, 1)
74+
## For sufficiently small nA, use lower order Padé-Approximations
75+
if (nA <= 2.1)
76+
A2 = A*A
77+
if nA > 0.95
78+
U = @evalpoly(A2, T(8821612800)*I, T(302702400)*I, T(2162160)*I, T(3960)*I, T(1)*I)
79+
U = A*U
80+
V = @evalpoly(A2, T(17643225600)*I, T(2075673600)*I, T(30270240)*I, T(110880)*I, T(90)*I)
81+
elseif nA > 0.25
82+
U = @evalpoly(A2, T(8648640)*I, T(277200)*I, T(1512)*I, T(1)*I)
83+
U = A*U
84+
V = @evalpoly(A2, T(17297280)*I, T(1995840)*I, T(25200)*I, T(56)*I)
85+
elseif nA > 0.015
86+
U = @evalpoly(A2, T(15120)*I, T(420)*I, T(1)*I)
87+
U = A*U
88+
V = @evalpoly(A2, T(30240)*I, T(3360)*I, T(30)*I)
89+
else
90+
U = @evalpoly(A2, T(60)*I, T(1)*I)
91+
U = A*U
92+
V = @evalpoly(A2, T(120)*I, T(12)*I)
93+
end
94+
expA = (V - U) \ (V + U)
95+
else
96+
s = log2(nA/5.4) # power of 2 later reversed by squaring
97+
if s > 0
98+
si = ceil(Int,s)
99+
A = A / T(2^si)
100+
end
101+
102+
A2 = A*A
103+
A4 = A2*A2
104+
A6 = A2*A4
105+
106+
U = A6*(T(1)*A6 + T(16380)*A4 + T(40840800)*A2) +
107+
(T(33522128640)*A6 + T(10559470521600)*A4 + T(1187353796428800)*A2) +
108+
T(32382376266240000)*I
109+
U = A*U
110+
V = A6*(T(182)*A6 + T(960960)*A4 + T(1323241920)*A2) +
111+
(T(670442572800)*A6 + T(129060195264000)*A4 + T(7771770303897600)*A2) +
112+
T(64764752532480000)*I
113+
expA = (V - U) \ (V + U)
114+
115+
if s > 0 # squaring to reverse dividing by power of 2
116+
for t=1:si
117+
expA = expA*expA
118+
end
119+
end
120+
end
121+
122+
expA
74123
end

‎src/inv.jl

Lines changed: 11 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -55,6 +55,15 @@ end
5555
return similar_type(A)(B)
5656
end
5757

58-
@inline function _inv(::Size, A)
59-
similar_type(A)(inv(Matrix(A)))
58+
@generated function _inv(::Size{S}, A) where S
59+
LinearAlgebra.checksquare(A)
60+
if prod(S) ≤ 14*14
61+
quote
62+
@_inline_meta
63+
L, U, p = lu(A)
64+
U \ (L \ eye(A)[p,:])
65+
end
66+
else
67+
:(@_inline_meta; similar_type(A)(inv(Matrix(A))))
68+
end
6069
end

‎src/linalg.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -525,7 +525,7 @@ end
525525
@inline LinearAlgebra.checksquare(::SM) where {SM<:StaticMatrix} = _checksquare(Size(SM))
526526
@inline LinearAlgebra.checksquare(::Type{SM}) where {SM<:StaticMatrix} = _checksquare(Size(SM))
527527

528-
@pure _checksquare(::Size{S}) where {S} = (S[1] == S[2] || error("marix must be square"); S[1])
528+
@pure _checksquare(::Size{S}) where {S} = (S[1] == S[2] || throw(DimensionMismatch("matrix is not square: dimensions are $S")); S[1])
529529

530530
@inline LinearAlgebra.Symmetric(A::StaticMatrix, uplo::Char='U') = (LinearAlgebra.checksquare(A);Symmetric{eltype(A),typeof(A)}(A, uplo))
531531
@inline LinearAlgebra.Hermitian(A::StaticMatrix, uplo::Char='U') = (LinearAlgebra.checksquare(A);Hermitian{eltype(A),typeof(A)}(A, uplo))

‎src/solve.jl

Lines changed: 21 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,4 @@
1-
@inline (\)(a::StaticMatrix, b::StaticVector) = solve(Size(a), Size(b), a, b)
2-
3-
# TODO: Ineffecient but requires some infrastructure (e.g. LU or QR) to make efficient so we fall back on inv for now
4-
@inline solve(::Size, ::Size, a, b) = inv(a) * b
1+
@inline (\)(a::StaticMatrix, b::StaticVecOrMat) = solve(Size(a), Size(b), a, b)
52

63
@inline function solve(::Size{(1,1)}, ::Size{(1,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb}) where {Ta, Tb}
74
@inbounds return similar_type(b, typeof(a[1] \ b[1]))(a[1] \ b[1])
@@ -28,3 +25,23 @@ end
2825
(a[1,2]*a[3,1] - a[1,1]*a[3,2])*b[2] +
2926
(a[1,1]*a[2,2] - a[1,2]*a[2,1])*b[3]) / d )
3027
end
28+
29+
@generated function solve(::Size{Sa}, ::Size{Sb}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVecOrMat{Tb}) where {Sa, Sb, Ta, Tb}
30+
if Sa[end] != Sb[1]
31+
throw(DimensionMismatch("right hand side B needs first dimension of size $(Sa[end]), has size $Sb"))
32+
end
33+
LinearAlgebra.checksquare(a)
34+
if prod(Sa) ≤ 14*14
35+
quote
36+
@_inline_meta
37+
L, U, p = lu(a)
38+
U \ (L \ $(length(Sb) > 1 ? :(b[p,:]) : :(b[p])))
39+
end
40+
else
41+
quote
42+
@_inline_meta
43+
T = typeof((one(Ta)*zero(Tb) + one(Ta)*zero(Tb))/one(Ta))
44+
similar_type(b, T)(Matrix(a) \ b)
45+
end
46+
end
47+
end

‎test/det.jl

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,5 +24,17 @@
2424
@test det(Mtag) == det(Array(Mtag))
2525
end
2626

27+
# lu-based (sz in 5:14) and fallback (sz > 15)
28+
for sz in (5, 14, 15, 50), typ in (Float64, Complex{Float64})
29+
A = rand(typ, sz, sz)
30+
SA = SMatrix{sz,sz,typ}(A)
31+
@test det(A) ≈ det(SA)
32+
if typ == Float64 && det(A) < 0
33+
A[:,1], A[:,2] = A[:,2], A[:,1]
34+
SA = SMatrix{sz,sz,typ}(A)
35+
end
36+
@test logdet(A) ≈ logdet(SA)
37+
end
38+
2739
@test_throws DimensionMismatch det(@SMatrix [0; 1])
2840
end

‎test/expm.jl

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,5 +4,15 @@
44
@test expm(@SMatrix [4 2; -2 1])::SMatrix ≈ expm([4 2; -2 1])
55
@test expm(@SMatrix [4 2; 2 1])::SMatrix ≈ expm([4 2; 2 1])
66
@test expm(@SMatrix [ -3+1im -1+2im;-4-5im 5-2im])::SMatrix ≈ expm(Complex{Float64}[ -3+1im -1+2im;-4-5im 5-2im])
7-
@test expm(@SMatrix [1 2 0; 2 1 0; 0 0 1])::SizedArray{Tuple{3,3}} ≈ expm([1 2 0; 2 1 0; 0 0 1])
7+
@test expm(@SMatrix [1 2 0; 2 1 0; 0 0 1]) ≈ expm([1 2 0; 2 1 0; 0 0 1])
8+
9+
for sz in (3,4), typ in (Float64, Complex{Float64})
10+
A = rand(typ, sz, sz)
11+
nA = norm(A, 1)
12+
for nB in (0.005, 0.1, 0.5, 3.0, 20.0)
13+
B = A*nB/nA
14+
SB = SMatrix{sz,sz,typ}(B)
15+
@test expm(B) ≈ expm(SB)
16+
end
17+
end
818
end

‎test/inv.jl

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -90,6 +90,14 @@ end
9090
@test inv(SMatrix{5,5}(m))::StaticMatrix ≈ inv(m)
9191
end
9292

93+
@testset "Matrix inverse NxN, N > 5" begin
94+
for sz in (5, 14, 15, 50), typ in (Float64, Complex{Float64})
95+
A = rand(typ, sz, sz)
96+
SA = SMatrix{sz,sz,typ}(A)
97+
@test inv(A) ≈ inv(SA)
98+
end
99+
end
100+
93101
#-------------------------------------------------------------------------------
94102
# More comprehensive but qualitiative testing for inv() accuracy
95103
#=

0 commit comments

Comments
 (0)