# MArray or SizedArray in StaticArrays package

I am a bit confused about the distinction between a mutable array and a sized array in the `StaticArrays` package. My application is block diagonal matrices with all of the blocks having the same, usually small, size, `k × k`.

One block diagonal matrix, `A`, is symmetric and positive semi-definite. Given a lower-triangular matrix, `T`, of size `k×k`, I fill out the diagonal blocks of a similar, lower-triangular matrix, `L` according to

``````for i in 1:eachindex(A.diag)
_chol(copy!(L.diag[i], T'A.diag[i]*T + I), Val{:L})
end
``````

where `_chol` is some suitable in-place lower Cholesky factorization method (the matrix multiplications and addition of `I` will be done in-place)

`A` could be `Diagonal{SMatrix{k,k,T}}` but the diagonal blocks of `L` must be mutable. Should I use `Diagonal{MMatrix{k,k,T})` or `Diagonal{SizedArray{(k,k),T,k,k}}` or something else?

Hi @dmbates,

Two things:

1. The difference between `SizedArray` and `MArray` is just the internal representation.
• `SizedArray` is a wrapper of `Array` and is useful whenever you already have an `Array` and you happen to know its size, and want to use the accelerated methods in `StaticArrays`. It is fast to construct given an `Array` but slightly slower from other types. It is also better if the element type is not an `isbits` immutable.

• `MArray` is the mutable form `SArray` using a `Tuple` internally. They both have the same functionality for `isbits` immutable element types.

1. Mutability of the static matrices in this case might not be optimal. You already have a mutable `Diagonal` (backed by `Vector`), so `Diagonal{SMatrix{k,k,T,k*k}}` might actually be ideal (I use this type of structure all the time).

This is an example of something which shouldn’t allocate memory (sorry if this doesn’t work perfectly - I don’t have access to the REPL right now). It does immutable operations on your small `k` x `k` chunks and saves them back to the mutable `Diagonal`.

``````v = [(tmp = rand(SMatrix{k,k}); tmp*tmp') for i = 1:N]
d = Diagonal(v)
map!(x -> chol(Symmetric(x)), d, d)
``````

It gets tricky to reuse memory and change the element type (e.g. go from `Symmetric` to `LowerTriangular`), but perhaps these wrappers can go outside the `Diagonal` or not used at all.

EDIT: Sorry `map!` doesn’t seem to work well with `Diagonal`… you can just mutate the diagonal entries manually, though.

1 Like

(what is `k`? I think `chol` might only be optimized in `StaticArrays` up to 3x3, but I don’t remember at the moment).

@dmbates - how did you go with this?

Sorry for the delay in responding. New equipment arrived and I got distracted.

I have spent some time working with MMatrix types in this application and reached the conclusion that what I want to do is perhaps too complicated at present, mostly because I want in-place operations.

In most applications `k` is small (single digit) but the number of blocks on the diagonal can be very large (thousands). I don’t want to keep allocating storage in the multiplication operations. Calling out to BLAS.trmm! to evaluate these products for small matrices is not the best but each of these calls takes a few milliseconds at most and I can live with that.

Thanks for responding to my question.

OK.

Just another idea is to keep some small, `k` * `k` arrays about as temporaries, where you can copy data from the block-diagonal structure and back again.

You can even make these into (mutable) `SizedArray`. You will most likely get better performance from a full `SizedArray` multiplication then a triangular multiplication using BLAS when `k` is small (use `A_mul_B!` to avoid allocations (but don’t alias the output with the input)).