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?
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.
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.
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.
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)).