How to preallocate a static array for use with `lu!`

General Usage
1

If you call lu on an SMatrix, you get a static LU factorization—cool:

using StaticArrays
using LinearAlgebra

A = SMatrix{3, 3}(rand(3, 3))
lu(A)
# StaticArrays.LU
# L factor:
# 3×3 LowerTriangular{Float64, SMatrix{3, 3, Float64, 9}}:
#  1.0        ⋅         ⋅ 
#  0.706739  1.0        ⋅ 
#  0.648963  0.100059  1.0
# U factor:
# 3×3 UpperTriangular{Float64, SMatrix{3, 3, Float64, 9}}:
#  0.807758  0.146236  0.681894
#   ⋅        0.76478   0.338519
#   ⋅         ⋅        0.412485

But sometimes, you want to preallocate space for the L and U factors. In this case, I would like L to be a LowerTriangular based on an MMatrix (needs to be mutable, obviously). How can I initialize an empty such LU factorization to be used as follows?

lu!(what_goes_here, A)

My attempt:

my_LU = StaticArrays.LU(
    LowerTriangular{Float64, zeros(MMatrix{3, 3, Float64})},
    UpperTriangular{Float64, zeros(MMatrix{3, 3, Float64})}
)
lu!(my_LU, A)
# TypeError: in Type, in parameter, expected Type, got a value of type MMatrix{3, 3, Float64, 9}

(And why isn’t MMatrix a subtype of AbstractMatrix?)

2

Preallocation doesn’t make sense for SMatrix — there is no “allocation” cost in the first place. So what is your motivation?

3

For small matrices, won’t my_LU \ b benefit from having my_LU as an SArray for the same reason that A \ b does? (Or is the latter assumption also incorrect?)

4

lu(smatrix) already returns L and U factors that are backed by SMatrix. What more do you want?

(Note that StaticArrays already has optimized methods for triangular solves with SMatrix storage.)

5

2 posts were split to a new topic: How to preallocate LU factorization