Compute inverse of Cholesky with a single solve (#1326)

Alternative to #1258 that instead of averaging out the errors of the two
solves computes the inverse with a single solve and a symmetric matrix
multiplication.

Checks of symmetry in the tests pass locally, and the example in the
tests is even slightly faster:

## This PR

```julia
julia> using StaticArrays, LinearAlgebra, Chairmarks

julia> @be cholesky((x -> x * x' + I)(@smatrix(randn(3,3)))) inv
Benchmark: 3702 samples with 4741 evaluations
 min    4.561 ns
 median 5.018 ns
 mean   5.419 ns
 max    11.680 ns
```

## master

```julia
julia> using StaticArrays, LinearAlgebra, Chairmarks

julia> @be cholesky((x -> x * x' + I)(@smatrix(randn(3,3)))) inv
Benchmark: 3756 samples with 3276 evaluations
 min    6.983 ns
 median 7.326 ns
 mean   7.712 ns
 max    53.406 ns
```

---------

Co-authored-by: Mateusz Baran <mateuszbaran89@gmail.com>

Author: devmotion

Project.toml

@@ -1,6 +1,6 @@
 name = "StaticArrays"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.15"
+version = "1.9.16"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/cholesky.jl

@@ -70,8 +70,9 @@ end
 LinearAlgebra.hermitian_type(::Type{SA}) where {T, S, SA<:SArray{S,T}} = Hermitian{T,SA}
 
 function inv(A::Cholesky{T,<:StaticMatrix{N,N,T}}) where {N,T}
-    B = A.U \ (A.U' \ SDiagonal{N}(I))
-    return (B .+ B') ./ T(2)
+    # We can't use M = inv(A.U) here since it returns an `UpperTriangular` of a `Matrix`
+    M = UpperTriangular(A.U \ SDiagonal{N}(I))
+    return M * M'
 end
 
 function Base.:\(A::Cholesky{T,<:StaticMatrix{N,N,T}}, B::StaticVecOrMatLike) where {N,T}