I moved this discussion to https://discourse.julialang.org where it is easier to format the code chunks.
On Monday, October 31, 2016 at 2:04:02 PM UTC-5, Douglas Bates wrote: > > I am encountering an unexpected amount of storage allocation in the > cfactor!{A::HBlkDiag) method in the MixedModels package. See > https://github.com/dmbates/MixedModels.jl/blob/master/src/cfactor.jl for > the code. > > An HBlkDiag matrix is a homogeneous block diagonal matrix where > "homogeneous" refers to the fact that all the diagonal blocks are square > and of the same size. Because of this homogeneity the blocks can be stored > in an r by r by k array where r is the size of each of the square block and > k is the number of such blocks. > > On entry to this method the blocks are symmetric, positive semi-definite. > I want to overwrite the upper triangle of each of these blocks with its > Cholesky factor. I call LAPACK.potrf! directly because I don't want > cholfact! to throw a non-positive-definite error. The strange thing to me > is that when I monitor the storage allocation, I get a huge amount of > storage being allocated in the line with that call. This may be because > LAPACK.potrf! returns a tuple of the original matrix and an Int (the info > code) but I'm not sure. > > To see an example of this unusual amount of allocation try the following > code with julia --track-allocation=user > > using Feather, MixedModels > cd(Pkg.dir("MixedModels", "test", "data")) > sleepstudy = Feather.read("sleepstudy.feather", nullable=false) > fm1 = fit!(lmm(Reaction ~ 1 + Days + (1 + Days | Subject), sleepstudy)) > Profile.clear_malloc_data() > devs, vars, betas, thetas = bootstrap(10_000, fm1) > > I get > > - function cfactor!{T}(A::HBlkDiag{T}) > - Aa = A.arr > 0 r, s, t = size(Aa) > 0 if r ≠ s > 0 throw(ArgumentError("HBlkDiag matrix A must be square")) > - end > 94428000 scm = Array(T, (r, r)) > 0 for k in 1 : t # FIXME: Lots of allocations in this loop > 0 for j in 1 : r, i in 1 : j > 0 scm[i, j] = Aa[i, j, k] > - end > 566568000 LAPACK.potrf!('U', scm) > 0 for j in 1 : r, i in 1 : j > 0 Aa[i, j, k] = scm[i, j] > - end > - end > 10492000 UpperTriangular(A) > - end > > In this case the HBlkDiag matrix being decomposed has is 36 by 36 > consisting of 18 2 by 2 diagonal blocks. scm is a scratch 2 by 2 > matrixthat is overwritten in sequence by the upper triangle of each of the > original 2 by 2 blocks and passed to LAPACK.potrf! >