Report index
StatProfilerHTML.jl report
Generated on Mon, 01 Apr 2024 21:01:18
File source code
Line Exclusive Inclusive Code
1 import Base: +, -, *, /, \
2
3 #--------------------------------------------------
4 # Vector space algebra
5
6 # Unary ops
7 @inline +(a::StaticArray) = map(+, a)
8 @inline -(a::StaticArray) = map(-, a)
9
10 # Binary ops
11 # Between arrays
12 @inline +(a::StaticArray, b::StaticArray) = map(+, a, b)
13 @inline +(a::AbstractArray, b::StaticArray) = map(+, a, b)
14 @inline +(a::StaticArray, b::AbstractArray) = map(+, a, b)
15
16 @inline -(a::StaticArray, b::StaticArray) = map(-, a, b)
17 @inline -(a::AbstractArray, b::StaticArray) = map(-, a, b)
18 @inline -(a::StaticArray, b::AbstractArray) = map(-, a, b)
19
20 # Scalar-array
21 @inline *(a::Number, b::StaticArray) = map(c->a*c, b)
22 @inline *(a::StaticArray, b::Number) = map(c->c*b, a)
23
24 @inline /(a::StaticArray, b::Number) = map(c->c/b, a)
25 @inline \(a::Number, b::StaticArray) = map(c->a\c, b)
26
27
28 # With UniformScaling
29 @inline +(a::StaticMatrix, b::UniformScaling) = _plus_uniform(Size(a), a, b.λ)
30 @inline +(a::UniformScaling, b::StaticMatrix) = _plus_uniform(Size(b), b, a.λ)
31 @inline -(a::StaticMatrix, b::UniformScaling) = _plus_uniform(Size(a), a, -b.λ)
32 @inline -(a::UniformScaling, b::StaticMatrix) = _plus_uniform(Size(b), -b, a.λ)
33
34 @generated function _plus_uniform(::Size{S}, a::StaticMatrix, λ) where {S}
35 n = checksquare(a)
36 exprs = [i == j ? :(a[$(LinearIndices(S)[i, j])] + λ) : :(a[$(LinearIndices(S)[i, j])]) for i = 1:n, j = 1:n]
37 return quote
38 $(Expr(:meta, :inline))
39 @inbounds return similar_type(a, promote_type(eltype(a), typeof(λ)))(tuple($(exprs...)))
40 end
41 end
42
43 @inline *(a::UniformScaling, b::Union{StaticMatrix,StaticVector}) = a.λ * b
44 @inline *(a::StaticMatrix, b::UniformScaling) = a * b.λ
45 @inline \(a::UniformScaling, b::Union{StaticMatrix,StaticVector}) = a.λ \ b
46 @inline /(a::StaticMatrix, b::UniformScaling) = a / b.λ
47
48
49 # Ternary ops
50 @inline Base.muladd(scalar::Number, a::StaticArray, b::StaticArray) = map((ai, bi) -> muladd(scalar, ai, bi), a, b)
51 @inline Base.muladd(a::StaticArray, scalar::Number, b::StaticArray) = map((ai, bi) -> muladd(ai, scalar, bi), a, b)
52
53
54 # @fastmath operators
55 @inline Base.FastMath.mul_fast(a::Number, b::StaticArray) = map(c -> Base.FastMath.mul_fast(a, c), b)
56 @inline Base.FastMath.mul_fast(a::StaticArray, b::Number) = map(c -> Base.FastMath.mul_fast(c, b), a)
57
58 @inline Base.FastMath.add_fast(a::StaticArray, b::StaticArray) = map(Base.FastMath.add_fast, a, b)
59 @inline Base.FastMath.sub_fast(a::StaticArray, b::StaticArray) = map(Base.FastMath.sub_fast, a, b)
60
61
62 #--------------------------------------------------
63 # Matrix algebra
64
65 # _adjointtype returns the eltype of the container when computing the adjoint/transpose
66 # of a static array. Using this method instead of calling `Base.promote_op` directly
67 # helps with type-inference, particularly for nested static arrays,
68 # where the adjoint is applied recursively.
69 @inline _adjointtype(f, ::Type{T}) where {T} = Base.promote_op(f, T)
70 for S in (:SMatrix, :MMatrix)
71 @eval @inline _adjointtype(f, ::Type{$S{M,N,T,L}}) where {M,N,T,L} = $S{N,M,_adjointtype(f, T),L}
72 end
73
74 # Transpose, etc
75 @inline transpose(m::StaticMatrix) = _transpose(Size(m), m)
76 # note: transpose of StaticVector is a Transpose, handled by Base
77 @inline transpose(a::Transpose{<:Any,<:Union{StaticVector,StaticMatrix}}) = a.parent
78 @inline transpose(a::Adjoint{<:Any,<:Union{StaticVector,StaticMatrix}}) = conj(a.parent)
79 @inline transpose(a::Adjoint{<:Real,<:Union{StaticVector,StaticMatrix}}) = a.parent
80
81 @generated function _transpose(::Size{S}, m::StaticMatrix{n1, n2, T}) where {n1, n2, S, T}
82 exprs = [:(transpose(m[$(LinearIndices(S)[j1, j2])])) for j2 in 1:n2, j1 in 1:n1]
83 return quote
84 $(Expr(:meta, :inline))
85 elements = tuple($(exprs...))
86 @inbounds return similar_type($m, _adjointtype(transpose, T), Size($(n2,n1)))(elements)
87 end
88 end
89
90 @inline adjoint(m::StaticMatrix) = _adjoint(Size(m), m)
91 @inline adjoint(a::Transpose{<:Any,<:Union{StaticVector,StaticMatrix}}) = conj(a.parent)
92 @inline adjoint(a::Transpose{<:Real,<:Union{StaticVector,StaticMatrix}}) = a.parent
93 @inline adjoint(a::Adjoint{<:Any,<:Union{StaticVector,StaticMatrix}}) = a.parent
94
95 @generated function _adjoint(::Size{S}, m::StaticMatrix{n1, n2, T}) where {n1, n2, S, T}
96 exprs = [:(adjoint(m[$(LinearIndices(S)[j1, j2])])) for j2 in 1:n2, j1 in 1:n1]
97 return quote
98 $(Expr(:meta, :inline))
99 elements = tuple($(exprs...))
100 @inbounds return similar_type($m, _adjointtype(adjoint, T), Size($(n2,n1)))(elements)
101 end
102 end
103
104 @inline Base.zero(a::SA) where {SA <: StaticArray} = zeros(SA)
105 @inline Base.zero(a::Type{SA}) where {SA <: StaticArray} = zeros(SA)
106
107 @inline _construct_sametype(a::Type, elements) = a(elements)
108 @inline _construct_sametype(a, elements) = typeof(a)(elements)
109
110 @inline one(m::StaticMatrixLike) = _one(Size(m), m)
111 @inline one(::Type{SM}) where {SM<:StaticMatrixLike}= _one(Size(SM), SM)
112 function _one(s::Size, m_or_SM)
113 if (length(s) != 2) || (s[1] != s[2])
114 throw(DimensionMismatch("multiplicative identity defined only for square matrices"))
115 end
116 λ = one(_eltype_or(m_or_SM, Float64))
117 _construct_sametype(m_or_SM, _scalar_matrix_elements(s, λ))
118 # TODO: Bring back the use of _construct_similar here:
119 # _construct_similar(m_or_SM, s, _scalar_matrix_elements(s, λ))
120 #
121 # However, because _construct_similar uses `similar_type`, it will be
122 # breaking for some StaticMatrix types (in particular Rotations.RotMatrix)
123 # which must have similar_type return a general type able to hold any
124 # matrix in the full general linear group.
125 #
126 # (Generally we're stuck here and things like RotMatrix really need to
127 # override one() and zero() themselves: on the one hand, one(RotMatrix)
128 # should return a RotMatrix, but zero(RotMatrix) can not be a RotMatrix.
129 # The best StaticArrays can do is to use similar_type to return an SMatrix
130 # for both of these, and let the more specialized library define the
131 # correct algebraic properties.)
132 end
133
134 # StaticMatrix(I::UniformScaling)
135 (::Type{SM})(I::UniformScaling) where {SM<:StaticMatrix} =
136 SM(_scalar_matrix_elements(Size(SM), I.λ))
137 # The following oddity is needed if we want `SArray{Tuple{2,3}}(I)` to work
138 # because we do not have `SArray{Tuple{2,3}} <: StaticMatrix`.
139 (::Type{SM})(I::UniformScaling) where {SM<:(StaticArray{Tuple{N,M}} where {N,M})} =
140 SM(_scalar_matrix_elements(Size(SM), I.λ))
141
142 # Construct a matrix with the scalar λ on the diagonal and zeros off the
143 # diagonal. The matrix can be non-square.
144 @generated function _scalar_matrix_elements(s::Size{S}, λ) where {S}
145 elements = Symbol[i == j ? :λ : :λzero for i in 1:S[1], j in 1:S[2]]
146 return quote
147 $(Expr(:meta, :inline))
148 λzero = zero(λ)
149 tuple($(elements...))
150 end
151 end
152
153 @generated function diagm(kv1::Pair{<:Val,<:StaticVector}, other_kvs::Pair{<:Val,<:StaticVector}...)
154 kvs = (kv1, other_kvs...)
155 diag_ind_and_length = [(kv.parameters[1].parameters[1], length(kv.parameters[2])) for kv in kvs]
156 N = maximum(abs(di) + dl for (di,dl) in diag_ind_and_length)
157 vs = [Symbol("v$i") for i=1:length(kvs)]
158 vs_exprs = [:(@inbounds $(vs[i]) = kvs[$i].second) for i=eachindex(kvs)]
159 element_exprs = Any[false for _=1:N*N]
160 for (i, (di, dl)) in enumerate(diag_ind_and_length)
161 diaginds = diagind(N, N, di)
162 for n = 1:dl
163 element_exprs[diaginds[n]] = :($(vs_exprs[i])[$n])
164 end
165 end
166 return quote
167 $(Expr(:meta, :inline))
168 kvs = (kv1, other_kvs...)
169 $(vs_exprs...)
170 @inbounds elements = tuple($(element_exprs...))
171 T = promote_tuple_eltype(elements)
172 @inbounds return similar_type(v1, T, Size($N,$N))(elements)
173 end
174 end
175 @inline diagm(v::StaticVector) = diagm(Val(0)=>v)
176
177 @inline diag(m::StaticMatrix, k::Type{Val{D}}=Val{0}) where {D} = _diag(Size(m), m, k)
178 @generated function _diag(::Size{S}, m::StaticMatrix, ::Type{Val{D}}) where {S,D}
179 S1, S2 = S
180 rng = D ≤ 0 ? range(1-D, step=S1+1, length=min(S1+D, S2)) :
181 range(D*S1+1, step=S1+1, length=min(S1, S2-D))
182 Snew = length(rng)
183 T = eltype(m)
184 exprs = [:(m[$i]) for i = rng]
185 return quote
186 $(Expr(:meta, :inline))
187 @inbounds return similar_type($m, Size($Snew))(tuple($(exprs...)))
188 end
189 end
190
191 #--------------------------------------------------
192 # Vector products
193 @inline cross(a::StaticVector, b::StaticVector) = _cross(same_size(a, b), a, b)
194 _cross(::Size{S}, a::StaticVector, b::StaticVector) where {S} = error("Cross product not defined for $(S[1])-vectors")
195 @inline function _cross(::Size{(2,)}, a::StaticVector, b::StaticVector)
196 @inbounds return a[1]*b[2] - a[2]*b[1]
197 end
198 @inline function _cross(::Size{(3,)}, a::StaticVector, b::StaticVector)
199 @inbounds return similar_type(a, typeof(a[2]*b[3]-a[3]*b[2]))((a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1]))
200 end
201 @inline function _cross(::Size{(2,)}, a::StaticVector{<:Any, <:Unsigned}, b::StaticVector{<:Any, <:Unsigned})
202 @inbounds return Signed(a[1]*b[2]) - Signed(a[2]*b[1])
203 end
204 @inline function _cross(::Size{(3,)}, a::StaticVector{<:Any, <:Unsigned}, b::StaticVector{<:Any, <:Unsigned})
205 @inbounds return similar_type(a, typeof(Signed(a[2]*b[3])-Signed(a[3]*b[2])))(((Signed(a[2]*b[3])-Signed(a[3]*b[2]), Signed(a[3]*b[1])-Signed(a[1]*b[3]), Signed(a[1]*b[2])-Signed(a[2]*b[1]))))
206 end
207
208 1 (2 %)
1 (2 %) samples spent in dot
1 (100 %) (incl.) when called from calc_particle_forces! line 253
1 (100 %) samples spent calling _vecdot
@inline dot(a::StaticArray, b::StaticArray) = _vecdot(same_size(a, b), a, b, dot)
209 @inline bilinear_vecdot(a::StaticArray, b::StaticArray) = _vecdot(same_size(a, b), a, b, *)
210
211 @inline function _vecdot(::Size{S}, a::StaticArray, b::StaticArray, product) where {S}
212 if prod(S) == 0
213 za, zb = zero(eltype(a)), zero(eltype(b))
214 else
215 # Use an actual element if there is one, to support e.g. Vector{<:Number}
216 # element types for which runtime size information is required to construct
217 # a zero element.
218 za, zb = zero(a[1]), zero(b[1])
219 end
220 ret = product(za, zb) + product(za, zb)
221 1 (2 %)
1 (2 %) samples spent in _vecdot
1 (100 %) (incl.) when called from dot line 208
1 (100 %) samples spent calling macro expansion
@inbounds @simd for j = 1 : prod(S)
222 1 (2 %)
1 (2 %) samples spent in macro expansion
1 (100 %) (incl.) when called from macro expansion line 77
1 (100 %) samples spent calling getindex
ret += product(a[j], b[j])
223 end
224 return ret
225 end
226
227 #--------------------------------------------------
228 # Norms
229 _inner_eltype(v::AbstractArray) = isempty(v) ? eltype(v) : _inner_eltype(first(v))
230 _inner_eltype(x::Number) = typeof(x)
231 @inline _init_zero(v::AbstractArray) = float(norm(zero(_inner_eltype(v))))
232
233 @inline function LinearAlgebra.norm_sqr(v::StaticArray)
234 return mapreduce(LinearAlgebra.norm_sqr, +, v; init=_init_zero(v))
235 end
236
237 @inline maxabs_nested(a::Number) = abs(a)
238 @inline function maxabs_nested(a::AbstractArray)
239 prod(size(a)) == 0 && (return _init_zero(a))
240
241 m = maxabs_nested(a[1])
242 for j = 2:prod(size(a))
243 m = max(m, maxabs_nested(a[j]))
244 end
245
246 return m
247 end
248
249 @generated function _norm_scaled(::Size{S}, a::StaticArray) where {S}
250 expr = :(LinearAlgebra.norm_sqr(a[1]/scale))
251 for j = 2:prod(S)
252 expr = :($expr + LinearAlgebra.norm_sqr(a[$j]/scale))
253 end
254
255 return quote
256 $(Expr(:meta, :inline))
257 scale = maxabs_nested(a)
258 !isfinite(scale) && return scale
259
260 iszero(scale) && return _init_zero(a)
261 return @inbounds scale * sqrt($expr)
262 end
263 end
264
265 @inline norm(a::StaticArray) = _norm(Size(a), a)
266 @generated function _norm(::Size{S}, a::StaticArray) where {S}
267 prod(S) == 0 && return :(_init_zero(a))
268
269 expr = :(LinearAlgebra.norm_sqr(a[1]))
270 for j = 2:prod(S)
271 expr = :($expr + LinearAlgebra.norm_sqr(a[$j]))
272 end
273
274 return quote
275 $(Expr(:meta, :inline))
276 l = @inbounds sqrt($expr)
277
278 zero(l) < l && isfinite(l) && return l
279 return _norm_scaled(Size(a), a)
280 end
281 end
282
283 function _norm_p0(x)
284 T = _inner_eltype(x)
285 return float(norm(iszero(x) ? zero(T) : one(T)))
286 end
287
288 # Do not need to deal with p == 0, 2, Inf; see norm(a, p).
289 @generated function _norm_scaled(::Size{S}, a::StaticArray, p::Real) where {S}
290 expr = :(norm(a[1]/scale)^p)
291 for j = 2:prod(S)
292 expr = :($expr + norm(a[$j]/scale)^p)
293 end
294
295 expr_p1 = :(norm(a[1]/scale))
296 for j = 2:prod(S)
297 expr_p1 = :($expr_p1 + norm(a[$j]/scale))
298 end
299
300 return quote
301 $(Expr(:meta, :inline))
302 scale = maxabs_nested(a)
303
304 scale==0 && return _init_zero(a)
305 p == 1 && return @inbounds scale * $expr_p1
306 return @inbounds scale * ($expr)^(inv(p))
307 end
308 end
309
310 @inline norm(a::StaticArray, p::Real) = _norm(Size(a), a, p)
311 @generated function _norm(::Size{S}, a::StaticArray, p::Real) where {S}
312 prod(S) == 0 && return :(_init_zero(a))
313
314 expr = :(norm(a[1])^p)
315 for j = 2:prod(S)
316 expr = :($expr + norm(a[$j])^p)
317 end
318
319 expr_p1 = :(norm(a[1]))
320 for j = 2:prod(S)
321 expr_p1 = :($expr_p1 + norm(a[$j]))
322 end
323
324 return quote
325 $(Expr(:meta, :inline))
326 p == 0 && return mapreduce(_norm_p0, +, a) # no need for scaling
327 p == 2 && return norm(a) # norm(a) takes care of scaling
328 p == Inf && return mapreduce(norm, max, a) # no need for scaling
329
330 l = p==1 ? @inbounds($expr_p1) : @inbounds(($expr)^(inv(p)))
331 0<l<Inf && return l
332 return _norm_scaled(Size(a), a, p) # p != 0, 2, Inf
333 end
334 end
335
336 @inline normalize(a::StaticArray) = inv(norm(a))*a
337 @inline normalize(a::StaticArray, p::Real) = inv(norm(a, p))*a
338
339 @inline normalize!(a::StaticArray) = (a .*= inv(norm(a)); return a)
340 @inline normalize!(a::StaticArray, p::Real) = (a .*= inv(norm(a, p)); return a)
341
342 @inline tr(a::StaticMatrix) = _tr(Size(a), a)
343 @generated function _tr(::Size{S}, a::StaticMatrix) where {S}
344 checksquare(a)
345
346 if S[1] == 0
347 return :(zero(eltype(a)))
348 end
349
350 exprs = [:(a[$(LinearIndices(S)[i, i])]) for i = 1:S[1]]
351 total = reduce((ex1, ex2) -> :($ex1 + $ex2), exprs)
352
353 return quote
354 @_inline_meta
355 @inbounds return $total
356 end
357 end
358
359
360 #--------------------------------------------------
361 # Outer products
362
363 const _length_limit = Length(200)
364
365 @inline kron(a::StaticMatrix, b::StaticMatrix) = _kron(_length_limit, Size(a), Size(b), a, b)
366 @generated function _kron(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
367 outsize = SA .* SB
368 if prod(outsize) > length_limit
369 return :( SizedMatrix{$(outsize[1]),$(outsize[2])}( kron(drop_sdims(a), drop_sdims(b)) ) )
370 end
371
372 M = [ :(a[$ia,$ja] * b[$ib,$jb]) for ib in 1:SB[1], ia in 1:SA[1], jb in 1:SB[2], ja in 1:SA[2] ]
373
374 return quote
375 @_inline_meta
376 @inbounds return similar_type($a, promote_type(eltype(a),eltype(b)), Size($(outsize)))(tuple($(M...)))
377 end
378 end
379
380 @inline kron(a::StaticVector, b::StaticVector) = _kron_vec_x_vec(_length_limit, Size(a), Size(b), a, b)
381 @generated function _kron_vec_x_vec(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
382 outsize = SA .* SB
383 if prod(outsize) > length_limit
384 return :( SizedVector{$(outsize[1])}( kron(drop_sdims(a), drop_sdims(b)) ) )
385 end
386
387 m = [ :(a[$ia] * b[$ib]) for ib in 1:SB[1], ia in 1:SA[1]]
388
389 return quote
390 @_inline_meta
391 @inbounds return similar_type($a, promote_type(eltype(a),eltype(b)), Size($(outsize)))(tuple($(m...)))
392 end
393 end
394
395 @inline function kron(
396 a::Union{Transpose{<:Number,<:StaticVector}, Adjoint{<:Number,<:StaticVector}},
397 b::Union{Transpose{<:Number,<:StaticVector}, Adjoint{<:Number,<:StaticVector}})
398 _kron_tvec_x_tvec(_length_limit, Size(a), Size(b), a, b)
399 end
400 @generated function _kron_tvec_x_tvec(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
401 outsize = SA .* SB
402 if prod(outsize) > length_limit
403 return :( SizedMatrix{$(outsize[1]),$(outsize[2])}( kron(drop_sdims(a), drop_sdims(b)) ) )
404 end
405
406 m = [ :(a[$ia] * b[$ib]) for ib in 1:SB[2], ia in 1:SA[2]]
407
408 return quote
409 @_inline_meta
410 @inbounds return similar_type($a, promote_type(eltype(a),eltype(b)), Size($(outsize)))(tuple($(m...)))
411 end
412 end
413
414 @inline function kron(
415 a::Union{Transpose{<:Number,<:StaticVector}, Adjoint{<:Number,<:StaticVector}},
416 b::StaticVector)
417 _kron_tvec_x_vec(_length_limit, Size(a), Size(b), a, b)
418 end
419 @generated function _kron_tvec_x_vec(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
420 outsize = (SA[1] * SB[1], SA[2])
421 if prod(outsize) > length_limit
422 return :( SizedMatrix{$(outsize[1]),$(outsize[2])}( kron(drop_sdims(a), drop_sdims(b)) ) )
423 end
424
425 M = [ :(a[$ia] * b[$ib]) for ib in 1:SB[1], ia in 1:SA[2]]
426
427 return quote
428 @_inline_meta
429 @inbounds return similar_type($a, promote_type(eltype(a),eltype(b)), Size($(outsize)))(tuple($(M...)))
430 end
431 end
432
433 @inline function kron(
434 a::StaticVector,
435 b::Union{Transpose{<:Number,<:StaticVector}, Adjoint{<:Number,<:StaticVector}})
436 _kron_vec_x_tvec(_length_limit, Size(a), Size(b), a, b)
437 end
438 @generated function _kron_vec_x_tvec(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
439 outsize = (SA[1] * SB[1], SB[2])
440 if prod(outsize) > length_limit
441 return :( SizedMatrix{$(outsize[1]),$(outsize[2])}( kron(drop_sdims(a), drop_sdims(b)) ) )
442 end
443
444 M = [ :(a[$ia] * b[$ib]) for ia in 1:SA[1], ib in 1:SB[2]]
445
446 return quote
447 @_inline_meta
448 @inbounds return similar_type($a, promote_type(eltype(a),eltype(b)), Size($(outsize)))(tuple($(M...)))
449 end
450 end
451
452 @inline kron(a::StaticVector, b::StaticMatrix) = _kron_vec_x_mat(_length_limit, Size(a), Size(b), a, b)
453 @generated function _kron_vec_x_mat(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
454 outsize = (SA[1] * SB[1], SB[2])
455 if prod(outsize) > length_limit
456 return :( SizedMatrix{$(outsize[1]),$(outsize[2])}( kron(drop_sdims(a), drop_sdims(b)) ) )
457 end
458
459 M = [ :(a[$ia] * b[$ib,$jb]) for ib in 1:SB[1], ia in 1:SA[1], jb in 1:SB[2]]
460
461 return quote
462 @_inline_meta
463 @inbounds return similar_type($b, promote_type(eltype(a),eltype(b)), Size($(outsize)))(tuple($(M...)))
464 end
465 end
466
467 @inline kron(a::StaticMatrix, b::StaticVector) = _kron_mat_x_vec(_length_limit, Size(a), Size(b), a, b)
468 @generated function _kron_mat_x_vec(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
469 outsize = (SA[1] * SB[1], SA[2])
470 if prod(outsize) > length_limit
471 return :( SizedMatrix{$(outsize[1]),$(outsize[2])}( kron(drop_sdims(a), drop_sdims(b)) ) )
472 end
473
474 M = [ :(a[$ia,$ja] * b[$ib]) for ib in 1:SB[1], ia in 1:SA[1], ja in 1:SA[2] ]
475
476 return quote
477 @_inline_meta
478 @inbounds return similar_type($a, promote_type(eltype(a),eltype(b)), Size($(outsize)))(tuple($(M...)))
479 end
480 end
481
482 @inline function kron(
483 a::StaticMatrix,
484 b::Union{Transpose{<:Number,<:StaticVector}, Adjoint{<:Number,<:StaticVector}})
485 _kron_mat_x_tvec(_length_limit, Size(a), Size(b), a, b)
486 end
487 @generated function _kron_mat_x_tvec(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
488 outsize = SA .* SB
489 if prod(outsize) > length_limit
490 return :( SizedMatrix{$(outsize[1]),$(outsize[2])}( kron(drop_sdims(a), drop_sdims(b)) ) )
491 end
492
493 M = [ :(a[$ia,$ja] * b[$ib,$jb]) for ib in 1:SB[1], ia in 1:SA[1], jb in 1:SB[2], ja in 1:SA[2] ]
494
495 return quote
496 @_inline_meta
497 @inbounds return similar_type($a, promote_type(eltype(a),eltype(b)), Size($(outsize)))(tuple($(M...)))
498 end
499 end
500
501 @inline function kron(
502 a::Union{Transpose{<:Number,<:StaticVector}, Adjoint{<:Number,<:StaticVector}},
503 b::StaticMatrix)
504 _kron_tvec_x_mat(_length_limit, Size(a), Size(b), a, b)
505 end
506 @generated function _kron_tvec_x_mat(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
507 outsize = SA .* SB
508 if prod(outsize) > length_limit
509 return :( SizedMatrix{$(outsize[1]),$(outsize[2])}( kron(drop_sdims(a), drop_sdims(b)) ) )
510 end
511
512 M = [ :(a[$ia,$ja] * b[$ib,$jb]) for ib in 1:SB[1], ia in 1:SA[1], jb in 1:SB[2], ja in 1:SA[2] ]
513
514 return quote
515 @_inline_meta
516 @inbounds return similar_type($b, promote_type(eltype(a),eltype(b)), Size($(outsize)))(tuple($(M...)))
517 end
518 end
519
520
521 #--------------------------------------------------
522 # Some shimming for special linear algebra matrix types
523 @inline LinearAlgebra.Symmetric(A::StaticMatrix, uplo::Char='U') = (checksquare(A); Symmetric{eltype(A),typeof(A)}(A, uplo))
524 @inline LinearAlgebra.Hermitian(A::StaticMatrix, uplo::Char='U') = (checksquare(A); Hermitian{eltype(A),typeof(A)}(A, uplo))