Wrong results of matrix inverse multiplication with Float16 static arrays

[heliosdrm]

This is my configuration (using StaticArrays v1.9.15):

julia> versioninfo()
Julia Version 1.10.10
Commit 95f30e51f4 (2025-06-27 09:51 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 24 × Intel(R) Core(TM) Ultra 9 275HX
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, goldmont)
Threads: 24 default, 0 interactive, 12 GC (on 24 virtual cores)
Environment:
  JULIA_DEPOT_PATH = D:\.julia
  JULIA_EDITOR = code
  JULIA_VSCODE_REPL = 1

And this is a working example:

using StaticArrays

A = [0.0517 0.05106 0.05518;
     0.8870 0.88000 0.88040;
     1.0    1.0     1.0]

b = [0.04922, 0.889, 1.0]

x = A \ b     # correct

## Using StaticArrays:
Astatic = SMatrix{3,3}(A)
bstatic = SVector{3}(b)
xstatic = Astatic \ bstatic     # correct: xstatic ≈ x

## Changing to Float32 everything is still right:
A32 = Matrix{Float32}(A)
b32 = Vector{Float32}(b)
x32 = A32 \ b32     # correct

A32static = SMatrix{3,3,Float32}(A)
b32static = SVector{3,Float32}(b)
x32static = A32static \ b32static     # correct

## But things go wrong with Float16 and StaticArrays:
A16 = Matrix{Float16}(A)
b16 = Vector{Float16}(b)
x16 = A16 \ b16     # correct

A16static = SMatrix{3,3,Float16}(A)
b16static = SVector{3,Float16}(b)
x16static = A16static \ b16static     # WRONG: x ≈ [1.32, 0.33, -0.65]; x16static ≈ [1.33, 0.6665, 0.15]

The discrepancy of the last result is way beyond rounding errors.

[mateuszbaran]

This is basically a duplicate of #1285, though with Float16 under/overflow is even easier to get. As a workaround, you can use

x16static = lu(A16static) \ b16static     # correct again

By default \ doesn't do the LU decomposition because it's very costly (it would introduce roughly x10 slowdown for Float64), and sometimes the decomposition-free version is accurate enough.