diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 5392dd1f42..2cf1a7cdca 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -180,6 +180,10 @@ steps: key: unit_rmul_with_projection command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/rmul_with_projection.jl" + - label: "Unit: simple_symmetric" + key: unit_simple_symmetric + command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/unit_simple_symmetric.jl" + - group: "Unit: Meshes" steps: diff --git a/src/Geometry/Geometry.jl b/src/Geometry/Geometry.jl index 030c4551a4..c6dacba617 100644 --- a/src/Geometry/Geometry.jl +++ b/src/Geometry/Geometry.jl @@ -22,8 +22,7 @@ export Contravariant1Vector, Contravariant23Vector, Contravariant123Vector - - +include("simple_symmetric.jl") include("coordinates.jl") include("axistensors.jl") include("localgeometry.jl") diff --git a/src/Geometry/axistensors.jl b/src/Geometry/axistensors.jl index 534e628380..a7bae11023 100644 --- a/src/Geometry/axistensors.jl +++ b/src/Geometry/axistensors.jl @@ -136,7 +136,7 @@ struct AxisTensor{ T, N, A <: NTuple{N, AbstractAxis}, - S <: StaticArray{<:Tuple, T, N}, + S <: Union{SimpleSymmetric{N, T}, StaticArray{<:Tuple, T, N}}, } <: AbstractArray{T, N} axes::A components::S @@ -147,7 +147,7 @@ AxisTensor( components::S, ) where { A <: Tuple{Vararg{AbstractAxis}}, - S <: StaticArray{<:Tuple, T, N}, + S <: Union{SimpleSymmetric{N, T}, StaticArray{<:Tuple, T, N}}, } where {T, N} = AxisTensor{T, N, A, S}(axes, components) AxisTensor(axes::Tuple{Vararg{AbstractAxis}}, components) = @@ -396,8 +396,16 @@ end function Base.:*(x::AxisVector, y::AdjointAxisVector) AxisTensor((axes(x, 1), axes(y, 2)), components(x) * components(y)) end +import InteractiveUtils function Base.:*(A::Axis2TensorOrAdj, x::AxisVector) check_dual(axes(A, 2), axes(x, 1)) + c = components(A) * components(x) + # @show typeof(A) + # @show typeof(x) + # s = InteractiveUtils.@which components(A) * components(x) + # @show s + # # @show typeof(x) + # @show c return AxisVector(axes(A, 1), components(A) * components(x)) end function Base.:*(A::Axis2TensorOrAdj, B::Axis2TensorOrAdj) diff --git a/src/Geometry/localgeometry.jl b/src/Geometry/localgeometry.jl index 6a8c6671ed..b5ccaaad34 100644 --- a/src/Geometry/localgeometry.jl +++ b/src/Geometry/localgeometry.jl @@ -1,14 +1,24 @@ +import LinearAlgebra +import .Geometry: SimpleSymmetric + +import InteractiveUtils import LinearAlgebra: issymmetric isapproxsymmetric(A::AbstractMatrix{T}; rtol = 10 * eps(T)) where {T} = Base.isapprox(A, A'; rtol) +function SimpleSymmetric(x::AxisTensor{FT}) where {FT} + c = SimpleSymmetric(components(x)) + A = axes(x) + return AxisTensor{FT, ndims(x), typeof(A), typeof(c)}(A, c) +end + """ LocalGeometry The necessary local metric information defined at each node. """ -struct LocalGeometry{I, C <: AbstractPoint, FT, S} +struct LocalGeometry{I, C <: AbstractPoint, FT, S, N, L} "Coordinates of the current point" coordinates::C "Jacobian determinant of the transformation `ξ` to `x`" @@ -22,9 +32,17 @@ struct LocalGeometry{I, C <: AbstractPoint, FT, S} "Partial derivatives of the map from `x` to `ξ`: `∂ξ∂x[i,j]` is ∂ξⁱ/∂xʲ" ∂ξ∂x::Axis2Tensor{FT, Tuple{ContravariantAxis{I}, LocalAxis{I}}, S} "Contravariant metric tensor (inverse of gᵢⱼ), transforms covariant to contravariant vector components" - gⁱʲ::Axis2Tensor{FT, Tuple{ContravariantAxis{I}, ContravariantAxis{I}}, S} + gⁱʲ::Axis2Tensor{ + FT, + Tuple{ContravariantAxis{I}, ContravariantAxis{I}}, + SimpleSymmetric{N, FT, L}, + } "Covariant metric tensor (gᵢⱼ), transforms contravariant to covariant vector components" - gᵢⱼ::Axis2Tensor{FT, Tuple{CovariantAxis{I}, CovariantAxis{I}}, S} + gᵢⱼ::Axis2Tensor{ + FT, + Tuple{CovariantAxis{I}, CovariantAxis{I}}, + SimpleSymmetric{N, FT, L}, + } @inline function LocalGeometry( coordinates, J, @@ -34,15 +52,27 @@ struct LocalGeometry{I, C <: AbstractPoint, FT, S} ∂ξ∂x = inv(∂x∂ξ) C = typeof(coordinates) Jinv = inv(J) - gⁱʲ = ∂ξ∂x * ∂ξ∂x' - gᵢⱼ = ∂x∂ξ' * ∂x∂ξ - isapproxsymmetric(components(gⁱʲ)) || error("gⁱʲ is not symmetric.") - isapproxsymmetric(components(gᵢⱼ)) || error("gᵢⱼ is not symmetric.") - return new{I, C, FT, S}(coordinates, J, WJ, Jinv, ∂x∂ξ, ∂ξ∂x, gⁱʲ, gᵢⱼ) + gⁱʲ₀ = ∂ξ∂x * ∂ξ∂x' + gᵢⱼ₀ = ∂x∂ξ' * ∂x∂ξ + isapproxsymmetric(components(gⁱʲ₀)) || error("gⁱʲ is not symmetric.") + isapproxsymmetric(components(gᵢⱼ₀)) || error("gᵢⱼ is not symmetric.") + gⁱʲ = SimpleSymmetric(gⁱʲ₀) + gᵢⱼ = SimpleSymmetric(gᵢⱼ₀) + L = triangular_nonzeros(S) + N = size(components(gⁱʲ₀), 1) + return new{I, C, FT, S, N, L}( + coordinates, + J, + WJ, + Jinv, + ∂x∂ξ, + ∂ξ∂x, + gⁱʲ, + gᵢⱼ, + ) end end - """ SurfaceGeometry @@ -55,7 +85,7 @@ struct SurfaceGeometry{FT, N} normal::N end -undertype(::Type{LocalGeometry{I, C, FT, S}}) where {I, C, FT, S} = FT +undertype(::Type{<:LocalGeometry{I, C, FT}}) where {I, C, FT} = FT undertype(::Type{SurfaceGeometry{FT, N}}) where {FT, N} = FT """ diff --git a/src/Geometry/simple_symmetric.jl b/src/Geometry/simple_symmetric.jl new file mode 100644 index 0000000000..e654d8f188 --- /dev/null +++ b/src/Geometry/simple_symmetric.jl @@ -0,0 +1,326 @@ +#= +This was adapted from StaticArrays' `SHermitianCompact`: +https://github.com/JuliaArrays/StaticArrays.jl/blob/609aa343a0504c3f5bae82a6236c563a0fa72681/src/SHermitianCompact.jl +=# + +import LinearAlgebra +import StaticArrays: SVector, SMatrix + +""" + SimpleSymmetric{N, T, L} <: StaticMatrix{N, N, T} + +A [`StaticArray`](@ref) subtype that represents a symmetric matrix. Unlike +`LinearAlgebra.Symmetric`, `SimpleSymmetric` stores only the lower triangle +of the matrix (as an `SVector`). The lower triangle is stored in column-major order. +For example, for an `SimpleSymmetric{3}`, the indices of the stored elements +can be visualized as follows: + +``` +┌ 1 ⋅ ⋅ ┐ +| 2 4 ⋅ | +└ 3 5 6 ┘ +``` + +Type parameters: +* `N`: matrix dimension; +* `T`: element type for lower triangle; +* `L`: length of the `SVector` storing the lower triangular elements. + +Note that `L` is always the `N`th [triangular number](https://en.wikipedia.org/wiki/Triangular_number). + +A `SimpleSymmetric` may be constructed either: + +* from an `AbstractVector` containing the lower triangular elements; or +* from a `Tuple` containing both upper and lower triangular elements in column major order; or +* from another `StaticMatrix`. + +For the latter two cases, only the lower triangular elements are used; the upper triangular +elements are ignored. +""" +struct SimpleSymmetric{N, T, L} <: StaticMatrix{N, N, T} + lowertriangle::SVector{L, T} + + @inline function SimpleSymmetric{N, T, L}( + lowertriangle::SVector{L}, + ) where {N, T, L} + _check_simple_symmetric_parameters(Val(N), Val(L)) + new{N, T, L}(lowertriangle) + end +end +StaticArrays.check_parameters( + ::Type{SimpleSymmetric{N, T, L}}, +) where {N, T, L} = _check_simple_symmetric_parameters(Val(N), Val(L)) + +triangular_nonzeros(::SMatrix{N}) where {N} = Int(N * (N + 1) / 2) +triangular_nonzeros(::Type{<:SMatrix{N}}) where {N} = Int(N * (N + 1) / 2) +tail_params(::Type{S}) where {N, T, S <: SMatrix{N, N, T}} = + (T, S, N, triangular_nonzeros(S)) + +@generated function SimpleSymmetric(A::S) where {S <: SMatrix} + N = size(S, 1) + L = triangular_nonzeros(S) + _check_simple_symmetric_parameters(Val(N), Val(L)) + expr = Vector{Expr}(undef, L) + T = eltype(S) + i = 0 + for col in 1:N, row in 1:N + if col ≥ row + expr[i += 1] = :(A[$row, $col]) + end + end + quote + Base.@_inline_meta + @inbounds return SimpleSymmetric{$N, $T, $L}( + SVector{$L, $T}(tuple($(expr...))), + ) + end +end + +@inline function _check_simple_symmetric_parameters( + ::Val{N}, + ::Val{L}, +) where {N, L} + if 2 * L !== N * (N + 1) + throw( + ArgumentError( + "Size mismatch in SimpleSymmetric parameters. Got dimension $N and length $L.", + ), + ) + end +end + +triangularnumber(N::Int) = div(N * (N + 1), 2) +@generated function triangularroot(::Val{L}) where {L} + return div(isqrt(8 * L + 1) - 1, 2) # from quadratic formula +end + +lowertriangletype(::Type{SimpleSymmetric{N, T, L}}) where {N, T, L} = + SVector{L, T} +lowertriangletype(::Type{SimpleSymmetric{N, T}}) where {N, T} = + SVector{triangularnumber(N), T} +lowertriangletype(::Type{SimpleSymmetric{N}}) where {N} = + SVector{triangularnumber(N)} + +@inline SimpleSymmetric{N, T}(lowertriangle::SVector{L}) where {N, T, L} = + SimpleSymmetric{N, T, L}(lowertriangle) +@inline SimpleSymmetric{N}(lowertriangle::SVector{L, T}) where {N, T, L} = + SimpleSymmetric{N, T, L}(lowertriangle) + +@inline function SimpleSymmetric(lowertriangle::SVector{L, T}) where {T, L} + N = triangularroot(Val(L)) + SimpleSymmetric{N, T, L}(lowertriangle) +end + +@generated function SimpleSymmetric{N, T, L}(a::Tuple) where {N, T, L} + _check_simple_symmetric_parameters(Val(N), Val(L)) + expr = Vector{Expr}(undef, L) + i = 0 + for col in 1:N, row in col:N + index = N * (col - 1) + row + expr[i += 1] = :(a[$index]) + end + quote + Base.@_inline_meta + @inbounds return SimpleSymmetric{N, T, L}( + SVector{L, T}(tuple($(expr...))), + ) + end +end + +@inline function SimpleSymmetric{N, T}(a::Tuple) where {N, T} + L = triangularnumber(N) + SimpleSymmetric{N, T, L}(a) +end + +@inline (::Type{SSC})(a::SimpleSymmetric) where {SSC <: SimpleSymmetric} = + SSC(a.lowertriangle) + +@inline (::Type{SSC})(a::AbstractVector) where {SSC <: SimpleSymmetric} = + SSC(convert(lowertriangletype(SSC), a)) + +# disambiguation +@inline (::Type{SSC})( + a::StaticArray{<:Tuple, <:Any, 1}, +) where {SSC <: SimpleSymmetric} = SSC(convert(SVector, a)) + +@generated function _compact_indices(::Val{N}) where {N} + # Returns a Tuple{Pair{Int, Bool}} I such that for linear index i, + # * I[i][1] is the index into the lowertriangle field of a SimpleSymmetric{N}; + # * I[i][2] is true iff i is an index into the lower triangle of an N × N matrix. + indexmat = Matrix{Int}(undef, N, N) + i = 0 + for col in 1:N, row in 1:N + indexmat[row, col] = if row >= col + i += 1 + else + indexmat[col, row] + end + end + quote + Base.@_inline_meta + return $(tuple(indexmat...)) + end +end + +Base.@propagate_inbounds function Base.getindex( + a::SimpleSymmetric{N}, + i::Int, +) where {N} + I = _compact_indices(Val(N)) + j = I[i] + @inbounds value = a.lowertriangle[j] + return value +end + +Base.@propagate_inbounds function Base.setindex( + a::SimpleSymmetric{N, T, L}, + x, + i::Int, +) where {N, T, L} + I = _compact_indices(Val(N)) + j = I[i] + value = x + return SimpleSymmetric{N}(setindex(a.lowertriangle, value, j)) +end + +# needed because it is used in convert.jl and the generic fallback is slow +@generated function Base.Tuple(a::SimpleSymmetric{N}) where {N} + exprs = [:(a[$i]) for i in 1:(N^2)] + quote + Base.@_inline_meta + tuple($(exprs...)) + end +end + +LinearAlgebra.issymmetric(a::SimpleSymmetric) = true + +# TODO: factorize? + +@inline Base.:(==)(a::SimpleSymmetric, b::SimpleSymmetric) = + a.lowertriangle == b.lowertriangle + +@inline function Base.map(f, a1::SimpleSymmetric, as::AbstractArray...) + _map(f, a1, as...) +end + +@generated function _map(f, a::SimpleSymmetric...) + S = Size(a[1]) + N = S[1] + L = triangularnumber(N) + exprs = Vector{Expr}(undef, L) + for i in 1:L + tmp = [:(a[$j].lowertriangle[$i]) for j in 1:length(a)] + exprs[i] = :(f($(tmp...))) + end + return quote + Base.@_inline_meta + same_size(a...) + @inbounds return SimpleSymmetric(SVector(tuple($(exprs...)))) + end +end + +@inline Base.:*(a::Number, b::SimpleSymmetric) = + SimpleSymmetric(a * b.lowertriangle) +@inline Base.:*(a::SimpleSymmetric, b::Number) = + SimpleSymmetric(a.lowertriangle * b) + +@inline Base.:/(a::SimpleSymmetric, b::Number) = + SimpleSymmetric(a.lowertriangle / b) +@inline Base.:\(a::Number, b::SimpleSymmetric) = + SimpleSymmetric(a \ b.lowertriangle) + +@generated function _plus_uniform( + ::Size{S}, + a::SimpleSymmetric{N, T, L}, + λ, +) where {S, N, T, L} + @assert S[1] == N + @assert S[2] == N + exprs = Vector{Expr}(undef, L) + i = 0 + for col in 1:N, row in col:N + i += 1 + exprs[i] = + row == col ? :(a.lowertriangle[$i] + λ) : :(a.lowertriangle[$i]) + end + return quote + Base.@_inline_meta + R = promote_type(eltype(a), typeof(λ)) + SimpleSymmetric{N, R, L}(SVector{L, R}(tuple($(exprs...)))) + end +end + +LinearAlgebra.transpose(a::SimpleSymmetric) = a + +@generated function _one( + ::Size{S}, + ::Type{SSC}, +) where {S, SSC <: SimpleSymmetric} + N = S[1] + L = triangularnumber(N) + T = eltype(SSC) + if T == Any + T = Float64 + end + exprs = Vector{Expr}(undef, L) + i = 0 + for col in 1:N, row in col:N + exprs[i += 1] = row == col ? :(one($T)) : :(zero($T)) + end + quote + Base.@_inline_meta + return SimpleSymmetric(SVector(tuple($(exprs...)))) + end +end + +@inline _scalar_matrix( + s::Size{S}, + t::Type{SSC}, +) where {S, SSC <: SimpleSymmetric} = _one(s, t) + +# _fill covers fill, zeros, and ones: +@generated function _fill( + val, + ::Size{s}, + ::Type{SSC}, +) where {s, SSC <: SimpleSymmetric} + N = s[1] + L = triangularnumber(N) + v = [:val for i in 1:L] + return quote + Base.@_inline_meta + $SSC(SVector(tuple($(v...)))) + end +end + +# import Random +# import Random: AbstractRNG +# @generated function _rand( +# randfun, +# rng::AbstractRNG, +# ::Type{SSC}, +# ) where {N, SSC <: SimpleSymmetric{N}} +# T = eltype(SSC) +# if T == Any +# T = Float64 +# end +# L = triangularnumber(N) +# v = [:(randfun(rng, $T)) for i in 1:L] +# return quote +# Base.@_inline_meta +# $SSC(SVector(tuple($(v...)))) +# end +# end + +# @inline Random.rand( +# rng::AbstractRNG, +# ::Type{SSC}, +# ) where {SSC <: SimpleSymmetric} = _rand(rand, rng, SSC) +# @inline Random.randn( +# rng::AbstractRNG, +# ::Type{SSC}, +# ) where {SSC <: SimpleSymmetric} = _rand(randn, rng, SSC) +# @inline Random.randexp( +# rng::AbstractRNG, +# ::Type{SSC}, +# ) where {SSC <: SimpleSymmetric} = _rand(randexp, rng, SSC) diff --git a/src/Grids/finitedifference.jl b/src/Grids/finitedifference.jl index aa7f739551..12ab69c447 100644 --- a/src/Grids/finitedifference.jl +++ b/src/Grids/finitedifference.jl @@ -79,7 +79,8 @@ function fd_geometry_data( ) where {FT, periodic} CT = Geometry.ZPoint{FT} AIdx = (3,) - LG = Geometry.LocalGeometry{AIdx, CT, FT, SMatrix{1, 1, FT, 1}} + S = SMatrix{1, 1, FT, 1} + LG = Geometry.LocalGeometry{AIdx, CT, Geometry.tail_params(S)...} (Ni, Nj, Nk, Nv, Nh) = size(face_coordinates) Nv_face = Nv - periodic Nv_cent = Nv - 1 diff --git a/src/Grids/spectralelement.jl b/src/Grids/spectralelement.jl index 47681d3a6f..a2211bfde6 100644 --- a/src/Grids/spectralelement.jl +++ b/src/Grids/spectralelement.jl @@ -50,7 +50,8 @@ function _SpectralElementGrid1D( nelements = Topologies.nlocalelems(topology) Nq = Quadratures.degrees_of_freedom(quadrature_style) - LG = Geometry.LocalGeometry{AIdx, CoordType, FT, SMatrix{1, 1, FT, 1}} + S = SMatrix{1, 1, FT, 1} + LG = Geometry.LocalGeometry{AIdx, CoordType, Geometry.tail_params(S)...} local_geometry = DataLayouts.IFH{LG, Nq}(Array{FT}, nelements) quad_points, quad_weights = Quadratures.quadrature_points(FT, quadrature_style) @@ -218,7 +219,8 @@ function _SpectralElementGrid2D( high_order_quadrature_style = Quadratures.GLL{Nq * 2}() high_order_Nq = Quadratures.degrees_of_freedom(high_order_quadrature_style) - LG = Geometry.LocalGeometry{AIdx, CoordType2D, FT, SMatrix{2, 2, FT, 4}} + S = SMatrix{2, 2, FT, 4} + LG = Geometry.LocalGeometry{AIdx, CoordType2D, Geometry.tail_params(S)...} local_geometry = DataLayouts.IJFH{LG, Nq}(Array{FT}, nlelems) diff --git a/test/DataLayouts/opt_similar.jl b/test/DataLayouts/opt_similar.jl index f69383a34c..8dd43af565 100644 --- a/test/DataLayouts/opt_similar.jl +++ b/test/DataLayouts/opt_similar.jl @@ -16,7 +16,8 @@ function test_similar!(data) FT = eltype(parent(data)) CT = Geometry.ZPoint{FT} AIdx = (3,) - LG = Geometry.LocalGeometry{AIdx, CT, FT, SMatrix{1, 1, FT, 1}} + S = SMatrix{1, 1, FT, 1} + LG = Geometry.LocalGeometry{AIdx, CT, Geometry.tail_params(S)...} (_, _, _, Nv, _) = size(data) similar(data, LG, Val(Nv)) @test_opt similar(data, LG, Val(Nv)) diff --git a/test/Geometry/axistensor_conversion_benchmarks.jl b/test/Geometry/axistensor_conversion_benchmarks.jl index 9188561363..716f681042 100644 --- a/test/Geometry/axistensor_conversion_benchmarks.jl +++ b/test/Geometry/axistensor_conversion_benchmarks.jl @@ -1,3 +1,7 @@ +#= +julia --project +using Revise; include(joinpath("test", "Geometry", "axistensor_conversion_benchmarks.jl")) +=# using Test, StaticArrays #! format: off import Random, BenchmarkTools, StatsBase, @@ -109,6 +113,11 @@ function benchmark_func(args, key, f, flops, ::Type{FT}; print_method_info) wher print("Time (opt, ref): ($(opt.t_pretty), $(ref.t_pretty)). Key: $key_str\n") # end end + correctness = compare(components(opt.result), components(ref.result)) # test correctness + # @show correctness + @show components(opt.result) + @show components(ref.result) + @show correctness bm = (; opt, ref, @@ -117,7 +126,7 @@ function benchmark_func(args, key, f, flops, ::Type{FT}; print_method_info) wher flops, # current flops computed_flops, reduced_flops, - correctness = compare(opt.result, ref.result), # test correctness + correctness, # test correctness perf_pass = (opt.time - ref.time)/ref.time*100 < -100, # test performance ) return bm @@ -145,6 +154,7 @@ components(x::T) where {T <: Real} = x components(x) = Geometry.components(x) compare(x::T, y::T) where {T<: Real} = x ≈ y || (x < eps(T)/100 && y < eps(T)/100) compare(x::T, y::T) where {T <: SMatrix} = all(compare.(x, y)) +compare(x::T, y::T) where {T <: Geometry.SimpleSymmetric} = all(compare.(x.lowertriangle, y.lowertriangle)) compare(x::T, y::T) where {T <: SVector} = all(compare.(x, y)) compare(x::T, y::T) where {T <: AxisTensor} = compare(components(x), components(y)) @@ -164,6 +174,9 @@ function test_optimized_functions(::Type{FT}; print_method_info=false) where {FT end for key in keys(benchmarks) + if !(benchmarks[key].correctness) + @show key + end @test benchmarks[key].correctness # test correctness @test benchmarks[key].Δflops ≤ 0 # Don't regress # @test_broken benchmarks[key].Δflops < 0 # Error on improvements. TODO: fix, this is somehow flakey diff --git a/test/Geometry/unit_simple_symmetric.jl b/test/Geometry/unit_simple_symmetric.jl new file mode 100644 index 0000000000..629083b335 --- /dev/null +++ b/test/Geometry/unit_simple_symmetric.jl @@ -0,0 +1,39 @@ +#= +julia --project=.buildkite +using Revise; include(joinpath("test", "Geometry", "unit_simple_symmetric.jl")) +=# +using Test +using StaticArrays +using ClimaCore.Geometry: SimpleSymmetric +import ClimaCore.Geometry +using JET +simple_symmetric(A::Matrix) = SimpleSymmetric(SMatrix{size(A)..., eltype(A)}(A)) + +@testset "SimpleSymmetric" begin + A = simple_symmetric([1 2; 3 4]) # pass in non-symmetric matrix + @test A[1, 1] == 1 + @test A[1, 2] == 2 + @test A[2, 1] == 2 # Do not access below diagonal + @test A[2, 2] == 4 + + B₀ = @SMatrix [1 2; 2 4] + A = SimpleSymmetric(B₀) + C = A * B₀ + @test C isa SMatrix + @test C == B₀ * B₀ + + A = @SMatrix [1 2; 2 4] + @test SimpleSymmetric(A) * 2 === SimpleSymmetric(A * 2) + + A = @SMatrix [1 2; 2 4] + @test SimpleSymmetric(A) / 2 === SimpleSymmetric(A / 2) + @test_opt SimpleSymmetric(A) + @test Geometry.tail_params(typeof(@SMatrix Float32[1 2; 2 4])) == + (Float32, SMatrix{2, 2, Float32, 4}, 2, 3) +end + +@testset "sizs" begin + for N in (1, 2, 3, 5, 8, 10) + simple_symmetric(rand(N, N)) # pass in non-symmetric matrix + end +end